jochym / elastic

A module for ASE for elastic constants calculation.

Elastic

Elastic is a set of python routines for calculation of elastic properties of crystals (elastic constants, equation of state, sound velocities, etc.). It is a third version of the in-house code I have written over few years and is implemented as a extension to the ASE system. The code was a basis for some of my publications and was described briefly in these papers. The code was available to anyone, presented at the Workshop on ab initio Calculations in Geosciences and used by some of my co-workers. This code is a re-implementation of elastic as a module for the ASE.

The on-line documentation is placed on ReadTheDocs or Elastic website. You can obtain the documentation as a PDF file as well.

Installation

The installation is simple if you are using conda package menager:

conda install -c conda-forge elastic

If you need a more elaborate installation instruction for computing environment to work with ASE - I have written just such a document. It is a first in the series of tutorials of this subject and you can find it under nbviewer link above.

The project is open and I welcome patches, ideas and other feedback.

Acknowledgments

In the past the project was partially supported by:

• Department of Computational Material Science, Institute of Nuclear Physics, PAN, Poland
• Department of Engineering, University of Saskatchewan, Canada

elastic's Issues

Birch-Murnaghan equation of state

Thank you so much for sharing your code. I met some issues in the code. I think the equation used in file 'elastic.py' function 'BMEOS' is not Birch-Murnaghan equation.
It is Murnaghan equation: https://en.wikipedia.org/wiki/Murnaghan_equation_of_state
The correct equation is: https://en.wikipedia.org/wiki/Birch%E2%80%93Murnaghan_equation_of_state

ImportError: cannot import name 'Crystal' from 'elastic'

Hi again, I installed the elastic package, and I used the test case in 'tests' folder. I ran it, but it had an error 'ImportError: cannot import name 'Crystal' from 'elastic''
I also didn't find this class in 'elastic.py' code, could you please tell me where this class is?

__version__.txt is missing

During installation via pip this error happens:

Complete output from command python setup.py egg_info:
Traceback (most recent call last):
File "", line 1, in
File "/tmp/pip-build-DPwm2Z/elastic/setup.py", line 5, in
with open("version.txt") as w:
IOError: [Errno 2] No such file or directory: 'version.txt'

Seems that version.txt is missing within the release package version 4.0.28.

Wrong coefficients in `hexagonal`

Apart from #64, I would also like to report wrong coefficients in hexagonal.

Here, we have unnecessary coefficients of 2. This is inconsistent with the Elastic document itself as well as my own derivation (main.pdf).

To demonstrate this is really an issue, I have conducted a test calculation for graphite (graphite_gaussian.zip). In experiments (Ref. 11 in this):

• C11 = 1060 GPa, C12 = 180 GPa, C13 = 15 GPa, C33 = 36.5 GPa, C44 = 4 GPa, C66 = (C11-C12)/2 = 440 GPa

Fo DFT (PBE-D3(BJ)), using VASP IBRION=6, I obtain (without ionic-relaxation contribution)

• C11 = 1075 GPa, C12 = 237 GPa, C13 = -4 GPa, C33 = 36 GPa, C44 = -5 GPa, C66 = (C11-C12)/2 = 419 GPa

Using Elastic with get_elastic_tensor (which calls hexagonal) I get

• C11 = 948 GPa, C12 = 69 GPa C13 = -3 GPa, C33 = 36 GPa, C44 = -5 GPa, C66 = (C11-C12)/2 = 440 GPa

Which is substantially deviated from the IBRION=6 values for C11 and C12, because of the coefficient error in hexagonal.

