Engines¶

Dispersion corrections¶

For engines that do not natively support Grimme’s D3 and D4 dispersion corrections, the D3Dispersion and D4Dispersion engine add-ons can be used to add such corrections.

EngineAddons
   D3Dispersion
      Damping [BJ | Zero]
      Enabled Yes/No
      Functional string
      a1 float
      a2 float
      s6 float
      s8 float
      sr6 float
   End
   D4Dispersion
      Enabled Yes/No
      Functional [...]
      Verbosity [Silent | Normal | Verbose | VeryVerbose]
      a1 float
      a2 float
      s6 float
      s8 float
      s9 float
   End
End

EngineAddons

D3Dispersion

Type:	Block
Description:	This block configures the add-on that adds the Grimme D3 dispersion correction to the engine’s energy, gradients, and stress tensor.

Damping

Type:	Multiple Choice
Default value:	BJ
Options:	[BJ, Zero]
Description:	Type of damping: BJ (Becke-Johnson) or Zero. BJ is recommended for most applications.

Enabled

Type:	Bool
Default value:	No
Description:	Enables the D3 dispersion correction addon.

Functional

Type:	String
Default value:	PBE
Description:	Use the D3 parameterization by Grimme for a given xc-functional. Accepts the same values as the –func command line option of the official dftd3 program. Note: the naming convention is different from elsewhere in the AMS suite. For example, BLYP should be called b-lyp.

a1

Type:	Float
Description:	The a1 parameter. Only used if Damping is set to BJ. If set, it overwrites the a1 value for the chosen functional.

a2

Type:	Float
Description:	The a2 parameter. Only used if Damping is set to BJ. If set, it overwrites the a2 value for the chosen functional.

s6

Type:	Float
Description:	The s6 parameter, global scaling parameter. If set, it overwrites the s6 value for the chosen functional.

s8

Type:	Float
Description:	The s8 parameter. If set, it overwrites the s8 value for the chosen functional.

sr6

Type:	Float
Description:	The sr6 parameter. Only used if Damping is set to Zero. If set, it overwrites the sr6 value for the chosen functional.

D4Dispersion

Type:	Block
Description:	This block configures the addon that adds the Grimme D4(EEQ) dispersion correction to the engine’s energy and gradients.

Enabled

Type:	Bool
Default value:	No
Description:	Enables the D4 dispersion correction addon.

Functional

Type:	Multiple Choice
Default value:	PBE
Options:	[HF, BLYP, BPBE, BP86, BPW, LB94, MPWLYP, MPWPW91, OLYP, OPBE, PBE, RPBE, REVPBE, PW86PBE, RPW86PBE, PW91, PW91P86, XLYP, B97, TPSS, REVTPSS, SCAN, B1LYP, B3LYP, BHLYP, B1P86, B3P86, B1PW91, B3PW91, O3LYP, REVPBE0, REVPBE38, PBE0, PWP1, PW1PW, MPW1PW91, MPW1LYP, PW6B95, TPSSH, TPSS0, X3LYP, M06L, M06, OMEGAB97, OMEGAB97X, CAM-B3LYP, LC-BLYP, LH07TSVWN, LH07SSVWN, LH12CTSSIRPW92, LH12CTSSIFPW92, LH14TCALPBE, B2PLYP, B2GPPLYP, MPW2PLYP, PWPB95, DSDBLYP, DSDPBE, DSDPBEB95, DSDPBEP86, DSDSVWN, DODBLYP, DODPBE, DODPBEB95, DODPBEP86, DODSVWN, PBE02, PBE0DH, B1B95, MPWB1K, REVTPSSH, GLYP, REVPBE0DH, REVTPSS0]
Description:	Use the D4 parameterization by Grimme for a given xc-functional.

Verbosity

Type:	Multiple Choice
Default value:	Silent
Options:	[Silent, Normal, Verbose, VeryVerbose]
Description:	Controls the verbosity of the dftd4 code. Equivalent to the –silent, –verbose, and –very-verbose command line switches of the official dftd4 program.

a1

Type:	Float
Description:	The a1 parameter, see D4 article. The physically reasonable range for a1 is [0.0,1.0]. If set, it overwrites the a1 value for the chosen functional.

a2

Type:	Float
Description:	The a2 parameter, see D4 article. The physically reasonable range for a2 is [0.0,7.0]. If set, it overwrites the a2 value for the chosen functional.

s6

Type:	Float
Description:	The s6 parameter, see D4 article. The physically reasonable range for s6 is [0.0,1.0]. If set, it overwrites the s6 value for the chosen functional.

s8

Type:	Float
Description:	The s8 parameter, see D4 article. The physically reasonable range for s8 is [0.0,3.0]. If set, it overwrites the s8 value for the chosen functional.

s9

Type:	Float
Description:	The s9 parameter, see D4 article. If set, it overwrites the s9 value for the chosen functional.

Pressure¶

Pressure can be included using the following keyword (this can only be used for 3D periodic systems):

EngineAddons
   Pressure float
End

EngineAddons

Pressure

