17. Structure relaxation calculation

In this section, we show the rudiments of relaxation calculations.

17.1. Example: Si crystal bond length

To perform a structure relaxation, set the keyword info.calculationType to relaxation. Let’s consider a silicon primitive cell in which one atom is out of place, a calculation which boils down to computing the Si crystal bond length. Create a text file si_real_rlx containing the following keyword definitions.

info.savepath = 'results/si_real_rlx'
info.calculationType = 'relaxation'
atom.element = [1 1]
atom.fracxyz = [0 0 0;0.26 0.27 0.28]
domain.latvec = 10.25*[0.0 0.5 0.5;0.5 0.0 0.5;0.5 0.5 0.0]
domain.lowres = 0.4
kpoint.gridn = [3 3 3]
element(1).species = 'Si'
element(1).path = './Si_TM_LDA.mat'
functional.list = {'XC_LDA_X','XC_LDA_C_PW'}
sr.method = 'qn'
sr.tol = 2e-4
sr.relaxedstructure = 'si2_relaxed.xyz'
eigensolver.emptyBand = 1

Note that keywords setting the relaxation engine and other relaxation related parameters begin with sr.

We described the rationale behind the above input file.

  • info.calculationType is set to ’relaxation’ so that RESCU optimizes the structure;

  • atom.fracxyz contains the atoms’ fractional coordinates. Here, the second atom is displaced by [0.01 0.02 0.03] from it’s equilibrium position at [0.25 0.25 0.25];

  • sr.method is a string determining the algorithm used to predict the equilibrium positions. The allowed values are ’sd’, ’cg’ and ’qn’, which stand for steepest descent, conjugate gradient and quasi-Newton respectively. The default value is ’qn’. One may also use ’fire’ to use the FIRE algorithm, a variation of damped molecular dynamics;

  • sr.tol sets the tolerance for the structure relaxation engine in atomic units. More precisely, RESCU stops the structure optimization when all force coordinates fall below the tolerance. The default value is 1e-4, which corresponds to roughly 0.005 eV/Å;

  • sr.relaxedstructure is a string containing the path to the file which contains the relaxed atomic coordinates. The default value is ./relaxedstructure.xyz;

  • eigensolver.emptyBand controls the number of unoccupied bands included in a calculation. Although the ground state density includes only occupied bands, the inclusion of unoccupied bands is advisable (e.g. for numerical stability) or crucial (e.g. for metallic materials). The default value depends on system size, but it is rather ‘safe’, i.e. rather large. Since relaxation calculations are computationally demanding, it is best to use the minimal amount of bands possible. For a simple silicon crystal, a single extra band will do.

Various additional parameters are described by the function inputDescription. One parameter worth mentioning is sr.moveableAtomList, which controls which atom coordinates are included in the relaxation’s engine degrees of freedom. When sr.moveableAtomList is not set, RESCU relaxes all atoms except the first one, taking advantage of translational invariance to eliminate three degrees of freedom.

The calculation is executed as follows.

rescu -i si_real_rlx

During the relaxation procedure, RESCU performs several self-consistent calculations whose output goes according to info.outfile as usual. The relaxation procedure itself may be monitored via the file relaxationprocess.txt, which looks as follows. If the positions or forces vary too wildly, it is thus possible to stop the execution and prevent waste of computational resources.

Relaxation Iteration: 000
Total energy = -7.905e+00, normi(force) = +2.559e-02
Atomi Spec:       xCrd       yCrd       zCrd       xFrc       yFrc       zFrc
    1   Si: +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00
    2   Si: +2.818e+00 +2.767e+00 +2.716e+00 -2.559e-02 -1.781e-02 -8.477e-03


Line search iteration: 001
Total energy = -7.912e+00, normi(force) = +1.712e-04
sk.force = -1.834e-04 (norm(sk) = +5.559e-03, alpha = +1.000e+00)
Atomi Spec:       xCrd       yCrd       zCrd       xFrc       yFrc       zFrc
    1   Si: +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00
    2   Si: +2.563e+00 +2.562e+00 +2.563e+00 -1.581e-04 -1.747e-05 -1.712e-04
Structure relaxation calculation finished!
Total time: 9e+01 seconds

At each iteration of the relaxation engine, RESCU prints the iteration number, the total energy and norm of the forces in atomic units. The forces associated with frozen degrees of freedom (see section 1.2) are set to zero. The norm of the forces is the infinity norm, which is essentially the force coordinate with the largest magnitude. RESCU also prints the atom index, the atom species, the coordinates and the forces (in Cartesian coordinates in atomic units) for every atom. The relaxation engines proceed by so-called line searches (except for the FIRE algorithm). For each iteration of a line search, RESCU prints similar information, but it will additionally print the force projected on the line search direction (norm(sk)) and the current line search step (alpha).

When the relaxation procedure is completed, RESCU prints the resulting atomic coordinates to an XYZ file. For instance,


AtomType                      X                      Y                      Z
      Si +0.000000000000000e+00 +0.000000000000000e+00 +0.000000000000000e+00
      Si +1.356583441021671e+00 +1.356077852198500e+00 +1.356628976586721e+00

The atomic positions are Cartesian coordinates in units of Angstrom.

17.2. Constrained relaxation

It is often advantageous to choose which degrees of freedom should be relaxed as opposed to which atoms. In RESCU, the three Cartesian coordinates of every atom are the independent degrees of freedom. Which is to be relax is again specified via the keyword sr.moveableAtomList, only that this time it is a \(n_{\text{atom}}\times 3\) Boolean array. It is usually inconvenient to write such an array in the main input file. Like the atomic coordinates, which are specified in a file given by atom.xyz, the active degrees of freedom may be specified in an XYZ file given by sr.moveableAtomList. RESCU reads columns 8 to 10 in that file and fixes the value of the degrees of freedom marked 0. Note that the first column contains the element label, the columns 2 to 4 the atomic coordinates and the columns 5 to 7 the initial magnetic moment, and hence the need for 10 columns. Depending on the spin treatment, the initial magnetic moment may be ignored, but it must be included nevertheless. For example, the following xyz-file will minimize \(f_x^1\), \(f_y^2\) and \(f_z^3\) (where the superscript denotes the atom index)

#element x y z sx sy sz fx fy fz
Si 0.26           0         0 0 0 0 1 0 0
Si 0         4.9052    5.1052 0 0 0 0 1 0
Si 5.1052         0    5.3052 0 0 0 0 0 1
Si 5.1052    5.1052         0 0 0 0 0 0 0
Si 2.5526    2.5526    2.5526 0 0 0 0 0 0
Si 2.5526    7.6578    7.6578 0 0 0 0 0 0
Si 7.6578    2.5526    7.6578 0 0 0 0 0 0
Si 7.6578    7.6578    2.5526 0 0 0 0 0 0