If we use the one-axis approach (i.e., without hexagonal) (with solving the issue in #64), I get

• C11 = 1090 GPa, C12 = 225 GPa, C13 = -4 GPa, C33 = 36 GPa, C44 = -6 GPa, C66 = (C11-C12)/2 = 432 GPa

Which is consistent with IBRION=6.

If we will fix the coefficients in hexagonal, then get_elastic_tensor gives

• C11 = 1077 GPa, C12 = 240 GPa C13 = -3 GPa, C33 = 36 GPa, C44 = -5 GPa, C66 = (C11-C12)/2 = 419 GPa

which is consistent with IBRION=6.

Wrong matrix for the tetragonal system?

@jochym

I think there maybe something wrong with your code for the tetragonal systems.

def tetragonal(u):
    '''
    Equation matrix generation for the tetragonal lattice.
    The order of constants is as follows:
    .. math::
       C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{66}
    :param u: vector of deformations:
        [ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ]
    :returns: Symmetry defined stress-strain equation matrix
    '''

    uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5]
    return array(
                [[uxx,   0,    uyy,  uzz,      0,      0],
                 [uyy,   0,    uxx,  uzz,      0,      0],
                 [0,     uzz,  0,    uxx+uyy,  0,      0],
                 [0,     0,    0,    0,        2*uyz,  0],
                 [0,     0,    0,    0,        2*uxz,  0],
                 [0,     0,    0,    0,        0,  2*uxy]])

Also,

Elastic/elastic/elastic.py

Line 297 in 8daae37

4: ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_14'),

The last independent elastic constant should be C_66 not C_14.

            4: ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_66'), #error in the original code

Could you please check if my fix is correct?

Thanks,

Adding ASE LAMMPSlib Calculator

Was wondering if you have any guidance on how might be able to proceed with using ASE LAMMPSlib. I tried to simply pass the ASE LAMMPSlib calculator object but got deepcopy errors.

Heres the source to the ASE LAMMPSlib implementation:

https://wiki.fysik.dtu.dk/ase/_modules/ase/calculators/lammpslib.html#LAMMPSlib

Regarding QHA calculation

Dear Dr Jochym

i was going through the tutorial of QHA in your website (http://nbviewer.ipython.org/gist/jochym/334b658cc8b3f6864c23) and found the following warning.

/home/jochym/nipy/local/lib/python2.7/site-packages/qeutil/analyzers.py:293: RuntimeWarning: overflow encountered in square
Cv=ndf_k_B_simps(e_e_dos/sinh(e)2,nu)
/home/jochym/nipy/local/lib/python2.7/site-packages/qeutil/analyzers.py:294: RuntimeWarning: overflow encountered in exp
S=k_B(simps(2_dos_e/(exp(2_e)-1),nu)-simps(dos_(1-exp(-2e)),nu))

I got a similar warning when i was trying

Can you please let me know how to avoid this

With regards

Linu

Elastic constant calculation for zirconia

Dear Dr. Jochym

I was trying to calculate the elastic constant for Zirconia and was facing these this error.

Can you please help in solving this

In [1]:

Import the basic libraries

ASE system

import ase

from ase import Atom, Atoms

from ase import io

from ase.lattice.spacegroup import crystal

Spacegroup/symmetry library

from pyspglib import spglib

iPython utility function

from IPython.core.display import Image

Import the qe-util package

from qeutil import RemoteQE

Access info

import hostj

In [2]:

Import additional library for elastic constants calculations

import elastic

In [3]:

Stup the SiC crystal

Create a cubic crystal with a spacegroup F-43m (216)

a=5.010767

cryst = crystal(['Zr','O','O'],

            [(0, 0, 0), (0.25, 0.25, 0.25),(-0.25,-0.25,-0.25)],

            spacegroup=225,

            cellpar=[a, a, a, 90, 90, 90])

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/ase/lattice/spacegroup/spacegroup.py:380: UserWarning: scaled_positions 1 and 2 are equivalent
'are equivalent'%(kinds[ind], kind))

In [4]:

Check the spacegroup (symmetry) of our creation

spglib.get_spacegroup(cryst)

Out[4]:

'Fm-3m (225)'

In [5]:

Create a Quantum Espresso calculator for our work.

This object encapsulates all parameters of the calculation,

not the system we are investigating.

qe=RemoteQE(label='Zirconia',

           kpts=[4,4,4],

           xc='pz',         # Exchange functional type in the name of the pseudopotentials

           pp_type='hgh',   # Variant of the pseudopotential

           pp_format='UPF', # Format of the pseudopotential files

           ecutwfc=70,

           pseudo_dir='../pspot',

           use_symmetry=True,

           procs=8)         # Use 8 cores for the calculation

