vasp 5.2.2 Elastic Constants

[dpsi]

I am looking at the output for the elastic constants but I'm getting some inconsistencies.
I have set ISIF=3 and IBRION=6 in my INCAR file in order to calculate elastic constants of a pre-relaxed structure (relaxed cell parameters and ionic coords).

I noticed that by default, just two strains for each unit-cell coordinate are applied on either side of the equilibrium. I tried to change that by setting NFREE=4 and now there are four strains on either side of equilibrium but POTIM should be halved in order to have the max/min strains the same as before.

The problem is that the results do not agree with NFREE=4 and the default (no NFREE specified).

Here's an idea of the disagreement (rigid ions, i.e. SYMMETRIZED ELASTIC MODULI) for a slightly tetragonal cell:

NFREE =4:
unrelaxed ions: C_11=C_33=2.18, C_22=2.24, C_12=1.02, C_13=0.68,
C_44=C_55=1.42, C_66=1.09

no NFREE specified
C_11=C_33=2.70, C_22=2.25, C_12=1.36,C_13=1.14,
C_44=C_55=1.25, C_66=0.95

The results for relaxed ions are also in severe disagreement.

The case of no NFREE specified agrees with my independent calculations for this system, obtained by a script which applies strains. Why does NFREE=4 give such strange results?

Many thanks

[dpsi]

A quick follow-up to the posting above.

Here's a more concrete comparison with a well-defined system.

The system is bcc-Fe using the Fe_pv.PAW.GGA pseudopotential
for which I got a lattice constant of 2.834 Ang. I used 12x12x12 k-points and 400 eV as a quick estimate, although the problem persists with more accurate calculations.

Experiments give C_11=2.34, C_12=1.36, C_44=1.18 Mbar
Phys Rev 122, 1714 (1961)
while another ab initio rsult from PRB 66, 024113 (2002) gives
C_11=2.85, C_12=1.39, C_44=1.00 Mbar.

Independently, by fitting to E, I get
C_11=2.84, C_12=1.52, C_44=1.04
and fitting to stress,
C_11=2.74, C_12=1.42, C_44=1.05.

Now, with vasp 5.2 and NFREE unspecified, I get reasonable agreement with the above:
C_11=2.91, C_12=1.39, C_44=0.92. POTIM was set equal to 0.01.

However, with NFREE=4 I get :
C_11=3.54, C_12=2.01, C_44=1.32. POTIM was changed to 0.005
to have the same strains as above.

Any clue as to why I'm getting such different values?



<span class='smallblacktext'>[ Edited Mon Jul 20 2009, 02:38PM ]</span>

[dpsi]

An update:
I fitted the stress vs strain data from the OUTCAR file on my own and the results do not agree with the elastic constants in the NFREE=4 case. They agree with the other values obtained either with NFREE=2 or separately.

[tholme]

I am also having a problem calculating the elastic constants of bcc Fe. I've tried a reduced unit cell (1 atom), 1 unit cell (2 atoms), and 2x2x2 unit cells (16 atoms), all with fine k-point meshes (17x17x17 for the first two cases and 9x9x9 for the third supercell). When I do a volume relaxation (ISIF=7) with PREC=HIGH, and tight convergence criteria (EDIFF=1e-6, EDIFFG=1e-7), I get a nice convergence to a stress=0 lattice constant. However, when I do a series of manual calculations, changing my lattice constant and calculating energy at each step to find the minimum energy vs. lattice constant, I find a different lattice constant than the ISIF=7 calculation, for all 3 choices of supercell.

Even more disturbingly, either lattice constant I pick, for all 3 supercell choices, show the following problem. I do a series of manual changes to the x-component of the supercell dimension and calculate the energy versus delta x-perturbation. The energy is not at a minimum for my delta=0, even though I had minimized energy versus volume. This is true for all 3 supercells, and both lattice constants (from ISIF=7 and the manual relaxation). Instead of seeing the energy at a minimum for delta=0, I'll see a changing energy that is minimized for some delta as high as 2.5%.

Any ideas for why this could be?

NB: I'm using the PBE GGA POTCARs, both the Fe_pv and Fe POTCAR files. I'm using ISPIN=2 for a magnetic iron, and I'm using GGA+U with LDAUJ=1 and LDAUU=6.