Type:	Float
Default value:	0.0
Unit:	GPa
Description:	Add a hydrostatic pressure term to the engine’s energy and stress tensor. Can only be used for 3D periodic boundary conditions.

The engine’s energy will include the following extra term: \(P \cdot V\), where \(V\) is the volume of the unit cell. The engine’s stress tensor will include the following extra term: \(P \cdot \mathbb{1}\), where \(\mathbb{1}\) is the identity matrix. This can only be used for 3D periodic boundary conditions (i.e. bulk). The energy and stress tensor printed in output and written to the .rkf binary file will include these extra terms (i.e. the printed stress tensor is the sum of the internal stress and the external stress do to pressure).

When studying the effect that pressure has on the structure and properties of your system, one should generally start by optimizing the structure including the lattice vectors under pressure. Properties such as phonons or elastic tensor can then be computed at the relaxed geometry. If you are investigating phase transitions under pressure (or if you simply expect some symmetry breaking) you should disable symmetry and/or perturb the initial geometry of your system.

Non-isotropic external stress¶

An non-isotropic external stress tensor can be included by using the following keywords:

EngineAddons
   ExternalStress
      StressTensorVoigt float_list
      UpdateReferenceCell Yes/No
   End
End

EngineAddons

ExternalStress

Type:	Block
Description:	This block configures the addon that adds external stress term to the engine’s energy and stress tensor.

StressTensorVoigt

Type:	Float List
Unit:	a.u.
GUI name:	External stress tensor
Description:	The elements of the external stress tensor in Voigt notation. One should specify 6 numbers for 3D periodic system (order: xx,yy,zz,yz,xz,xy), 3 numbers for 2D periodic systems (order: xx,yy,xy) or 1 number for 1D periodic systems.

UpdateReferenceCell

Type:	Bool
Default value:	No
Description:	Whether ot not the reference cell should be updated every time the system changes (see documentation).

The energy and stress tensor printed in output and written to the .rkf binary file already include the corresponding extra terms, i.e. the printed stress tensor is the sum of the internal stress and the input-specified external stress.

When studying the effect that an external stress has on the structure and properties of your system, one should generally start by optimizing the structure including the lattice vectors under the external stress (depending on the magnitute of the applied external stress, you might have to adjust the stress energy per atom convergence threshold). If your system is symmetric, you should disable symmetry when optimizing structures under external stress. Be aware that the geometry optimization might go completely astray (e.g. the material will break apart) if you apply a) too large shear stress or b) too large tension stress (too large negative stresses for the diagonal values).

Following this paper [4] from Parrinello and Rahman, the extra energy term due to a non-isotropic external stress [5] is defined with respect to a reference unit cell, which in our case is the unit cell at the beginning of the simulation. If, during the simulation, large deformations of the unit cell occur, the above mentioned energy expression is only approximately consistent with the stress tensor. This affects the calculation differently depending on whether the stress tensor is computed by the engine or by AMS via numerical differentiation: a) If the stress tensor is computed directly by the engine, for large unit-cell deformations the energy might increase during the optimization; b) If the stress tensor is computed by AMS via numerical differentiation, the actual final value of the stress tensor might not match perfectly the external stress specified in the input (for large cell deformations, this error can be in the order of 10%).

After the optimization under external stress is converged, it is therefore good practice to validate the results. To do this, you should compute the stress tensor at the optimized geometry by performing a single point calculation without applying the external stress. The values in the computed stress tensor should have opposite sign compare to the external stress applied during the optimization. If the values differ too much, then you can run a second geometry optimization under external stress starting from the optimized geometry.

An alternative option is to set the key ExternalStress%UpdateReferenceCell to True; this will update the reference unit cell at every optimization step, effectively changing the definition of the energy expression at every geometry optimization step. The energy might not consistently go down during the optimization, but the resulting internal stress will match much better the applied external stress. This option should only be used during geometry optimizations (i.e. it should not be used when computing properties such as the elastic tensor).

Atom energies¶

This add-on adds an element-dependent energy for each atom in a system, effectively applying an energy shift. It does not affect any other property.

EngineAddons
   AtomEnergies # Non-standard block. See details.
      ...
   End
End

EngineAddons

Type:	Block
Description:	This block configures all the engine add-ons.

AtomEnergies

Type:	Non-standard block
Description:	Add an element-dependent energy per atom. On each line, give the chemical element followed by the energy (in atomic units).

The below example will add -0.5 Hartree for every H atom and -5.0 Hartree for every C atom in the system to the total energy. The add-on will ignore any atoms of other elements:

EngineAddons
   AtomEnergies
      H -0.5
      C -5.0
   End
End

Restraints¶

A restraint is potential energy function (a spring) attached to a certain coordinate, for example, an interatomic distance, the minimum of the potential energy being at the specified optimal value. A restraint can have one or two parameters: a ForceConstant and, for some types, a F(Inf) value. The ForceConstant parameter corresponds to second derivative of the restraint potential energy \(\frac {d^{2}V(x)}{dx^2}\) for any \(x\) (harmonic restraints) or only at \(x = 0\) (other restraints). Here, \(x\) is a deviation from the restraint’s optimal value.

