JDFTx  1.1.2
Geometry optimization

Lattice relaxation

We will first look at the lattice relaxation of silicon.

Here is the input file. It is essentially the same as the one in the Density of States tutorial page.

#---------------- Lattice Constant of Silicon ------------------
lattice \
    0   0.5 0.5 \
    0.5 0   0.5 \
    0.5 0.5 0

latt-scale 10.77 10.77 10.77

ion-species Si.fhi

coords-type lattice    
ion Si  0.250  0.250  0.250 0
ion Si  0.000  0.000  0.000 0


kpoint-folding 4 4 4

elec-ex-corr gga

dump-name silicon.$VAR

lattice-minimize \
    nIterations 10 \
    energyDiffThreshold 1e-6

electronic-minimize \
    energyDiffThreshold 1e-7

The experimental (cubic) lattice constant of Silicon is 5.43 Angstroms, or roughly 10.26 bohrs. We run the calculation, starting from a 5 % strained initial state. (jdftx -i silicon.in -o silicon.out).

Just before the electronic minimization starts, you should see something like this:

--------- Lattice Minimization ---------
Minimization of dimension 1 over strains spanned by:
[  0.57735  0  0  ]
[  0  0.57735  0  ]
[  0  0  0.57735  ]

This means that jdftx has detected the symmetries of your crystal and realized that there is only one lattice parameter. The basis matrices for the strain tensor are expressed in lattice coordinates, so a diagonal strain tensor corresponds to scaling all three fcc lattice vectors by the same constant. The normalization of the basis matrices is such that sum of all squares always add to 1 (i.e. 3 * 0.57734^2 = 1.0)

Calculations with lattice minimization tend to produce rather long output files, so it can be useful to use search tools (such as UNIX's grep) to look at them.

Lets grep the output file for lattice minimization iterations. grep LatticeMinimize silicon.out produces the following output:

LatticeMinimize: Iter:   0  Etot: -7.927814070374968e+00  |grad|_K:  1.940e-01  alpha:  1.000e+00
LatticeMinimize: Iter:   1  Etot: -7.937043278195619e+00  |grad|_K:  3.265e-02  alpha:  3.966e-01, linmin = -1.000e+00, cgtest =  1.000e+00
LatticeMinimize: Encountered beta<0, resetting CG.
LatticeMinimize: Iter:   2  Etot: -7.937240053671994e+00  |grad|_K:  2.902e-04  alpha:  3.685e-01, linmin =  1.000e+00, cgtest = -1.000e+00
LatticeMinimize:    Bad step direction: g.d > 0.
LatticeMinimize:    Undoing step.
LatticeMinimize:    Step failed: resetting search direction.
LatticeMinimize: Iter:   3  Etot: -7.937240076434499e+00  |grad|_K:  2.798e-04  alpha:  0.000e+00
LatticeMinimize: Iter:   4  Etot: -7.937240091707153e+00  |grad|_K:  1.913e-04  alpha:  6.715e-01, linmin =  1.000e+00, cgtest = -1.000e+00
LatticeMinimize: Converged (|Delta Etot|<1.000000e-06 for 2 iters).

We see from that the lattice minimizer converged to 1 microHartrees in just 4 iterations.

Lets now grep for the strain tensor, i.e grep Strain -A 3 silicon.out. -A 3 just tells UNIX to report the 3 lines after finding the word Strain (case sensitive).

Strain Tensor = 
[          0          0          0  ]
[          0          0          0  ]
[          0          0          0  ]
--
Strain Tensor = 
[ -0.0444153          0          0  ]
[          0 -0.0444153          0  ]
[          0          0 -0.0444153  ]
--
Strain Tensor = 
[ -0.0513619          0          0  ]
[          0 -0.0513619          0  ]
[          0          0 -0.0513619  ]
--
Strain Tensor = 
[ -0.0513619          0          0  ]
[          0 -0.0513619          0  ]
[          0          0 -0.0513619  ]
--
Strain Tensor = 
[ -0.0512534          0          0  ]
[          0 -0.0512534          0  ]
[          0          0 -0.0512534  ]

It seems that the crystal strained down by roughly 5.1 % from the initial size, i.e. almost back to the exprimental value. Lets look at the lattice vectors, grep "R = " -A 3 silicon.out

R = 
[          0      5.385      5.385  ]
[      5.385          0      5.385  ]
[      5.385      5.385          0  ]
--
R = 
[          0      5.385      5.385  ]
[      5.385          0      5.385  ]
[      5.385      5.385          0  ]
--
R = 
[          0    5.14582    5.14582  ]
[    5.14582          0    5.14582  ]
[    5.14582    5.14582          0  ]
--
R = 
[          0    5.10842    5.10842  ]
[    5.10842          0    5.10842  ]
[    5.10842    5.10842          0  ]
--
R = 
[          0    5.10842    5.10842  ]
[    5.10842          0    5.10842  ]
[    5.10842    5.10842          0  ]
--
R = 
[          0      5.109      5.109  ]
[      5.109          0      5.109  ]
[      5.109      5.109          0  ]

The calculated lattice constant is 10.218 Bohrs (5.109x2), or 5.41 Angstroms, roughly 0.02 Angstorms off the experimental value.

Tips

  • Whether you are using lattice or cartesian coordinates becomes very important when using the lattice minimizer. Cartesian coordinates for ions are allowed in lattice minimizer, but not recommended.
  • As the unit cell is strained, JDFTx keeps the number of basis functions the same. However, because the reciprocal lattice vectors are also changing, this leads to changes in energy. If pulay files for the pseudopotentials are given, JDFTx can partially correct for this on the fly. Pulay files contain information about how the total energy of an atom changes as the volume of the G-sphere (sphere containing all basis functions) changes. Pulay files can be generated using the calcPulay script provided with the source code. The default setting in JDFTx is that it looks for a *.pulay* file in the same directory as the pseudopotential file, but this behaviour can be overwritten.
  • For very large strains of the unit cell, the pulay correction may become insufficient. In such cases, it is best to restart the calculation with the converged lattice vectors. There is no easy rule of thumb for this, but anything around 10% strain of the initial lattice vectors should be under suspect and restarted. If the pulay corrections were sufficient, the lattice vectors will not change and the calculation will exit after few iterations.