# Engines¶

The engines are the core of the Amsterdam Modeling Suite: While the AMS driver
steers the calculation over the potential energy surface in e.g. a
geometry optimization or molecular
dynamics calculation, the engines calculate energies and
gradients and in this way *define* the PES on which the driver works.

```
Engine header
```

`Engine`

Type: Block Description: The input for the computational engine. The header of the block determines the type of the engine.

The engine used for an AMS calculation is selected and configured with the special
`Engine`

block in the AMS input:

```
Engine DFTB
... input for the DFTB engine ...
EndEngine
```

Here the type of engine, e.g. `DFTB`

as in the example above, is specified on
the line that opens the block. Note that the `Engine`

block ends with
`EndEngine`

, and is in this way different from all the other blocks in the AMS
input, which close just with `End`

. The content of the engine block is what we
call the “engine input”. Generally the engine input consists of a series of
blocks and keywords, and looks just like the AMS driver input. However, many
engines have a lot of options and keywords, which are documented in a separate
engine manual. In other words: This AMS driver manual documents all the keywords
outside of the `Engine`

block, while the individual engine manuals document
the contents of the engine block.

## Available engines¶

The following engines are available in the 2020 release of the Amsterdam Modeling Suite:

- ADF
- A DFT engine particularly strong in understanding and predicting structure, reactivity, and spectra of molecules.

- BAND
- An atomic-orbital based DFT engine aimed at periodic systems (crystals, slabs, chains) but supporting also molecular systems.

- DFTB
- An engine implementing Density Functional based Tight-Binding, a fast approximation to DFT.

- ReaxFF
- An engine for modeling chemical reactions with atomistic potentials based on the reactive force field approach.

- MLPotential
- Machine learning potentials (machine learning force fields).

- ForceField
- An engine implementing classical non-reactive force fields such as UFF (a non-reactive force field covering the entire periodic table).

- Hybrid
- Hybrid engine, for embedding and QM/MM calculations that combine multiple engines.

- MOPAC
- An engine wrapping the MOPAC code, a general-purpose semiempirical molecular orbital package for the study of solid state and molecular structures and reactions.

- External
- A flexible scripting interface that allows advanced users to use external modeling programs as engines in AMS.

- LennardJones
- A simple toy engine implementing a Lennard-Jones potential.

## Summary of engine capabilities¶

Some options/properties can only be computed in combination with some of the engines (i.e., not all engines support all features). These tables summarize the capabilities of the engines of the Amsterdam Modeling Suite:

Feature | ADF | BAND | DFTB | ReaxFF | MLPot. | ForceField | MOPAC |
---|---|---|---|---|---|---|---|

All elements available | ✓ | ✓ | ✓ [1] | ✓ [2] | ✓ | ||

Non-periodic systems | ✓ | ✓ | ✓ | ✓ | ✓ [3] | ✓ | ✓ |

1D periodic systems | ✓ | ✓ | ✓ | ✓ | |||

2D periodic systems | ✓ | ✓ | ✓ | ✓ | |||

3D periodic systems | ✓ | ✓ | ✓ | ✓ [3] | ✓ | ✓ | |

Charged molecular systems | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | |

External electric field | ✓ | ✓ | ✓ | ✓ | |||

External point charges | ✓ | ✓ | ✓ |

Feature | ADF | BAND | DFTB | ReaxFF | MLPot. | ForceField | MOPAC |
---|---|---|---|---|---|---|---|

Atomic charges | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | |

Bond orders | ✓ | ✓ | ✓ | ✓ | ✓ | ||

Dipole gradients | ✓ | ✓ | ✓ | ✓ | |||

Dipole moment | ✓ | ✓ | ✓ | ✓ | |||

Elastic tensor | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | |

Gradients / Forces | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |

Hessian | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |

Infrared spectra | ✓ | ✓ | ✓ | ✓ | |||

Molecule detection | ✓ | ✓ | ✓ | ✓ | ✓ | ||

Normal modes | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |

PES point character | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |

Phonons | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | |

Polarizability | ✓ | ||||||

Raman | ✓ | ||||||

Stress tensor | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | |

Thermodynamic properties | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |

VCD | ✓ | ✓ | ✓ | ✓ | |||