Geometry restraints can be added using the following input block:

Restraints
   Profile [Harmonic | Hyperbolic | Erf]
   fInfinity float
   Distance atomIdx1 atomIdx2 OptValue {Profile} {fInfinity}
   Angle
   Dihedral
   SumDist
   DifDist
End

The shape of the restraint potential function is defined by its profile type (Harmonic, Hyperbolic, or Erf) and the appropriate parameter(s). The Harmonic profile is most suitable for geometry optimizations but may result is very large forces that can be problematic in molecular dynamics. For MD simulations the Hyperbolic or Erf may be more suitable because the restraint force is bounded by a user-defined value. A Harmonic restraint is defined as \(V(x)=C \frac {x^2}{2}\). A Hyperbolic restraint is defined as \(V(x) = a (\sqrt{1+(b x)^2}-1)\). An Erf restraint is defined as \(V(x) = a (b x \cdot erf(b x) + \frac {e^{-(b x)^2}-1}{\sqrt{\pi}})\), which corresponds to an indefinite integral of the error function \(\frac {dV(x)}{dx} = a b \cdot erf(b x)\). The \(a\) and \(b\) parameters are calculated from the FInfinity and ForceConstant values. The Erf profile is very similar to Hyperbolic however its derivative converges to F(Inf) much faster.

The following keys define global parameters for all restraints.

Restraints

Profile

Type:	Multiple Choice
Default value:	Harmonic
Options:	[Harmonic, Hyperbolic, Erf]
GUI name:	Default restraint profile
Description:	Select the default type of restraint profile. The harmonic profile is most suitable for geometry optimizations but may result is very large forces that can be problematic in molecular dynamic. For MD simulations the Hyperbolic or Erf may be more suitable because the restraint force is bounded by a user-defined value. A per-restraint profile type can be specified after the ForceConstant value on the corresponding restraint line.

FInfinity

Type:	Float
Default value:	1.0
GUI name:	Default F(inf)
Description:	Specify the default asymptotic value for the restraint force for the Hyperbolic and Erf profiles, in Hartree/Bohr or Hartree/radian. A per-restraint value can be specified after the profile type on the corresponding restraint line.

The following keys define individual restraints. The force constant, profile type and the F(Inf) value can be set for each restraint individually.

Restraints

Distance

Type:	String
Recurring:	True
Description:	Specify two atom indices followed by the distance in Angstrom and, optionally, by the ForceConstant (default is 1.0 in a.u.), profile type and F(Inf) (in a.u.). This restraint will try to keep the distance between the two specified atoms at the given value. For periodic systems this restraint follows the minimum image convention.

Angle

Type:	String
Recurring:	True
Description:	Specify three atom indices i j k followed by an angle in degrees and, optionally, by the ForceConstant (default is 0.3 in a.u.), profile type and F(Inf) (in a.u.). This restraint will try to keep the i-j-k angle at the given value. For periodic systems this restraint follows the minimum image convention.

Dihedral

Type:	String
Recurring:	True
Description:	Specify four atom indices i j k l followed by an angle in degrees and, optionally, by the ForceConstant (default is 0.1 in a.u.), profile type and F(Inf) (in a.u.). This restraint will try to keep the i-j-k-l dihedral angle at the given value. For periodic systems this restraint follows the minimum image convention.

DifDist

Type:	String
Recurring:	True
Description:	Specify four atom indices i j k l followed by the distance in Angstrom and, optionally, by the ForceConstant (default is 1.0 in a.u.), profile type and F(Inf) (in a.u.). This restraint will try to keep the difference R(ij)-R(kl) at the given value. For periodic systems this restraint follows the minimum image convention.

SumDist

Type:	String
Recurring:	True
Description:	Specify four atom indices i j k l followed by the distance in Angstrom and, optionally, by the ForceConstant (default is 1.0 in a.u.), profile type and F(Inf) (in a.u.). This restraint will try to keep the sum R(ij)+R(kl) at the given value. For periodic systems this restraint follows the minimum image convention.

A default value for the force constant depends on the restraint type: 1.0 Hartree/Bohr for Distance, SumDist and DifDist, 0.3 Hartree/radian for Angle and 0.1 Hartree/radian for Dihedral.

Below is an example if the Restraints block.

Restraints
   Profile Hyperbolic     # Change the default profile type
   fInfinity 10.0         # Change the asymptotic value for the restraint force away from optimum
#  Type     Atoms    OptValue  FC   Profile  F(Inf)
   Distance 1 2        5.0     1.0  Erf      1.0
   Angle    1 2 3     90.0
   Dihedral 4 1 2 3  180.0     0.1
#  The next two together are equivalent to imposing two distance restraints: R(23) to 0.85 Angstrom and R(14) to 0.65 Angstrom
   SumDist  1 4 2 3    1.5     # Keep R(14)+R(23) as close to 1.5 Angstrom as possible
   DifDist  1 4 2 3    0.2     # Keep R(14)-R(23) as close to 0.2 Angstrom as possible
End