Elastic properties#
See this reference cite:PhysRevB.65.104104.
We seek the elastic constant tensor that relates stress (σ) and strain (ε) via \(\sigma = c \epsilon \). The stress and strain are six component vectors, so \(c\) will be a 6 × 6 symmetric matrix.
Fe elastic properties#
As with molecular vibrations, we need a groundstate geometry. Let us get one for BCC Fe.
from vasp import Vasp
from ase.lattice.cubic import BodyCenteredCubic
atoms = BodyCenteredCubic(symbol='Fe')
for atom in atoms:
atom.magmom = 3.0
# from vasp.vasprc import VASPRC # Deprecated in new vasp
# # VASPRC['mode'] = None # Deprecated in new vasp # Deprecated in new vasp
import logging
log = logging.getLogger('Vasp')
#log.setLevel(logging.DEBUG)
calc = Vasp(label='bulk/Fe-bulk',
xc='PBE',
kpts=(6, 6, 6),
encut=350,
ispin=2,
isif=3,
nsw=30,
ibrion=1,
atoms=atoms)
print(atoms.get_potential_energy())
print(atoms.get_stress())
-16.45852388
[ 0.04574869 0.04574869 0.04574869 0. -0. -0. ]
Ok, now with a relaxed geometry at hand, we proceed with the elastic constants. This is accomplished with IBRION=6 and ISIF > 3 in VASP. See this reference (from the VASP page) Y. Le Page and P. Saxe, Phys. Rev. B 65, 104104 (2002) https://doi.org/10.1103/PhysRevB.65.104104 for more details.
from vasp import Vasp
calc = Vasp(label='bulk/Fe-bulk').clone('bulk/Fe-elastic')
calc.set(ibrion=6, #
isif=3, # gets elastic constants
potim=0.05, # displacements
nsw=1,
nfree=2)
print(calc.potential_energy)
-22.39158131
Now, the results are written out to the OUTCAR file. Actually, three sets of moduli are written out 1) the elastic tensor for rigid ions, 2) the contribution from allowing the atoms to relax, and 3) the total elastic modulus, all in kBar.
SYMMETRIZED ELASTIC MODULI (kBar)
Direction XX YY ZZ XY YZ ZX
--------------------------------------------------------------------------------
XX 2803.5081 1622.6085 1622.6085 0.0000 0.0000 0.0000
YY 1622.6085 2803.5081 1622.6085 0.0000 0.0000 0.0000
ZZ 1622.6085 1622.6085 2803.5081 0.0000 0.0000 0.0000
XY 0.0000 0.0000 0.0000 866.8792 0.0000 0.0000
YZ 0.0000 0.0000 0.0000 0.0000 866.8792 0.0000
ZX 0.0000 0.0000 0.0000 0.0000 0.0000 866.8792
--------------------------------------------------------------------------------
and
ELASTIC MODULI CONTR FROM IONIC RELAXATION (kBar)
Direction XX YY ZZ XY YZ ZX
--------------------------------------------------------------------------------
XX 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
YY 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
ZZ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
XY 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
YZ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
ZX 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
--------------------------------------------------------------------------------
TOTAL ELASTIC MODULI (kBar)
Direction XX YY ZZ XY YZ ZX
--------------------------------------------------------------------------------
XX 2803.5081 1622.6085 1622.6085 0.0000 0.0000 0.0000
YY 1622.6085 2803.5081 1622.6085 0.0000 0.0000 0.0000
ZZ 1622.6085 1622.6085 2803.5081 0.0000 0.0000 0.0000
XY 0.0000 0.0000 0.0000 866.8792 0.0000 0.0000
YZ 0.0000 0.0000 0.0000 0.0000 866.8792 0.0000
ZX 0.0000 0.0000 0.0000 0.0000 0.0000 866.8792
--------------------------------------------------------------------------------
Let us write a small code here to extract the Total elastic moduli from the OUTCAR file. First we get the line where the total elastic moduli start, then take the six lines that start three lines after that. Finally we parse out the matrix elements and cast them as floats.
import numpy as np
EM = []
with open('bulk/Fe-elastic/OUTCAR') as f:
lines = f.readlines()
for i, line in enumerate(lines):
if line.startswith(' TOTAL ELASTIC MODULI (kBar)'):
j = i + 3
data = lines[j:j+6]
break
for line in data:
EM += [[float(x) for x in line.split()[1:]]]
print(np.array(EM))
[[2119.0655 1022.4304 1022.4304 -0. -0. -0. ]
[1022.4304 2119.0655 1022.4304 0. 0. 0. ]
[1022.4304 1022.4304 2119.0655 0. -0. 0. ]
[ -0. 0. 0. 2564.0346 0. 0. ]
[ -0. 0. -0. 0. 2564.0346 0. ]
[ -0. 0. 0. 0. 0. 2564.0346]]
Fe is in a BCC crystal structure, which is high in symmetry. Consequently, many of the elements in the matrix are equal to zero.
See http://www.nist.gov/data/PDFfiles/jpcrd34.pdffor a lot of detail about Fe-Ni alloys and general theory about elastic constants. In the next section, we show how the code above is integrated into mod:Vasp.
Al elastic properties#
First, the relaxed geometry.
from vasp import Vasp
from ase.lattice.cubic import FaceCenteredCubic
atoms = FaceCenteredCubic(symbol='Al')
calc = Vasp(label='bulk/Al-bulk',
xc='PBE',
kpts=(12, 12, 12),
encut=350,
prec='High',
isif=3,
nsw=30,
ibrion=1,
atoms=atoms)
print(calc.potential_energy)
-14.97542738
Ok, now with a relaxed geometry at hand, we proceed with the elastic constants. This is accomplished with IBRION=6 and ISIF ≥ 3 in VASP.
from vasp import Vasp
calc = Vasp(label='bulk/Al-bulk').clone('bulk/Al-elastic')
calc.set(ibrion=6, #
isif=3, # gets elastic constants
potim=0.015, # displacements
nsw=1,
nfree=2)
calc.calculate()
EM = calc.get_elastic_moduli()
print(EM)
c11 = EM[0, 0]
c12 = EM[0, 1]
B = (c11 + 2 * c12) / 3.0
print(B)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[18], line 13
5 calc.set(ibrion=6, #
6 isif=3, # gets elastic constants
7 potim=0.015, # displacements
8 nsw=1,
9 nfree=2)
11 calc.calculate()
---> 13 EM = calc.get_elastic_moduli()
15 print(EM)
17 c11 = EM[0, 0]
File ~/vasp/vasp/mixins/analysis.py:297, in AnalysisMixin.get_elastic_moduli(self)
294 parts = line.split()
295 if len(parts) >= 7:
296 # Skip first column (XX, YY, etc.)
--> 297 row = [float(x) for x in parts[1:7]]
298 matrix.append(row)
300 # Convert from kBar to GPa
File ~/vasp/vasp/mixins/analysis.py:297, in <listcomp>(.0)
294 parts = line.split()
295 if len(parts) >= 7:
296 # Skip first column (XX, YY, etc.)
--> 297 row = [float(x) for x in parts[1:7]]
298 matrix.append(row)
300 # Convert from kBar to GPa
ValueError: could not convert string to float: 'of'
This example shows the basic mechanics of getting the elastic constants. The C44 constant above is too low, and probably we need to check these constants for convergence with respect to kpoints, planewave cutoff, and maybe the value of POTIM.
Manual calculation of elastic constants#
It is possible to manually calculate single elastic constants; you just need to know what strain corresponds to the elastic constant.
For the C11 elastic constant in a cubic system, we simply strain the cell along the x-axis, and then evaluate the second derivative at the minimum to calculate C11 like this.
\(C_{11} = \frac{1}{V_C}\frac{\partial^2 E^{tot}}{\partial \gamma^2}\)
from vasp import Vasp
from ase.lattice.cubic import FaceCenteredCubic
import numpy as np
import matplotlib.pyplot as plt
DELTAS = np.linspace(-0.05, 0.05, 5)
calcs = []
volumes = []
for delta in DELTAS:
atoms = FaceCenteredCubic(symbol='Al')
cell = atoms.cell
T = np.array([[1 + delta, 0, 0],
[0,1, 0],
[0, 0, 1]])
newcell = np.dot(cell, T)
atoms.set_cell(newcell, scale_atoms=True)
volumes += [atoms.get_volume()]
calcs += [Vasp(label='bulk/Al-c11-{}'.format(delta),
xc='pbe',
kpts=(12, 12, 12),
encut=350,
atoms=atoms)]
Vasp.run()
energies = [calc.potential_energy for calc in calcs]
# fit a parabola
eos = np.polyfit(DELTAS, energies, 2)
# first derivative
d_eos = np.polyder(eos)
print(np.roots(d_eos))
xfit = np.linspace(min(DELTAS), max(DELTAS))
yfit = np.polyval(eos, xfit)
plt.plot(DELTAS, energies, 'bo', xfit, yfit, 'b-')
plt.xlabel('$\delta$')
plt.ylabel('Energy (eV)');
[ 0.00727102]