Implementation

Elastic is implemented as an extension module to ASE system

The Elastic package provides, basically, one main python module and one auxiliary module (Parallel Calculator Module) which can be useful outside of the scope of the main code. The Parallel Calculator Module is not distributed separately but can be just placed by itself somewhere in your python path and used with any part of the ASE. I hope it will be incorporated in the main project sometime in the future.

Modules

Parallel Calculator Module

Parallel calculator module is an extension of the standard ASE calculator working in the parallel cluster environment. It is very useful in all situations where you need to run several, independent calculations and you have a large cluster of machines at your disposal (probably with some queuing system).

This implementation uses VASP but the code can be easily adapted for use with other ASE calculators with minor changes. The final goal is to provide a universal module for parallel calculator execution in the cluster environment.

The SIESTA code by Georgios Tritsaris <gtritsaris@seas.harvard.edu> Not fully tested after merge.

Class description

class parcalc.parcalc.ClusterSiesta(nodes=1, ppn=8, **kwargs)

Siesta calculator. Not fully tested by me - so this should be considered beta quality. Nevertheless it is based on working implementation

class parcalc.parcalc.ClusterVasp(nodes=1, ppn=8, block=True, ncl=False, **kwargs)

Adaptation of VASP calculator to the cluster environment where you often have to make some preparations before job submission. You can easily adapt this class to your particular environment. It is also easy to use this as a template for other type of calculator.

calc_finished()

Check if the lockfile is in the calculation directory. It is removed by the script at the end regardless of the success of the calculation. This is totally tied to implementation and you need to implement your own scheme!

calculate(atoms)

Blocking/Non-blocking calculate method

If we are in blocking mode we just run, wait for the job to end and read in the results. Easy ...

The non-blocking mode is a little tricky. We need to start the job and guard against it reading back possible old data from the directory - the queuing system may not even started the job when we get control back from the starting script. Thus anything we read after invocation is potentially garbage - even if it is a converged calculation data.

We handle it by custom run function above which raises an exception after submitting the job. This skips post-run processing in the calculator, preserves the state of the data and signals here that we need to wait for results.

prepare_calc_dir()

Prepare the calculation directory for VASP execution. This needs to be re-implemented for each local setup. The following code reflects just my particular setup.

run()

Blocking/Non-blocing run method. In blocking mode it just runs parent run method. In non-blocking mode it raises the __NonBlockingRunException to bail out of the processing of standard calculate method (or any other method in fact) and signal that the data is not ready to b collected.

parcalc.parcalc.ParCalculate(systems, calc, cleanup=True, block=True, prefix='Calc_')

Run calculators in parallel for all systems. Calculators are executed in isolated processes and directories. The resulting objects are returned in the list (one per input system).

parcalc.parcalc.work_dir(*args, **kwds)

Context menager for executing commands in some working directory. Returns to the previous wd when finished.

Usage: >>> with work_dir(path): ... subprocess.call(‘git status’)

Elastic Module

This module depends on Parallel Calculator Module for parallelisation of independent calculations.

Elastic is a module for calculation of \(C_{ij}\) components of elastic tensor from the strain-stress relation.

The strain components here are ordered in standard way which is different to ordering in previous versions of the code.

The ordering is: \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\).

The general ordering of \(C_{ij}\) components is (except for triclinic symmetry and taking into account customary names of constants - e.g. \(C_{16} \rightarrow C_{14}\)):

\[C_{11}, C_{22}, C_{33}, C_{12}, C_{13}, C_{23}, C_{44}, C_{55}, C_{66}, C_{16}, C_{26}, C_{36}, C_{45}\]

The functions outside of the Crystal class define the symmetry of the \(C_{ij}\) matrix. The matrix is N columns by 6 rows where the columns corespond to independent elastic constants of the given crystal, while the rows corespond to the canonical deformations of a crystal. The elements are the second partial derivatives of the free energy formula for the crystal written down as a quadratic form of the deformations with respect to elastic constant and deformation.

Note: The elements for deformations \(u_{xy}, u_{xz}, u_{yz}\) have to be divided by 2 to properly match the usual definition of elastic constants.

See: [LL] L.D. Landau, E.M. Lifszyc, “Theory of elasticity”

There is some usefull summary also at: ScienceWorld

Class description

class elastic.elastic.Crystal(symbols=None, positions=None, numbers=None, tags=None, momenta=None, masses=None, magmoms=None, charges=None, scaled_positions=None, cell=None, pbc=None, celldisp=None, constraint=None, calculator=None, info=None)

Backward compatibility class. To be removed later.

class elastic.elastic.ElasticCrystal

Mixin extension of standard ASE Atoms class designed to handle specifics of the crystalline materials. This code should, in principle, be folded into the Atoms class in the future. At this moment it is too early to think about it. Additionally there are some aspects of this code which may be difficult to harmonize with the principles of the Atoms class. I am sure it is better, for now to leave this as a separate extension class.

Basically, this class provides set of functions concerned with derivation of elastic properties using “finite deformation approach” (see the documentation for physics background information).

get_BM_EOS(n=5, lo=0.98, hi=1.02, recalc=False, cleanup=True, mode='full', data=None)

Calculate Birch-Murnaghan Equation of State for the crystal:

\[P(V)= \frac{B_0}{B'_0}\left[ \left({\frac{V}{V_0}}\right)^{-B'_0} - 1 \right]\]

