How to do uniaxial/biaxial compressive/tensile strain engineering calculations using DFT (VASP, Quantum ESPRESSO, etc.)

Strain engineering has proven quite exciting in exploring quantum phenomena in condensed matter physics, materials science, and other fields. For instance, one could investigate the effect of applying a tensile (i.e., stretching) or compressive (i.e., squeezing) force on the electronic properties (band structure, magnetic exchange, conductivity, etc.) of a 2D material.

To do strain calculations using ab initio calculations using software like VASP and Quantum ESPRESSO, all one needs to do is:

1. Get initial structure for the unstrained material. You can use The Materials Project or other source for this. It is important the lattice vectors are reasonable. So, either use values from the literature (preferably experimental, e.g., from Single-Crystal X-ray Diffraction), or do a full structure optimization by also relaxing lattice vectors (e.g., using ISIF = 3 in VASP). If using literature values, I recommend you relax atomic positions too before proceeding, as that could help subsequent relaxations happen faster!

2. Use the above lattice vectors to get lattice vectors of the strained structure. For unstrained lattice vectors a=(a_1,a_2,a_3),b=...,c=... written as a 3\times 3 matrix V:

V=\begin{pmatrix} a_1 & a_2 & a_3 \\ b_1 & b_2 & b_3 \\ c_1 & c_2 & c_3 \end{pmatrix},

you can get the strained lattice vectors using:

\begin{pmatrix} 1+x & 0 & 0 \\ 0 & 1+y & 0 \\ 0 & 0 & 1+z \end{pmatrix}\cdot V

Above, x,y,z are the percentages by which to strain the x, y, and z components of the lattice vectors by. For example, for 3% compressive biaxial strain in the xy directions, you would have x=y=-3/100, z=0. For 1% tensile uniaxial strain in the y direction, you would have x=0, y=+1/100, z=0.
At the end of this post, I have simple Mathematica code that exemplifies these calculations.

3. Relax atomic positions for each strained structure. For this, change the lattice vectors of the relaxed unstrained structure in 1 above, leaving atomic positions unchanged. In VASP, you would directly modify the lattice vectors in the POSCAR file for this, as highlighted in the screenshot below. Then relax the atomic positions only, keeping lattice vectors fixed (e.g., using ISIF = 2 in VASP).

4. Use the above relaxed strained structure to calculate the properties you desire! For instance, you could do scf calculations for a range of strain percentages, and plot the total energy vs strain.

  • If doing this, as a sanity check, you should expect the unstrained structure (i.e., 0%) to have the lowest energy, as it should be naturally preferred. However, this might not always be the case due to the software’s numerical idiosyncrasies. For instance, the different choices of basis expansion sets might cause 1% biaxial strain to have the lowest energy, instead of the expected 0%-minimum case depicted below:

Above, although I used VASP to exemplify, similar steps hold for other software. With Quantum ESPRESSO, for instance, use calculation = relax to optimize the structure by relaxing atomic positions only. Use calculation = vc-relax to additionally relax lattice vectors.


Mathematica code for strained lattice vector calculation:

strainMyVectors[x_, y_, z_] = {{1 + x, 0, 0}, {0, 1 + y, 0}, {0, 0, 1 + z}};
unstrainedVectors = {{7.600, 0.000, 0.000}, {-3.800, 6.582, 0.000}, {0.000, 0.000, 15.000}}; (* = V in discussion above *)

(* 3% biaxial compressive strain in the xy directions *)
strainMyVectors[-3/100, -3/100, 0].unstrainedVectors

(* 1% tensile strain in the y direction *)
strainMyVectors[0, +1/100, 0].unstrainedVectors

Leave a Comment