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 parcalc.parcalc.ClusterAims(nodes=1, ppn=8, **kwargs)

Encapsulating Aims calculator for the cluster environment.

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.

clean()

Method which cleans up after a calculation.

The default files generated by Vasp will be deleted IF this method is called.

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 be 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).

class parcalc.parcalc.RemoteCalculator(restart=None, ignore_bad_restart_file=False, label=None, atoms=None, calc=None, block=False, **kwargs)

Remote calculator based on ASE calculator class. This class is only involved with the machanics of remotly executing the software and transporting the data. The calculation is delegated to the actual calculator class.

classmethod ParallelCalculate(syslst, properties=['energy'], system_changes=['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms'])

Run a series of calculations in parallel using (implicitely) some remote machine/cluster. The function returns the list of systems ready for the extraction of calculated properties.

read_results()

Read energy, forces, … from output file(s).

run_calculation(atoms=None, properties=['energy'], system_changes=['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms'])

Internal calculation executor. We cannot use FileIOCalculator directly since we need to support remote execution.

This calculator is different from others. It prepares the directory, launches the remote process and raises the exception to signal that we need to come back for results when the job is finished.

write_input(atoms=None, properties=['energy'], system_changes=['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms'])

Write input file(s).

parcalc.parcalc.work_dir(path)

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

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 (up to 4.0). 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 with the name of bravais lattices 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

---

elastic.elastic.get_BM_EOS(cryst, systems)

Calculate Birch-Murnaghan Equation of State for the crystal.

The B-M equation of state is defined by:

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

It’s coefficients are estimated using n single-point structures ganerated from the crystal (cryst) by the scan_volumes function between two 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 raised. The data from the calculation and fit is stored in the bm_eos and pv members of cryst for future reference. You have to provide properly optimized structures in cryst and systems list.

Parameters:

• cryst – Atoms object, basic structure
• systems – A list of calculated structures

Returns:

tuple of EOS parameters \(V_0, B_0, B_0'\).

elastic.elastic.get_bulk_modulus(cryst)

Calculate bulk modulus using the Birch-Murnaghan equation of state.

The EOS must be previously calculated by get_BM_EOS routine. The returned bulk modulus is a \(B_0\) coefficient of the B-M EOS. The units of the result are defined by ASE. To get the result in any particular units (e.g. GPa) you need to divide it by ase.units.<unit name>:

get_bulk_modulus(cryst)/ase.units.GPa

Parameters:cryst – ASE Atoms object Returns:float, bulk modulus \(B_0\) in ASE units.

elastic.elastic.get_cart_deformed_cell(base_cryst, axis=0, size=1)

Return the cell deformed along one of the cartesian directions

Creates new deformed structure. The deformation is based on the base structure and is performed along single axis. The axis is specified as follows: 0,1,2 = x,y,z ; sheers: 3,4,5 = yz, xz, xy. The size of the deformation is in percent and degrees, respectively.

Parameters:

• base_cryst – structure to be deformed
• axis – direction of deformation
• size – size of the deformation

Returns:

new, deformed structure

elastic.elastic.get_cij_order(cryst)

Give order of of elastic constants for the structure

Parameters:cryst – ASE Atoms object Returns:Order of elastic constants as a tuple of strings: C_ij

elastic.elastic.get_deformed_cell(base_cryst, 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).

elastic.elastic.get_elastic_tensor(cryst, systems)

Calculate elastic tensor of the crystal.

The elastic tensor is calculated from the stress-strain relation and derived by fitting this relation to the set of linear equations build from the symmetry of the crystal and strains and stresses of the set of elementary deformations of the unit cell.

It is assumed that the crystal is converged and optimized under intended pressure/stress. The geometry and stress on the cryst is taken as the reference point. No additional optimization will be run. Structures in cryst and systems list must have calculated stresses. The function returns tuple of \(C_{ij}\) elastic tensor, raw Birch coefficients \(B_{ij}\) and fitting results: residuals, solution rank, singular values returned by numpy.linalg.lstsq.

Parameters:

• cryst – Atoms object, basic structure
• systems – list of Atoms object with calculated deformed structures

Returns:

tuple(\(C_{ij}\) float vector, tuple(\(B_{ij}\) float vector, residuals, solution rank, singular values))

elastic.elastic.get_elementary_deformations(cryst, n=5, d=2)

Generate elementary deformations for elastic tensor calculation.

The deformations are created based on the symmetry of the crystal and are limited to the non-equivalet axes of the crystal.

Parameters:

• cryst – Atoms object, basic structure
• n – integer, number of deformations per non-equivalent axis
• d – float, size of the maximum deformation in percent and degrees

Returns:

list of deformed structures

elastic.elastic.get_lattice_type(cryst)

Find the symmetry of the crystal using spglib symmetry finder.

Derive name of the space group and its number extracted from the result. Based on the group number identify also the lattice type and the Bravais lattice of the crystal. The lattice type numbers are (the numbering starts from 1):

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

Parameters:cryst – ASE Atoms object Returns:tuple (lattice type number (1-7), lattice name, space group name, space group number)

elastic.elastic.get_pressure(s)

Return external isotropic (hydrostatic) pressure in ASE units.

If the pressure is positive the system is under external pressure. This is a convenience function to convert output of get_stress function into external pressure.

Parameters:cryst – stress tensor in Voight (vector) notation as returned by the get_stress() method. Returns:float, external hydrostatic pressure in ASE units.

elastic.elastic.get_strain(cryst, refcell=None)

Calculate strain tensor in the Voight notation

Computes 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.

Parameters:

• cryst – deformed structure
• refcell – reference, undeformed structure

Returns:

6-vector of strain tensor in the Voight notation

elastic.elastic.get_vecang_cell(cryst, uc=None)

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

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}\]

Parameters:u – vector of deformations: [ \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\) ] Returns:Symmetry defined stress-strain equation matrix

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}\]

Parameters:u – vector of deformations: [ \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\) ] Returns:Symmetry defined stress-strain equation matrix

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}\]

Parameters:u – vector of deformations: [ \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\) ] Returns:Symmetry defined stress-strain equation matrix

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}\]

Parameters:u – vector of deformations: [ \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\) ] Returns:Symmetry defined stress-strain equation matrix

elastic.elastic.scan_pressures(cryst, lo, hi, n=5, eos=None)

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

elastic.elastic.scan_volumes(cryst, lo=0.98, hi=1.02, n=5, scale_volumes=True)

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==False the scalling is applied to lattice vectors instead of volumes.

Parameters:

• lo – lower bound of the V/V_0 in the scan
• hi – upper bound of the V/V_0 in the scan
• n – number of volume sample points
• scale_volumes – If True scale the unit cell volume or, if False, scale the length of lattice axes.

Returns:

a list of deformed systems

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}\]

Parameters:u – vector of deformations: [ \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\) ] Returns:Symmetry defined stress-strain equation matrix

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}\]

Parameters:u – vector of deformations: [ \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\) ] Returns:Symmetry defined stress-strain equation matrix

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}\]

Parameters:u – vector of deformations: [ \(u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}\) ] Returns:Symmetry defined stress-strain equation matrix