Check where the calculation files will reside on the local machine.

print qe.directory

calc/Zirconia.MeFc1O

In [6]:

Assign the calculator to our system

cryst.set_calculator(qe)

In [7]:

Run the calculation to get stress tensor (in Voigt notation, GPa) and pressure (in GPa)

print "Stress tensor (GPa):", cryst.get_stress()/ase.units.GPa

print "External pressure (GPa):", cryst.get_pressure()/ase.units.GPa

Stress tensor (GPa): [ 0.208 0.208 0.208 -0. -0. -0. ]
External pressure (GPa): -0.208

In [8]:

fit=cryst.get_BM_EOS(lo=0.96, # lower bound of the volumes range

                 hi=1.04,    # higher bound of the volumes range

                 n=5)        # number of volume points used in the fit

print "\nA0=%.4f A ; B0=%.1f GPa ; B0'=%.2f " % (fit[0]**(1.0/3), fit[1]/ase.units.GPa, fit[2])

Launching: 1 2 3 4 5
Done: 1 2 3 4 5

A0=5.0085 A ; B0=245.8 GPa ; B0'=3.78

In [9]:

deformations=cryst.get_elastic_tensor(mode='deformations')

In [10]:

elastic.ParCalculate(deformations,qe);

Launching:

---

NotImplementedError Traceback (most recent call last)
in ()
----> 1 elastic.ParCalculate(deformations,qe);

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/elastic/init.pyc in ParCalculate(systems, calc)
240 for s in systems:
241 s.set_calculator(calc.copy())
--> 242 calc.ParallelCalculate(systems,properties=['stress'])
243 return systems
244

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/qeutil/init.pyc in ParallelCalculate(cls, syslst, properties, system_changes)
483 try :
484 s.calc.block=False
--> 485 s.calc.calculate(atoms=s,properties=properties,system_changes=system_changes)
486 except CalcNotReadyError:
487 s.calc.block=True

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/qeutil/init.pyc in calculate(self, atoms, properties, system_changes)
165
166 self.command=self.build_command(properties,self.parameters)
--> 167 self.run_calculation(atoms, properties, system_changes)
168
169 # if {'energy','stress'} & set(properties) :

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/qeutil/init.pyc in run_calculation(self, atoms, properties, system_changes)
403 '''
404 Calculator.calculate(self, atoms, properties, system_changes)
--> 405 self.write_input(self.atoms, properties, system_changes)
406 if self.command is None:
407 raise RuntimeError('Please configure RemoteQE calculator!')

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/qeutil/init.pyc in write_input(self, atoms, properties, system_changes)
363 def write_input(self, atoms=None, properties=['energy'], system_changes=all_changes):
364 '''Write input file(s).'''
--> 365 QuantumEspresso.write_input(self, atoms, properties, system_changes)
366 self.write_pbs_in(properties)
367 subprocess.call(self.copy_out_cmd % {

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/qeutil/init.pyc in write_input(self, atoms, properties, system_changes)
196
197 if set(['energy','stress','forces','bands','edos']) & set(properties) :
--> 198 write_pw_in(self.directory, atoms, self.parameters)
199 if 'bands' in properties :
200 write_bands_in(self.directory, atoms, self.parameters)

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/qeutil/writers.pyc in write_pw_in(d, a, p)
207 # Need to use symmetry properly
208 # create a dummy atoms object for primitive cell
--> 209 primcell=write_cell_params(fh,a,p)
210 if primcell :
211 # primitive cell has been found - let us use it

/home/lim520/miniconda/envs/nipy/lib/python2.7/site-packages/qeutil/writers.pyc in write_cell_params(fh, a, p)
164 fh.write(' C = %f,\n' % (C,))
165 elif sgn >= 75 :
--> 166 raise NotImplementedError
167 elif sgn ==1 :
168 # P1 symmetry - no special primitive cell signal to the caller

NotImplementedError:

With regards

Linu

Inconsistent treatment of strain throughout the code

First of all, thank you so much @jochym for sharing your Elastic code with the community.

Here I would like to report my concern about the convention of strain in the code.

To my best knowledge, in the common definition of the Voigt notation for strain, the strain values $\epsilon_1, \cdots, \epsilon_6$ are related to the tensor components as follows,

$$ \begin{pmatrix} \epsilon_{11} & \epsilon_{12} & \epsilon_{13} \\ \epsilon_{21} & \epsilon_{22} & \epsilon_{23} \\ \epsilon_{31} & \epsilon_{32} & \epsilon_{33} \end{pmatrix} = \begin{pmatrix} \epsilon_{1} & \epsilon_{6} / 2 & \epsilon_{5} / 2 \\ \epsilon_{6} / 2 & \epsilon_{2} & \epsilon_{4} / 2 \\ \epsilon_{5} / 2 & \epsilon_{4} / 2 & \epsilon_{3} \end{pmatrix}, $$

i.e., $\epsilon_4 = 2 \epsilon_{23}$, where we have the coefficient 2.
This is different from the convention for the stress (without the coefficient 2).
Formally, the shear components in the Voigt notation are equal to “engineering” shear strains.
At least in the context of elastic properties, the Voigt notation is not just replacements of indices.
And, the common stress-strain relation like (for cubic):

$$ \begin{align} \begin{pmatrix} \sigma_{1} \\ \sigma_{2} \\ \sigma_{3} \\ \sigma_{4} \\ \sigma_{5} \\ \sigma_{6} \\ \end{pmatrix} &= \begin{pmatrix} B_{11} & B_{12} & B_{12} & 0 & 0 & 0 \\ B_{12} & B_{11} & B_{12} & 0 & 0 & 0 \\ B_{12} & B_{12} & B_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & B_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & B_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & B_{44} \\ \end{pmatrix} \begin{pmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \epsilon_{3} \\ \epsilon_{4} \\ \epsilon_{5} \\ \epsilon_{6} \\ \end{pmatrix} \end{align} $$

is written using the above common Voigt notation.
I confirmed these points with several textbooks (Nye and Grimvall).

My concern is that the strain treatment is not well consistent through the code. Specifically,

• get_strain: It returns stress, not in the common Voigt notation, but actually just as an array of tensor components.
• regular, hexagonal, …: They must take the tensor-component strains rather than those in the Voigt notation.
• get_cart_deformed_cell (and thus get_elementary_deformations): here, the shear strains are engineering strains and thus the strain must be given in the Voigt notation.
• The $s$ in the document are not in the Voigt notation but in just tensor components. The least-square stress-strain relation in the document is also given, not for the common Voigt notation, but for the array of tensor components.

Particularly for the non-Voigt convention of get_strain may cause confusion for readers in realistic cases:

• If users see an 6-component array returned, they would naturally believe this is given in the Voigt notation, while presently this is not true.
• Further, if users would like to check the correctness of their obtained elastic constants via the "one-axis" approach (described in the last of this page), they may get twice larger values for shear elastic constants due to the usage of non-Voigt strain in get_strain.

To demonstrate the issue, I made a test using ASE EMT potential Cu_EMT.zip:

• In jochym0 I used get_elementary_deformations and get_elastic_tensor and obtained C11 = 172 GPa, C12 = 116 GPa, C44 = 90 GPa
• In jochym1 I used the one-axis approach using get_cart_deformed_cell and linear fitting and obtained C11 = 173 GPa, C12 = 115 GPa, C44 = 180 GPa, where C44 is twice larger obtained. This is because the strain is given not in the common Voigt notation.

To solve the issue, I wish to suggest:

• At least the get_strain method, the returned values should not be called “in the Voigt notation”.
• Ideally the get_strain method can be modified to return strain in the common Voigt notation. This also requires corresponding changes in other methods like regular, hexagonal, etc.

It would be very appreciated if you kindly think of these points above. Thank you so much for your help in advance.