Advanced relaxation parameters

Specifying the crystal structure

Inversion and symmetry

If there is a center of inversion, the code automatically tries to shift the atomic coordinates such that the crystal has inversion symmetry. If you want to inhibit that, put a line saying:

no restore_inversion

Finally, you can specify whether you want the symmetry-operations to be generated by the code, or if they should be read from the file SYMMETRYOPS. If number_symmetry_ops is followed by -1, then the symmetry operations are determined by the code itself. If number_symmetry_ops is $>0$, the matrices for the symmetry operation are read in from the file. The file format is the same as the printout in e.g. OUT.

Default: {\tt number\_symmetry\_ops -1}

If number_symmetry_ops is followed by -2, then the symmetry operations are determined by the code itself, but nonzero fractional translations are ignored. This could be useful for supercells.

Specifying an initial strain

You may optionally specify an initial deformation tensor using the format...

begin deformation
0.00 0.00 0.02
0.00 0.00 0.00
0.00 0.00 0.00
end deformation

This is the deformation tensor applied to your structure. It is defined by the following equation:

$$\mathbf{a} = (1 + \mathbf{D}) \cdot \mathbf{a_o}$$

where $\mathbf{a}$ are the deformed lattice vectors, and $\mathbf{a_o}$ are the lattice vectors specified in your input file. This can be an easy way to apply a shear strain to a cubic crystal, for instance.

There are times when it is useful to ensure that a given direction is held fixed during the relaxation process. This does not involve a constraint on the relaxation, as there are three rotational degrees of freedom. One can also arrange to maintain fixed a plane which contains the fixed direction.

To keep the $\mathbf{z}$ direction fixed, and also keep the zx plane from rotating, you could specify...

fixed_vector                0 0 1
other_vector_in_fixed_plane 1 0 0

Relaxing a structure

For the impact of the job variable on the relaxation process, see section 2.2.

Controlling the relaxation process

If job relax or job relax_from_checkpoint or job relax_recycle_H is chosen, the variable relax_what allows different ways to relax a structure:

relax_what force

or

relax_what stress

or

relax_what force_and_stress

Default: relax_what force_and_stress

Another variable allows to manipulate the way the relaxation is done:

Set

relax_how gradient

to simply follow along the forces and stress, without updating the inverse of the Hessian matrix.

Default: update the Hessian

With the variable relax_how the line-minimization can be altered also:

relax_how slow        # do 3 steps per line minimization always
relax_how normal      # do up to 3 steps, based on gradient prediction
relax_how fast        # do up to 2 steps, only if energy decreases

Default: relax_how normal

Note: This can be combined with the gradient option with a statement like relax_how normal_gradient

The parameter lambda_limit allows you to limit the size of relaxation steps to a certain number. Normally, lambda should be on the order of 1. Sometimes, especially in the beginning of the iteration, and if the H matrix is initialized to small, lambda can become to big, and the algorithm bombs out. It is then advisable to limit lambda with a line like that:

lambda_limit 4.0

Default: lambda_limit 1e6 (no limit)

The statement starting_distortion followed by nine floating point numbers allows to set an initial distortion of the lattice. This feature is only for debugging, and not recommended for use.

With version 5.1.6, there are two relaxation algorithms with the addition of a limited memory BFGS algorithm. The benefit is more in stability and in general better efficiency (only based on a small test suite of systems) than for memory. The input keyword is

relax_method   default=1   original BFGS quasi-netwon algorithm
                  2        limited memory BFGS

Pressure

An external pressure (units are GPa) can be applied with:

relax_pressure 10.0

If an entry like

relax_pressure adjust

is found, the pressure is chosen from the first iteration, taking it to be 1/3 of the trace of the stress.

Default: relax_pressure 0.0

Constraining the relaxation

Paratec offers the possibility to fix atoms to their input coordinates, meaning they will not move during the relaxation process. This is simply done by zeroing out the Hellman-Feynman forces. To fix atoms to their original positions, use the stay_put structure. The following example fixes the positions of atoms number 1, 3, and 10.

begin stay_put
stay_put 1
stay_put 3
stay_put 10
end stay_put

Default: all atoms are relaxed.

As of paratec 4.5, paratec also offers the option of fixing certain strain coordinates while allowing the others to change during the relaxation. This is done by using the keyword stay_put_lattice within the stay_put structure. The following example fixes the $\epsilon_{xx}$ and $\epsilon_{yy}$ strains, while allowing the z axis to expand or contract.

begin stay_put
stay_put_lattice xx
stay_put_lattice yy
end stay_put

Default: all strain coordinates are relaxed.

NOTE: The syntax of the stay_put structure has changed in version 5.0. The end at the end is now an end stay_put.

Initialization of the inverse of the Hessian matrix

The Broyden scheme used for the relaxation of the forces and the stress updates an inverse of the Hessian matrix, which essentially describes the enthalpy parabola around the minimum. There are some open parameters which determine the values to which the matrix is initialized. They basically determine the length of the trial steps during the first few iterations:

relax_eps_factor   0.1     initialization perfactor for the stress part
of H
relax_coord_factor 0.1     initialization perfactor for the force part
of H

Try with the values above for the first run. Adjust later such that the first few steps give a lambda of about 0.5 to 1.

Default:

relax_coord_factor 1.0
relax_eps_factor   assuming the bulk modulus is about 1 Mbar

Alternatively, you can specify the physical quantities to which the above two variables correspond. This is the recommended way:

estimated_bulk_modulus  1.0                 Bulk modulus in MBar
estimated_optical_phonon_frequency 12.0     TO phonon frequency f at
                                            Gamma in THz)

Diamond would have a bulk modulus of about 4.4 MBar and around 35 THz for the optical phonon frequency.

For molecular solids, the option

decouple_coordinates yes

causes a different initialization of the Hessian matrix, which decouples the coordinates of the molecule from the lattice vectors to first order. This means that if the lattice vectors change 10 percent, the cartesian coordinates will change by about 1 percent. IF THIS OPTION IS USED, THE SYMMETRY OF THE CRYSTAL IS NOT PRESERVED!.

Default:

decouple_coordinates no