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.

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 2019.3 release of the Amsterdam Modeling Suite:

  • 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.
  • UFF
    An implementation of the Universal Force Field, a simple non-reactive force field covering the entire periodic table.
  • 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.
  • ADF
    A wrapper around the standalone ADF program, allowing it to be used as an engine from within the AMS driver.
  • 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.

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 "$ADFBIN"/startpython /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 $ADFBIN/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 ADFSpectra 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 [True | False]
      Functional string
      a1 float
      a2 float
      s6 float
      s8 float
      sr6 float
   End
   D4Dispersion
      Enabled [True | False]
      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:False
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:False
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 [True | False]
   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:False
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 [1] from Parrinello and Rahman, the extra energy term due to a non-isotropic external stress [2] 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).

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

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 reference unit cell, \(h\) are the lattice vectors in matrix form, \(h_0\) are the lattice vectors of the reference unit cell in matrix form and \(\mathbb{1}\) is the identity matrix.