The previous tutorials have calculated properties of silicon at fixed geometry, which in the case of a crystal includes the ionic positions as well as the lattice vectors. To optimize the ionic positions at fixed lattice vectors, command ionic-minimize suffices as in the molecular cases (see the Geometry optimization tutorial for example). This tutorial introduces lattice optimization using the lattice-minimize command.

We'll start with the essentially the same input as the previous tutorials, but in order to test the minimizer, we'll artificially introduce a strain:

#Save the following to Si.lattice
lattice face-centered Cubic 11.3   #this is deliberately about 10% too large

and specify the ionic positions to a separate file as well:

#Save the following to Si.ionpos
ion Si 0.00 0.00 0.00  1
ion Si 0.30 0.30 0.30  1   #deliberately perturbed (should have been 0.25)

and the input file for lattice optimization:

#Save the following to latticeOpt.in
include Si.lattice
include Si.ionpos

kpoint-folding 8 8 8
ion-species GBRV/$ID_pbe_v1.2.uspp
ion-species GBRV/$ID_pbe_v1.uspp
elec-cutoff 20 100

electronic-SCF                   #Electronic state optimization
ionic-minimize nIterations 10    #Internal geometry optimization
lattice-minimize nIterations 10  #Lattice geometry optimization
dump-name Si.$VAR
dump Lattice State               #Output state every lattice step

For silicon, symmetries completely determine the internal geometry (fractional ion coordinates), and therefore only lattice optimization is necessary. But in order to illustrate a more general situation, the above input starts from a perturbed internal geometry as well. Therefore, we specify commands for both ionic and lattice minimization.

Run JDFTx on this input file a few times:

mpirun -n 4 jdftx -i latticeOpt.in | tee latticeOpt.out
mpirun -n 4 jdftx -i latticeOpt.in | tee -a latticeOpt.out

Note that each run updates both Si.lattice and Si.ionpos, so that the next one starts at the final geometry of the previous one. However, we deliberately do not read in the wavefunctions of the previous run (using initial-state) in the subsequent runs. Once the lattice changes, the plane-wave basis within the specified kinetic energy cutoff will change making the wavefunctions incompatible.

Within each run, the basis is kept constant because otherwise the energy would change discontinuously with varying lattice vectors, making it impossible to use an efficient minimization algorithm. However, this means that the basis is inconsistent with the original elec-cutoff once the lattice vectors have changed from the input value. This is why it is essential to rerun lattice optimizations a few times, especially if the initial strain is large, and stop when the lattice changes negligibly (around 0.1%) in a given run. (This issue can be mitigated partially using Pulay corrections - see script calcPulay - but it is important to continue the lattice minimization over a few runs regardless.)

Examine the output file. Note the usual ionic minimize output, which performs electronic-SCF repeatedly at several geometries till the forces are within threshold. But now additionally, there is a lattice minimizer which runs an ionic minimization at each lattice geometry. At the end of the ionic minimization at one lattice geometry, it calculates and reports the stress and strain tensors, following which it updates the lattice vectors and iterates. Note that in the present case, the first run ends with a strain of -8.5%, while the second run ends with a strain less than 0.001%, so we can stop after the second run; but often more runs will be necessary.

You can easily see this inner-outer minimization structure using "grep Minimize latticeOpt.out". In this case it only prints IonicMinimize and LatticeMinimize since we used electronic-SCF; if we didn't specify electronic-SCF, you would also see ElecMinimize in that output. You can also visualize the convergence of the energy using the plotConvergence script. By default "plotConvergence latticeOpt.out" will plot the convergence every electronic iteration, but you can use "iterType=Ionic plotConvergence latticeOpt.out" or "iterType=Lattice plotConvergence latticeOpt.out" to examine the convergence with ionic or lattice opimization steps.

Check the final lattice vectors and internal geometry (Si.lattice and Si.ionpos). The final lattice is face-centered Cubic with a cubic lattice constant of 10.33 bohrs, which is only 0.7% larger than the experimental value of 10.263 bohrs (typical accuracy for the PBE exchange-correlation functional). Note that the ionic positions have a net translation to the ones we specified in previous tutorials: the optimizer moved both Si atoms in the unit cell toward each other to fix the perturbation we introduced. Of course, with periodic boundary conditions and the plane-wave basis, this makes no difference at all. If we wanted to keep the first Si atom pinned to the origin, we could have zero'd its moveScale factor (last parameter to ion command).

Finally, we can export the variation of the geometry using

createXSF latticeOpt.out latticeOpt.axsf Animated

and examine the animated XSF file using XCrysDen (since VESTA doesn't support animated XSFs).