VROA | ✓ |

Note

The features/options in the following tables are **engine-specific**, and are described in the corresponding engine’s manual. The input for these options should be specified in the Engine block section of the input file.

Feature | ADF | BAND | DFTB | ReaxFF | MLPot. | ForceField | MOPAC |
---|---|---|---|---|---|---|---|

EFG | ✓ (doc) | ✓ (doc) | |||||

Electronic transport | ✓ (doc) | ✓ (doc) | ✓ (doc) | ||||

Energy decomposition analysis | ✓ (doc) | ✓ (doc) | |||||

ESR/EPR | ✓ (doc) | ||||||

Hyperpolarizabilities | ✓ (doc) | ||||||

NMR | ✓ (doc) | ||||||

Orbital energies | ✓ | ✓ | ✓ | ✓ | |||

pKa | ✓ (doc) | ||||||

Solvation models | ✓ (doc) | ✓ (doc) | ✓ (doc) | ✓ (doc) | |||

Spin-polarization | ✓ (doc) | ✓ (doc) | ✓ (doc) | ✓ (doc) | |||

UV-VIS optical spectra | ✓ (doc) | ✓ (doc) |

Feature | ADF | BAND | DFTB | ReaxFF | MLPot. | ForceField | MOPAC |
---|---|---|---|---|---|---|---|

Density of states (DOS) | ✓ (doc) | ✓ (doc) | |||||

Dielectric function | ✓ (doc) | ||||||

Effective mass | ✓ (doc) | ✓ (doc) | |||||

Electronic band structures | ✓ (doc) | ✓ (doc) | |||||

k-space sampling | ✓ (doc) | ✓ (doc) | |||||

X-ray form factors | ✓ (doc) |

[1] | All elements are available for the GFN1-xTB model, but other DFTB parameter sets are not parametrized for all elements. |

[2] | All elements are available with UFF (universal force field), but other force fields are not parametrized for all elements. |

[3] | (1, 2) MLPotential backends may support non-periodic and/or 3D-periodic calculations. |

## External programs as engines¶

The AMS driver allows running external programs as an engine. In this way users can combine the functionality in the AMS driver (tasks and PES point properties) with the energies and gradients of any molecular modeling program they have access to. For example, the graphical user interface supports setting up VASP as an external engine to the AMS driver.

In general, the interfacing between the AMS driver and the external program has to be done by the user in form of a small script, which allows users to hook up any external program without access to the source code of AMS. The graphical user interface of the Amsterdam Modeling Suite can then be used to analyze the results of these calculations.

An external engine is specified with `Engine External`

. The keyword
`Execute`

is used to specify the command that is run to execute the external program:

```
Engine External
Execute /path/to/my_interface_script.sh
EndEngine
```

The command can in principle be anything, as it will just be executed as is by the system shell. However, it should not use relative paths (e.g. to files in the directory where the input file is). We recommend writing the interfacing script in Python and using the Python interpreter that ships with AMS:

```
Engine External
Execute "$AMSBIN"/amspython /path/to/my_python_interface_script.py
EndEngine
```

AMS then starts running and for every geometry prepares a folder in which the
external engine is supposed to run. This is the folder in which the interface
script specified with the `Execute`

key is executed (so any relative paths
are relative to that folder). AMS puts two files into this folder:

```
system.xyz
request.json
```

The `system.xyz`

just contains the geometry AMS wants the external engine to
calculate. It is an extended format XYZ file with the
`VEC1`

, `VEC2`

, `VEC3`

extension at the end for periodic systems, e.g.
diamond would look like this:

```
2
C -0.51292147 -0.51292147 -0.51292147
C 0.51292147 0.51292147 0.51292147
VEC1 0.00000000 2.05168587 2.05168587
VEC2 2.05168587 -0.00000000 2.05168587
VEC3 2.05168587 2.05168587 0.00000000
```

The `request.json`

file is just a small JSON file that specifies what exactly
AMS wants the external engine to calculate:

```
{
"title": "GOStep28",
"quiet": false,
"gradients": true,
"stressTensor": false,
"hessian": true,
"dipoleMoment": false,
"properties": true,
"prevResults": "GOStep27"
}
```