using n single-point structures ganerated from the crystal (self) by the scan_volumes method between lo and hi relative volumes. The BM EOS is fitted to the computed points by least squares method. The returned value is a list of fitted parameters: \(V_0, B_0, B_0'\) if the fit succeded. If the fitting fails the RuntimeError(‘Calculation failed’) is reised. The data from the calculation and fit is stored in the bm_eos and pv members for future reference.

Note: For now you have to set up the calculator to properly optimize the structure without changing the volume at each point. There will be a way to specify basic types of the calculator minimization at the later stage.

get_bulk_modulus(n=5, lo=0.98, hi=1.02, recalc=False)

Calculate bulk modulus using the Birch-Murnaghan equation of state data calculated by get_BM_EOS routine (see). The returned bulk modulus is a \(B_0\) coefficient of the B-M EOS. The arguments are the same as in BM EOS function.

get_cart_deformed_cell(axis=0, size=1)

Return the cell (with atoms) deformed along one of the cartesian directions (0,1,2 = x,y,z ; sheers: 3,4,5 = yz, xz, xy) by size percent.

get_deformed_cell(axis=0, size=1)

Return the cell (with atoms) deformed along one cell parameter (0,1,2 = a,b,c ; 3,4,5 = alpha,beta,gamma) by size percent or size degrees (axis/angles).

get_elastic_tensor(n=5, d=2, mode='full', systems=None)

Calculate elastic tensor of the crystal. It is assumed that the crystal is converged and optimized under intended pressure/stress. The geometry and stress at the call point is taken as the reference point. No additional optimization will be run. It is also assumed that the calculator is set to pure IDOF optimization. The size of used finite deformation is passed in d parameter as a percentage relative deformation. The n parameter defines number of deformed structures used in the calculation.

get_lattice_type()

Find the symmetry of the crystal using spglib symmetry finder. Assign to sg_name i sg_nr members name of the space group and its number extracted from the result. Based on the group number identify also the lattice type (assigned to sg_type member) and the Bravais lattice of the crystal (assigned to bravais member). The returned value is the lattice type number. The lattice type numbers are (see also Crystal.ls, the numbering starts from 1):

Triclinic (1), Monoclinic (2), Orthorombic (3), Tetragonal (4) Trigonal (5), Hexagonal (6), Cubic (7)

get_pressure(s=None)

Return external isotropic (hydrostatic) pressure in ASE units. If the pressure is positive the system is under external pressure. This is a convenience function.

get_strain(refcell=None)

Return the strain tensor in the Voight notation as a conventional 6-vector. The calculation is done with respect to the crystal geometry passed in refcell parameter.

get_vecang_cell(uc=None)

Compute A,B,C, alpha,beta,gamma cell params from the unit cell matrix (uc) or self. Angles in radians.

scan_pressures(lo, hi, n=5)

Scan the pressure axis from lo to hi (inclusive) using B-M EOS as the volume predictor. Pressure (lo, hi) in GPa

scan_volumes(lo, hi, n, scale_volumes=False)

Provide set of crystals along volume axis from lo to hi (inclusive). No volume cell optimization is performed. Bounds are specified as fractions (1.10 = 10% increase). If scale_volumes==True the scalling is applied to volumes instead of lattice vectors.

elastic.elastic.hexagonal(u)

The matrix is constructed based on the approach from L&L using auxiliary coordinates: \(\xi=x+iy\), \(\eta=x-iy\). The components are calculated from free energy using formula introduced in Crystal symmetry and elastic matrix derivation with appropriate coordinate changes. The order of constants is as follows:

\[C_{11}, C_{33}, C_{12}, C_{13}, C_{44}\]

elastic.elastic.monoclinic(u)

Monoclinic group, the ordering of constants is:

\[C_{11}, C_{22}, C_{33}, C_{12}, C_{13}, C_{23}, C_{44}, C_{55}, C_{66}, C_{16}, C_{26}, C_{36}, C_{45}\]

elastic.elastic.orthorombic(u)

Equation matrix generation for the orthorombic lattice. The order of constants is as follows:

\[C_{11}, C_{22}, C_{33}, C_{12}, C_{13}, C_{23}, C_{44}, C_{55}, C_{66}\]

elastic.elastic.regular(u)

Equation matrix generation for the regular (cubic) lattice. The order of constants is as follows:

\[C_{11}, C_{12}, C_{44}\]

elastic.elastic.tetragonal(u)

Equation matrix generation for the tetragonal lattice. The order of constants is as follows:

\[C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{14}\]

elastic.elastic.triclinic(u)

Triclinic crystals.

Note: This was never tested on the real case. Beware!

The ordering of constants is:

\[C_{11}, C_{22}, C_{33}, C_{12}, C_{13}, C_{23}, C_{44}, C_{55}, C_{66}, C_{16}, C_{26}, C_{36}, C_{46}, C_{56}, C_{14}, C_{15}, C_{25}, C_{45}\]

elastic.elastic.trigonal(u)

The matrix is constructed based on the approach from L&L using auxiliary coordinates: \(\xi=x+iy\), \(\eta=x-iy\). The components are calculated from free energy using formula introduced in Crystal symmetry and elastic matrix derivation with appropriate coordinate changes. The order of constants is as follows:

\[C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{14}\]