The job of the interfacing script is now to read these files, run the external
program and convert its output into a format understood by AMS. Generally these
are simple text files with the name of the property and the extension `.txt`

.
The bare minimum the interfacing script needs to produce is the file
`energy.txt`

containing a single number, i.e. the total energy in atomic
units (Hartree). Other properties are optional, and it is easiest to go through
the `request.json`

entries one by one to see what AMS might request and what
the interfacing script could produce in response.

`title`

- Just a title for this particular engine run. It can be passed on to the external program if desired, or can just be ignored.
`quiet`

- Whether AMS wants the external engine to write to standard output. This can be ignored in principle, but that might lead to really incomprehensible text output files of AMS if the external engine has to be called many times, e.g. for numerical derivatives.
`gradients`

- Whether or not to calculate nuclear gradients. The interface script should
put the gradients in a file called
`gradients.txt`

with nAtoms lines of 3 real numbers each, in atomic units, i.e. Hartree/Bohr. Note that AMS wants the gradients, not forces (beware the - sign!). `stressTensor`

- Whether to calculate the stressTensor for periodic systems. Should be written
to
`stresstensor.txt`

in atomic units. `hessian`

- Whether to calculate the Hessian, that is just the second derivative of the
energy with respect to the nuclear coordinates,
*without*applying any mass weighing to it. If the Hessian has been calculated, it should be put in`hessian.txt`

as a 3 nAtoms x 3 nAtoms matrix in atomic units. `dipoleMoment`

- If true, calculate the dipole moment and put it in
`dipolemoment.txt`

in atomic units, in one line with three numbers. `properties`

- This is set to true if AMS considers this “geometry” important and wants the
engine to calculate further properties that the user might be interested in.
In practice this is set to “true” for e.g. the final converged step in a
geometry optimization, so that the user can then let the engine calculate
e.g. the band structure, which one would not want to do at all the steps
*during*the optimization. AMS can’t do anything with the properties that the engine might calculate, but the files will remain on disk for people to inspect them. `prevResults`

- This is the title of a previous similar calculation that the engine has
already performed. These results can be accessed in
`../$prevResults/`

, so for the example above`GOStep28`

can access the results from the previous step in the geometry optimization in`../GOStep27/`

. This is just the directory in which the interfacing script was run when the`title`

field was set to`GOStep27`

, so files that were written back then are still accessible. They can in principle be used to restart for example the SCF of the engine from step to step. Of course all of that has to be done by the interfacing script. The AMS driver does not know anything about how to restart the external program and can only point the interfacing script to the right location.

That is really all there is to the external engine: AMS prepares a folder with
`system.xyz`

and `request.json`

and runs the user’s interfacing script in
there, which has to take care of preparing the input for the external engine,
running it, and putting the results in the text files that AMS expects, e.g.
`gradients.txt`

.

Note for properties that are in one way or another derivatives of the energy, it
is generally ok if the external engine does not calculate what was requested by
the AMS driver in `request.json`

. If AMS requests, for example, the gradients
from the external engine, but then does not find the `gradients.txt`

in the
directory after the interfacing script has run, it will just assume that the
engine was not capable of calculating the gradient analytically. AMS will then
just do the gradient numerically by rerunning the external engine for displaced
geometries, reading only the energy from `energy.txt`

. In this sense it is
only absolutely required for the external engine to produce the energy, the rest
can be done numerically by AMS if required. It is of course best to let the
engine do as much as possible, especially if it implements analytical
derivatives. Note that currently AMS can not calculate the Hessian numerically
for engines that do not provide gradients. This is just a technical limitation,
as it is of course possible to do a second derivative numerically, but it is
just not implemented in AMS yet. (And it would also be a very slow way to
calculate a Hessian.)

In addition to the `Execute`

keyword that specifies the interfacing script,
it is also possible to use the `ExecuteAtEnd`

keyword to specify another
command to be run at the end of the calculation, for example after the last
step of a geometry optimization or molecular dynamics simulation. This command
is run in the results directory.

Moreover, the `Engine External`

block can also contain some information about the
capabilities of the external engines:

```
Engine External
Execute {...}
ExecuteAtEnd {...}
Supports
DipoleMoment {true|false}
PeriodicityNone {true|false}
PeriodicityChain {true|false}
PeriodicitySlab {true|false}
PeriodicityBulk {true|false}
End
EndEngine
```

The normal engines that come with AMS (e.g. DFTB and BAND) produce the engine
output files with extension `.rkf`

in the results directory, see
here. These files are also produced when an external
engine is used and the information on them (anything related to the shape of the
PES at that point, e.g. normal modes, phonons, …) can be visualized normally
with the graphical interface. In addition to each engine output `.rkf`

file,
external engines will also produce a correspondingly named folder per engine
file, which is just the working directory of the interfacing script for that
particular invocation of the external program. These folders just contain the
full output of the external program and anything that the interfacing script
might have produced. In this way users still have access to all results from the
external program, even if these results were not communicated back to the AMS
driver.

This last point is probably best illustrated with a simple example. Consider the following job that uses an external engine to do a linear transit calculation of ethane, going from the staggered to the eclipsed configuration, calculating normal modes at all converged points along the path:

```
AMS_JOBNAME=ethane_torsion $AMSBIN/ams << EOF
Task PESScan
System
Atoms
C 0.00000000 0.00000000 0.76576000
C 0.00000000 0.00000000 -0.76576000
H -0.88668938 0.51193036 1.16677000
H 0.88668938 0.51193036 1.16677000
H 0.00000000 -1.02386071 1.16677000
H 0.00000000 1.02386071 -1.16677000
H -0.88668938 -0.51193036 -1.16677000
H 0.88668938 -0.51193036 -1.16677000
End
End
PESScan
CalcPropertiesAtPESPoints True
ScanCoordinate
nPoints 5
Dihedral 3 1 2 6 60.0 0.0
End
End
Properties
NormalModes True
End
Engine External
...
EndEngine
EOF
```

If we run this job and look into the results folder, we will find the standard
`ams.log`

and `ams.rkf`

as well as the usual engine result files
`PESPoint(1).rkf`

to `PESPoint(5).rkf`

. Just as if we had used one of the
native AMS engines, like DFTB. Each of these files can be opened in AMSspectra
to visualize the normal modes for this particular point. For an external engine
we additionally have one folder per engine file, so for this example we would
have `PESPoint(1)/`

to `PESPoint(5)/`

. These are the folders in which the
interfacing script ran for these particular points, so they contain all the
native output files of the external program.

## Toy engines¶

The AMS driver comes with a simple built-in toy engine that implements a Lennard-Jones potential. This can sometimes be useful for testing, as many properties of the Lennard-Jones gas/liquid/solid can be calculated analytically and compared to the results from AMS. Note that the potential is exactly the same for all elements, i.e. the N-N bond has exactly the same strength as the He-He bond.

The Lennard-Jones engine only has three keywords, which define the shape of the potential:

```
Engine LennardJones
RMin float
Eps float
Cutoff float
EndEngine
```

`Cutoff`

Type: Float Default value: 15.0 Unit: Angstrom Description: The distance at which the interaction is truncated.

`Eps`

Type: Float Default value: 1.0 Unit: Hartree Description: The depth of the potential well.

`RMin`

Type: Float Default value: 1.0 Unit: Angstrom Description: The distance of the potential minimum.

## Engine add-ons¶

Engine add-ons can be used to augment the results returned from an engine.

### 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.

Warning

If you want to include pressure in molecular dynamics calculations, you should not use this engine addon, but use the MD-specific pressure option.

### 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).

[4] | M. Parrinello, and A. Rahman, Polymorphic transitions in single crystals: A new molecular dynamics method, Journal of Applied Physics 52, 7182 (1981) |

[5] | The energy term due to a non-isotropic external stress is:
\[E_\text{stress} = p_h (V - V_0) + \frac{1}{2} V_0 Tr[(\sigma- p_h \mathbb{1}) (h_0^{T-1} h^T h h_0^{-1} - \mathbb{1})]\]
where: \(\sigma\) is the external stress tensor, \(p_h\) is the hydrostatic pressure associated with \(\sigma\) (i.e. the average of the diagonal elements of \(\sigma\)), \(V\) is the volume of the unit cell (for 2D periodic systems this is the area of the cell, and for 1D periodic systems is the length of the cell), \(V_0\) is the volume of the |

### 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
```