3.6. FRC Forcefield Files in MedeA
Contents
| download: | pdf |
|---|
3.6.1. Introduction
As described in the previous section Forcefields for Materials Simulations,
forcefields, or interatomic potentials describe interatomic interatomics in classical molecular simualtions.
The .frc format in MedeA is used for forcefield parameter files,
which contain all the necessary information for using a specific forcefield for running a classical
molecular dynamics (MD) or Monte Carlo (MC) simulation.
This section briefly describes the contents of MedeA frc files.
3.6.2. Getting Started
Lines starting with an exclamation point ! are considered as comments and no action is taken for those ones.
The .frc formatted forcefield file consists of several sections,
which start with a # character.
The order with which these sections appear in the .frc file is irrelevant.
In the beginning of each section there might be some lines starting with a > character,
which are special comments related to that specific section.
3.6.3. Forcefield Type
The first line of an .frc formatted file is the identification of the forcefield:
!MD forcefield 1
Currently, the type must be 1
3.6.4. Versions and References
Forcefield files may be augmented with the addition of new parameters or the refinement of old ones. In order to keep track of the changes and be able to easily and quickly locate them, version numbers and references are used, placed in the first two columns of a forcefield parameter entry. Version numbers have the general form of a release number, followed by a revision number, e.g. 1.2 refers to version 1 and revision 2. For the same entry, the data of the entry with the highest version number are included and the rest all ignored, while they remain in the file to facilitate keeping track of parameter changes.
For example, an entry for the nonbond parameters might look like this:
#nonbond(9-6) pcff+ 200
@type r-eps
@combination sixth-power
!Ver Ref I r eps
!--- --- ----- ------------- -----------
1.0 1 c0 3.8540 0.0150
2.0 19 c0 3.7500 0.0070
which means that new nonbond parameters have been added in version 2.0 (ref. 19) for the atom type “c0” since version 1.0 (ref. 1). The new parameters (version 2.0) will be the ones used in the simulation but the old parameters remain in the forcefield file, for tracking changes. Why these parameters have changed and by whom is something that may be mentioned in the corresponding reference. The references have the sole purpose of providing comments and documentation to the user and are disregarded by when the forcefield file is being read-in. References have the following form:
#reference 30
@Author D. Rigby
@Date 10-Sep-2014
Original published compass molecular oxygen parameters re-optimized to fit density and DHvap along entire saturation curve.
In order for MedeA to use the latest (highest) version for all parameters, in the beginning of the frc file (after the forcefield type definition), the latest (highest) version number needs to be defined:
#version pcff+.frc 1.0 12-August-2010
#version pcff+.frc 1.1 7-September-2010
#version pcff+.frc 2.0 15-September-2010
#version pcff+.frc 3.1 26-August-2013
The dates mentioned in this section help to keep track of when the changes were applied and are not used somewhere. Each time a new set of parameters is introduced in the frc file, the user should add this version number in the beginning of the frc file, as shown above. In this example, the 3.1 version is the highest version that will be used in case of multiple entries for the same set of parameters.
3.6.5. Description
A short description of the forcefield and some information on its used may be provided in the frc file. Description sections have the following form:
#description
Finnis-Sinclair potential for Zr
These potentials were imported from the NIST Interatomic Potentials
Repository Project on 13 June 2015 Please see
http://www.ctcms.nist.gov/potentials/ for more details of the project,
and the web pages below for more information about these particular
potentials. NIST requests that if you find these potentials useful, and
publish results with them, that you please reference
C.A. Becker, et al., "Considerations for choosing and using force fields
and interatomic potentials in materials science and engineering,"
Current Opinion in Solid State and Materials Science, 17, 277-283
(2013). http://www.ctcms.nist.gov/potentials
as well as the appropriate reference for the potential itself:
M.I. Mendelev and G.J. Ackland, "Development of an interatomic potential
for the simulation of phase transformations
in zirconium," Phil. Mag. Lett., 87, 349-359 (2007). DOI:
10.1080/09500830701191393.
The file at NIST was 'Zr\_1.eam.fs'. Information from the file's header:
Source: Potential #1 from [M.I. Mendelev and G.J. Ackland, Phil. Mag.
Letters 87, 349-359 (2007).]
Contact information: mendelev@ameslab.gov
Thursday, Nov 29, 2007 The potential was taken from v2\_3\_hcp (in
C:\\SIMULATION.MD\\Zr\\Results\\v2\_3)
Please see the following web page for more detail:
http://www.ctcms.nist.gov/potentials/Zr.html
3.6.6. Inclusion of Other FRC File(s)
There may be different “flavors” of a given forcefield, i.e. different
versions of it that have been published as such or that result from
additions and extensions to a forcefield. In such cases, it is most
convenient to include one frc file into another. This allows the use of
parts of the included forcefield(s) without the need to duplicate the
data therein and facilitates handling of the forcefield.
To include an frc file, this needs to be declared in the beginning of the
frc file, after the versions and description sections:
#include pcff.frc
3.6.7. Definition of Forcefield
In one frc file, there will be one or more defined forcefields.
One table listing the interactions used in the forcefield is necessary for
each of the forcefields defined.
The section of the definition of each forcefield starts with a define header:
#define cvff\_nocross\_nomorse
followed by a table with the list of all the sections that need to be included in this forcefield, including versions, reference numbers, function, label(s):
!Ver Ref Function Label
!---- --- --------------------------------- ------
2.0 18 atom\_types cvff
1.0 1 equivalence cvff
2.0 18 auto\_equivalence cvff\_auto
1.0 1 hbond\_definition cvff
2.0 18 quadratic\_bond cvff cvff\_auto
2.0 18 quadratic\_angle cvff cvff\_auto
2.0 18 torsion\_1 cvff cvff\_auto
2.0 18 out\_of\_plane cvff cvff\_auto
1.0 1 nonbond(12-6) cvff
3.0 32 bond\_increments cvff
3.0 31 templates cvff
In the above example, the first column provides the version for each interaction type, the second column provides the corresponding reference, the third column provides the name of the function (see the following sections of this document), and the fourth and last column provides the label of the section(s) included. More than one label (each one must be corresponding to one and only one section in the frc file for the specific function) can be included in the label column. In the above example, there is only one label “cvff” for “atom_types” but there are two “cvff” and “cvff_auto” for “quadratic_angle”. The name that appears in the header of the forcefield definition table is the one that will be used by . In the above example, that will be: cvff_nocross_nomorse.
In case there are more than one forcefields defined in the frc file, when reading the forcefield in , all the forcefields will be read and the default one will be the one used, unless otherwise specified:
#define clayff-dioctahedral default
Forcefield files may define a type for a forcefield:
#force_field_type mesoscale
This type is used by MedeA to control what system a forcefield can be applied to.
3.6.7.1. Atom Types
The atom types that are utilized in the forcefield, are defined in the atom_types section. For each atom type, the version number, the reference, the type name, the mass and the number of connections are the fields that need to be filled in, while there is an extra field where a comment may optionally be introduced. It is reading-in this information and is using it to assign parameters to the atom types, as those are defined in the frc file
#atom\_types pcff+ 200
!Ver Ref Type Mass Element connection Comment
!--- --- ----- ----- ------- ---------- --------------------
For all-atom forcefields, each atom is assigned an atom type, for which all relevant forcefield parameters need to be assigned.
For united-atom forcefields, a united atom type is assigned to one of the atoms present in the system (e.g. for a CH3 group, a CH3 atom type is assigned to the carbon atom of this group) while the hydrogens, that are “lumped” together into the united atom particle, are assigned a “UnitedH” atom type and are disregarded from the simulation.
3.6.7.2. Equivalences
Similar atom types may have equivalent behavior, as far as certain interactions are concerned. To use the minimum number of parameters that need to be defined in the parameter file and avoid large replications that increase the size of the forcefield file and the difficulty of editing it and safely adding new parameters, an equivalence table is introduced in the frc file.
In the equivalences table, there is one column for each different type of interaction (e.g. nonbond, bond, angle, torsion, charges etc) of a given forcefield. For each atom type for which some of the interactions may be equivalenced with those of another atom type, the name of that atom type is defined. That is done for every different interaction term, as shown in the following example:
#equivalence pcff+ 200
! Equivalences
! ------------------------------------------
!Ver Ref Type NonB Bond Angle Torsion OOP
!--- --- ----- ----- ----- ----- ------- -----
1.0 1 c0 c0 c c c c
In this example, for the atom_type “c0” unique nonbond parameters will be used, i.e. the ones introduced specifically for this atom type in the corresponding section of the forcefield file. For the bond, angle, torsion and out-of-plane terms, the parameters that will be used for “c0” are those that are defined for the atom type “c”; therefore there is no need to define any of the interaction parameters for the bond, angle, torsion and out-of-plane terms for “c0”.
3.6.7.3. Scaling
Some forcefields apply scaling factors to interactions between atoms which also interact in some other way. Usually, non-bond interactions between atoms which are bonded directly or to a common third atom are not included in the forcefield energy.
The scaling table permits to explicitly and globally specify the scaling factors for 1,2-, 1,3- and 1,4-interactions for both electrostatic and Van der Waals energies:
#scaling Martini
!Ver Ref I J K L Coulomb VdW
!---- --- ---- ---- ---- ---- ------- -------
1.0 1 * * 0.0 0.0
1.0 1 * * * 1.0 1.0
1.0 1 * * * * 1.0 1.0
In this example the 1,2-interactions are excluded from the forcefield energy, while the 1,3- and 1,4-interactions are fully included.
3.6.7.4. Nonbond (van der Waals) Interactions
Nonbond interactions are describing the interactions between nonbonded atoms (inter- or intra-molecular interactions).
3.6.7.4.1. 6-9 Lennard-Jones
The “6-9” Lennard-Jones nonbond interaction is:
where i and j are the two non-bonded atoms (or united atoms) that are interacting, \(\varepsilon_{ij}\) is the potential well depth of the interaction and \(r_{ij,min}\) is the interatomic distance at which the energy is minimized (i.e. where the interatomic force is zero).
Equation (1) can be equivalently written as:
where \(A_{ij}\) and \(B_{ij}\) are the interaction parameters of atoms i and j and rij is the interatomic distance.
One can easily convert one of the above formulas (1), (2) into the other:
that has the ability to read parameters in either format, to facilitate the introduction of these parameters in the frc file, in the same way as they have been introduced and published by the forcefield authors. In this way, conversion errors are avoided and the values present in the frc file can very easily be compared to those that are published in the literature.
The interaction parameters are introduced for ii interactions, i.e. interaction of identical atom types. Hetero-atomic interactions are defined thereafter as arithmetic, geometric or sixth-power averages, i.e. using a specific “combining” rule.
3.6.7.4.2. 6-12 Lennard-Jones
The “6-12” Lennard-Jones nonbond interaction is:
where i and j are the two non-bonded atoms (or united atoms) that are interacting, \(\varepsilon_{ij}\) is the potential well depth of the interaction and \(r_{ij,min}\) is the interatomic distance at which the energy is minimized (i.e. where the interatomic force is zero).
Equation (4) can be equivalently written as:
where i and j are the two non-bonded atoms (or united atoms) that are interacting, \(\varepsilon_{ij}\) is the potential well depth of the interaction and \(r_{ij,0}\) is the interatomic distance at which the energy is zero (often mentioned as \(\sigma\) ).
Equation (5) can be equivalently written as:
where \(A_{ij}\) and \(B_{ij}\) are the interaction parameters of atoms i and j and \(r_{ij}\) is the interatomic distance.
One can easily convert one of the above formulas (4), (5), (6) into the other:
that has the ability to read parameters in either format (A-B or \(\varepsilon-r_{min}\) or \(\varepsilon-r_{0}\), to facilitate the introduction of these parameters in the frc file, in the same way as they have been introduced and published by the forcefield authors. In this way, conversion errors are avoided and the values present in the frc file can very easily be compared to those that are published in the literature.
The interaction parameters are introduced for ii interactions, i.e. interaction of identical atom types. Hetero-atomic interactions are defined thereafter as arithmetic (8), geometric (9) or sixth-power (10) averages, i.e. using a specific “combining” rule.
The type of the potential (r-eps for use of Eq. (4a) or r0-eps for use of Eq. (4b) or A-B for use of Eq. (4c)) and the combining rule are specifically defined in the beginning of the nonbond 6-12 section (arithmetic, or geometric or sixth-power), e.g.:
@combination arithmetic
@type r0-eps
The format of the 6-12 Lennard-Jones nonbond energy term section is:
#nonbond(12-6) AUA
> E = 4.0\*eps(ij) [(r0(ij)\*/r(ij))\*\*12 - (r0(ij)\*/r(ij))\*\*6]
> where r0(ij)\* = 0.5\*((r0(i)\*) + (r0(j)\*))
> and eps(ij) = sqrt(eps(i) \* eps(j))
@combination arithmetic
@type r0-eps
@units Sigma Ang
@units Epsilon K
!Ver Ref I r0 eps
!---- --- ---- --------- ---------
1.0 1 CH3-AUA 3.6072 120.15
All lines starting with “>” are comments. Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 for \(\varepsilon_{ii}\) and \({\mathring{\mathrm{A}}}\) for \(r_{ii,min}\) and \(r_{ii,0}\), kcal mol-1 \({\mathring{\mathrm{A}}}\)12 for A and kcal mol-1 \({\mathring{\mathrm{A}}}\)6 for B.
3.6.7.4.3. Arithmetic Combining Rule
3.6.7.4.4. Geometric Combining Rule
3.6.7.4.5. Sixth-Power Combining Rule
The type of the potential (r-eps for use of Eq. (1) or A-B for use of Eq. (2)) and the combining rule are specifically defined in the beginning of the nonbond 6-9 section (arithmetic, or geometric or sixth-power), e.g.:
@type r-eps
@combination sixth-power
The format of the 6-9 Lennard-Jones nonbond energy term section is:
#nonbond(9-6) pcff+ 200
> E = eps(ij) [2(r(ij)\*/r(ij))\*\*9 - 3(r(ij)\*/r(ij))\*\*6]
> where r(ij) = [(r(i)\*\*6 + r(j)\*\*6))/2]\*\*(1/6)
>
> eps(ij) = 2 sqrt(eps(i) \* eps(j)) \*
> r(i)^3 \* r(j)^3/[r(i)^6 + r(j)^6]
@type r-eps
@combination sixth-power
!Ver Ref I r eps
!--- --- ----- ------------- -----------
1.0 1 ar 3.8800 0.2000
1.0 13 Br 5.4135 0.07993
All lines starting with “>” are comments. Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 for \(\varepsilon_{ii}\) and \({\mathring{\mathrm{A}}}\) for \(\sigma_{ii}\), kcal mol-1 \({\mathring{\mathrm{A}}}\)9 for A and kcal mol-1 \({\mathring{\mathrm{A}}}\)6 for B.
3.6.7.5. Buckingham
The “Buckingham” nonbond interaction is:
where A, \(\rho\) and C are constants and \(r_{ij}\) is the interatomic distance. Equation (11) can be equivalently written as:
where A, B and C are constants and \(r_{ij}\) is the interatomic distance. One can easily convert one of the above formulas (11) and (12) into the other:
that has the ability to read parameters in either format (A-Rho-C or A-B-C), to facilitate the introduction of these parameters in the frc file, in the same way as they have been introduced and published by the forcefield authors. In this way, conversion errors are avoided and the values present in the frc file can very easily be compared to those that are published in the literature.
No combining rule is needed, as the interaction parameters (A, \(\rho\) and C) are provided in the frc file for pairs of atoms.
The type of the potential (A-Rho-C for use of Eq. (11) or A-B-C for use of Eq. (12)) is specifically defined in the beginning of the nonbond Buckingham section, e.g.:
@type A-Rho-C
The format of the Buckingham nonbond energy term section is:
#nonbond(exp-6) inorganic
> E = Aij\*exp(-r/Rhoij) - Cij/r^6
@type A-Rho-C
@units A eV
@units Rho Ang
@units C eV\*Ang^6
!Ver Ref I J A Rho C
!---- --- ---- ---- ----------- ----------- ------
1.0 1 Ag1+ O2- 962.197 0.3000 0.0
All lines starting with “>” are comments. Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 for \(A_{ii}\), \({\mathring{\mathrm{A}}}\) for \(\rho_{ij}\), \({\mathring{\mathrm{A}}}\)-1 for \(B_{ij}\), and kcal mol-1 \({\mathring{\mathrm{A}}}\)6 for C.
3.6.7.6. Shinoda-DeVane-Klein (sdk)
Mesoscale simulations with the SPICA forcefield use different functional forms for Van der Waals interactions depending on the bead types interacting. For pairs with water a 12-4 potential is used:
while for all other pairs a standard Lennard-Jones potential is applied:
To ease the use of this forcefield the sdk nonbond energy term section can be used in forcefield files:
#nonbond(sdk) SPICA
@type r0-eps
@combination explicit
> E = 27/4*eps*[(sigma/r)^9-(sigma/r)^6]
> E = 3*sqrt(3)/2*eps*[(sigma/r)^12-(sigma/r)^4] (pairs involving W beads)
@units eps kcal/mol
@units sigma Ang
!Ver Ref I J cgtype sigma eps
!---- --- ---- ---- ------ ----------- -----------
1.0 1 C2T C2T lj9_6 4.8115 0.4003
1.0 6 CM2 W lj12_4 4.1722 0.2788
Combining rules cannot be applied for the Shinoda-DeVane-Klein potential. Therefore, only “explicit” is allowed in the specification and all parameters need to be listed.
All lines starting with > are comments. Any change in the functional form will not
be reflected in the calculations in the code. The “cgtype” parameter is directly passed
on to MedeA LAMMPS. Possible values are: “lj9_6”, “lj12_4” or “lj12_6”.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 for \(\epsilon_{ij}\) and \({\mathring{\mathrm{A}}}\) for \(\sigma_{ij}\).
3.6.7.7. Electrostatic Interactions
3.6.7.7.1. Bond Increments
The bond increments section provides the bond increments that are used from certain forcefields to assign charges to atoms. The charge on an atom is calculated by:
where \(q_{i}\) is the charge on atom i and \(\delta_{ij}\) is the bond increment between an atom and the other atom that it is bonded with. The sum runs over all atoms bonded to atom i.
The format of the bond-increments section is:
#bond_increments pcff
!Ver Ref I J DeltaIJ DeltaJI
!--- --- ----- ----- ------- -------
2.1 11 Ag Ag 0.0000 0.0000
3.0 10 az oah 0.3013 -0.3013
where I and J are the atom types of the two atoms that are bonded and \(\delta_{ij}\) and \(\delta_{ji}\) are the bond increments that will be added to atoms I and J, respectively.
Tip
The charges are given in units of multiple of electron charge (1.0 is a proton).
3.6.7.7.2. Charges
The charges section provides the charges that correspond to specific atom types. This is an alternative way of using bond increments, or can be used in addition to the use of bond-increments. For example, to assign charges, bond increments can be used for covalently bonded species and charges can be used for e.g. ions.
The format of the charges section is:
#charge pcff+
! Formal charge on ions
!Ver Ref I q
!--- --- ---- -----------
3.2 45 Br -1.0
3.2 45 ca+ 2.0001
where I is the atom bearing a charge q.
Tip
The charges are given in units of multiple of electron charge (1.0 is a proton).
3.6.7.8. Bond Interactions
For bonded atoms a bond interaction term may be used, depending on the description of a forcefield. A bond is representing a chemical bond and may be rigid or flexible.
If a bond is rigid, the length of the bond is constant and equal to a pre-defined value.
If a bond is flexible, the length of the bond is allowed to fluctuate around an equilibrium value. The bond may be physically thought of as a spring that is connecting the two atoms. The length of the spring under no strain is the equilibrium length of the bond. The bond term may be described by different functional forms (depending on the forcefield) and defines how easily the bond length may be augmented or reduced or equivalently, how “stiff” the spring is.
3.6.7.8.1. Quadratic Bond
A quadratic bond energy term is:
where \(l\) is the bond length, \(l_{eq}\) is the equilibrium bond length and \(k_{2}\) is the force constant.
The format of the quadratic bond energy term section is:
#quadratic\_bond oplsaa
> E = K2 \* (R - R0)^2
!Ver Ref I J R0 K2
!---- --- ---- ---- ------- --------
1.1 6 C CA 1.4900 400.0000
The third and fourth columns of the quadratic bond section refer to atom types I and J, which correspond to the bonded atoms. The fifth and sixth columns (R0 and K2) refer to the equilibrium bond lengths and force constants of Equation (17), respectively.
All lines starting with “>” are comments and are only considered as such by . Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 \({\mathring{\mathrm{A}}}\)-2 for K2 and \({\mathring{\mathrm{A}}}\) for the equilibrium bond length.
3.6.7.8.2. Quartic Bond
A quartic bond energy term is:
where \(l\) is the bond length, \(l_{eq}\) is the equilibrium bond length and \(k_{2}\), \(k_{3}\), \(k_{4}\) are the force constants for the quadratic, cubic and quartic terms, respectively. The format of the quartic bond energy term section is:
#quartic\_bond pcff+ 200
> E = K2 \* (R - R0)^2 + K3 \* (R - R0)^3 + K4 \* (R - R0)^4
!Ver Ref I J R0 K2 K3 K4
!--- --- ----- ----- ------- -------- -------- --------
1.0 3 c n\_2 1.4432 319.1593 -586.3243 961.4143
The third and fourth columns of the quartic bond section refer to atom types I and J, which correspond to the bonded atoms. The fifth column (R0) refers to the equilibrium bond lengths of the bonded atoms. The sixth, seventh and eighth columns (K2, K3 and K4) refer to the force constants of Equation (18), respectively.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 \({\mathring{\mathrm{A}}}\)-2 for K2, kcal mol-1 \({\mathring{\mathrm{A}}}\)-3 for K3 and kcal mol-1 \({\mathring{\mathrm{A}}}\)-4 for K4 and \({\mathring{\mathrm{A}}}\) for the equilibrium bond length.
3.6.7.8.3. Morse Bond
A Morse bond energy term is:
where \(D_{b}\) is the bond dissociation energy, \(l_{ij}\) is the instantaneous bond distance, \(l_{0}\) is the equilibrium bond length and \(\alpha\) is the Morse anharmonicity parameter. The format of the Morse bond energy term section is:
#morse\_bond cvff
> E = D \* (1 - exp(-ALPHA\*(R - R0)))^2
!Ver Ref I J R0 D ALPHA
!---- --- ---- ---- ------ -------- -------
2.3 23 no o- 1.2178 140.2486 2.0000
The third and fourth columns of the Morse bond section refer to atom types I and J, which correspond to the bonded atoms. The fifth column (R0) refers to the equilibrium bond lengths of the bonded atoms. The sixth and seventh (D and ALPHA) refer to the constants of Equation (19).
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 for D, \({\mathring{\mathrm{A}}}\) for the equilibrium bond length R0 and \({\mathring{\mathrm{A}}}\)-1 for the Morse anharmonicity parameter ALPHA.
3.6.7.8.4. Rigid Bond
A rigid bond is a bond for which only the bonded atoms and the equilibrium bond length need to be specified. In the case of atoms bonded with a rigid bond, there is no stretching or shrinking of the bond length. The format of the rigid bond term section is:
#rigid\_bond EH
!Ver Ref I J R0
!---- --- ----------------- ----------------- -------
1.1 1 C-arom-TraPPE-EH C-arom-TraPPE-EH 1.4
The third and fourth columns of the rigid bond section refer to atom types I and J, which correspond to the bonded atoms. The fifth column (R0) refers to the bond lengths of the bonded atoms.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. \({\mathring{\mathrm{A}}}\) for the equilibrium bond length.
3.6.7.9. Angle-Bending Interaction
For bonded atoms an angle-bending interaction term may be used, depending on the description of a forcefield. Three atoms, participating in succeeding bonds are interacting through such a term. If an angle is rigid, the angle formed by the three successively bonded atoms is constant and equal to a pre-defined value. If an angle is flexible, the value of the angle is allowed to fluctuate around an equilibrium value. The angle-bending term may be described by different functional forms (depending on the forcefield) and defines how easily the angle may be augmented or reduced.
3.6.7.9.1. Quadratic angle-bending potential
A quadratic angle-bending energy term section is:
where \(\theta\), \(\theta_{0}\) is the reference bond angle and \(k_{2}\) is the quadratic force constant.
The format of the quadratic angle-bending term section is:
#quadratic\_angle oplsaa
> E = K2 \* (Theta - Theta0)^2
!Ver Ref I J K Theta0 K2
!---- --- ---- ---- ---- -------- -------
1.1 6 CA C O 120.4000 80.0000
The third, fourth and fifth columns of the quadratic angle-bending section refer to atom types I, J and K, which correspond to the successively bonded atoms. The sixth column (Theta0) refers to the equilibrium angle of the bonded atoms and the seventh column (K2) refers to the quadratic force constant.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. degrees for the equilibrium angle and kcal mol-1 rad-2 for the force constant K2.
3.6.7.9.2. Quadratic cosine angle-bending potential
A quadratic cosine angle-bending energy term section is:
where \(\theta\) is the bond angle, \(\theta_{0}\) is the reference bond angle and \(k_{2}\) is the quadratic cosine force constant.
The format of the quadratic cosine angle-bending term section is:
#quadratic\_cosine\_angle AUA
@type K/2
@units K2 K
@units Theta0 degree
> E = K \* (cosTheta - cosTheta0)^2
!Ver Ref I J K Theta0 K2
!---- --- ----------- ---------- ------------ -------- --------
1.0 1 C-arom-AUA O-ROR-AUA C-aliph-AUA 112.0 69000.00
The third, fourth and fifth columns of the quadratic cosine angle-bending section refer to atom types I, J and K, which correspond to the successively bonded atoms. The sixth column (Theta0) refers to the equilibrium angle of the bonded atoms and the seventh column (K2) refers to the quadratic cosine force constant.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. degrees for the equilibrium angle and kcal mol-1 for the force constant K2.
3.6.7.9.3. Quartic angle-bending potential
A quartic angle-bending energy term section is:
where \(\theta\) is the bond angle, \(\theta_{0}\) is the reference bond angle and \(k_{2}\), \(k_{3}\) and \(k_{4}\) are the quartic force constants.
The format of the quadratic angle-bending term section is:
#quartic\_angle pcff+ 200
> Delta = Theta - Theta0
> E = K2 \* Delta^2 + K3 \* Delta^3 + K4 \* Delta^4
!Ver Ref I J K Theta0 K2 K3 K4
!--- --- ----- ----- ----- -------- ------- ------ ------
1.0 7 c c ct 112.7000 58.3500 0.0000 0.0000
The third, fourth and fifth columns of the quartic angle-bending section refer to atom types I, J and K, which correspond to the successively bonded atoms. The sixth column (Theta0) refers to the equilibrium angle of the bonded atoms and the seventh, eighth and ninth columns (K2, K3, K4) refer to the quartic force constants.
All lines starting with “>” are comments. Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. degrees for the equilibrium angle and kcal mol-1 rad-2, kcal mol-1 rad-3 and kcal mol-1 rad-4 for the force constants K2, K3 and K4, respectively.
3.6.7.9.4. Rigid angle-bending
A rigid angle is an angle for which only the bonded atoms and the equilibrium angle need to be specified. In the case of a rigid angle, there is no fluctuating of the angle around its equilibrium value.
The format of the rigid angle term section is:
# rigid\_angle trappeUA
!Ver Ref I J K Theta
!---- --- ----------------- ----------------- ----------------- ------
1.1 1 CH-arom-TraPPE-UA CH-arom-TraPPE-UA CH-arom-TraPPE-UA 120.00
The third and fourth columns of the rigid angle section refer to atom types I and J, which correspond to the bonded atoms. The fifth column (R0) refers to the equilibrium angle of the bonded atoms.
Tip
If the units are not defined in the beginning of this section, the default units will be used i.e. degrees for the equilibrium angle.
Note
Rigid angles can be used in MedeA-GIBBS but not in MedeA-LAMMPS simulations.
3.6.7.10. Torsion Interaction
For bonded atoms a torsion interaction term may be used, depending on the description of a forcefield. Four atoms, participating in succeeding bonds are interacting through such a term. If a torsion angle is rigid, the angle formed by the four successively bonded atoms is constant and equal to a pre-defined value. If a torsion angle is flexible, the value of the angle is allowed to fluctuate around an equilibrium value. The torsion angle term may be described by different functional forms (depending on the forcefield) and defines how easily the torsion angle may be augmented or reduced.
3.6.7.10.1. One term (cosine) torsion
A (one term) cosine torsion energy term section is:
where I, J, K and L are the atoms that are bonded to each other (in that succeeding order) forming the torsion angle \(\phi\), \(\phi\) is the instantaneous torsion angle, \(\phi_{0}\) is the equilibrium torsion angle for this quadruplet of atoms, \(K_{\phi}\) is half the energy barrier height and \(n\) is the periodicity of the torsion.
The format of the (one term) cosine torsion energy term section is:
#torsion\_1 cvff
> E = Kphi \* [ 1 + cos(n\*Phi - Phi0) ]
!Ver Ref I J K L Kphi n Phi0
!---- --- ---- ---- ---- ---- ------- ---- --------
2.3 23 \* cp no \* 10.0000 2 180.0000
1.9 17 cp cp c cp 0.6750 4 0.0000
The third, fourth, fifth and sixth columns of the (one term) cosine torsion section refer to atom types I, J, K and L, which correspond to the successively bonded atoms. The sixth column (Kphi) refers to the energy constant, the seventh column (n) refers to the periodicity of the torsion and the eighth column (Phi0) refers to the reference torsion angle.
All lines starting with “>” are comments. Any change in the functional form will not be reflected in the calculations in the code.
Tip
Asterisks are recognized as wildcards, i.e. atoms matching any atom types, as shown in the first line in the example above. Neither or both end atom types (I and L) must be wildcards.
Tip
If the units are not defined in the beginning of this section the default units will be used, i.e. degrees for the equilibrium angle and kcal mol-1 for Kphi.
3.6.7.10.2. Three term (cosine) torsion
A (three term) cosine torsion energy term section is:
where I, J, K and L are the atoms that are bonded to each other (in that succeeding order) forming the torsion angle \(\phi\), \(\phi\) is the instantaneous torsion angle, \(\phi_{1}\), \(\phi_{2}\), and \(\phi_{3}\) are the reference torsion angles for this quadruplet of atoms and \(V_{1}\), \(V_{2}\), \(V_{3}\), are the energy barrier heights.
The format of the (three term) cosine torsion energy term section is:
#torsion\_3 pcff
> E = SUM(n=1,3) { V(n) \* [ 1 + cos(n\*Phi - Phi0(n)) ] }
!Ver Ref I J K L V(1) Phi1(0) V(2) Phi2(0) V(3) Phi3(0)
!--- --- ---- ---- ---- ---- ------ ------- ------- ------- ------ -------
3.0 10 oah az oah hoa 0.2821 0.0 -0.0644 0.0 0.0752 0.0
The third, fourth, fifth and sixth columns of the (one term) cosine torsion section refer to atom types I, J, K and L, which correspond to the successively bonded atoms. Columns 7, 9 and 11 (V(1), V(2) and V(3)) refer to the energy constants, while columns 8, 10 and 12 (Phi(1), Phi(2) and Phi(3)) refer to the equilibrium torsion angles. The equilibrium torsion angles are usually 0 or 180 degrees.
All lines starting with “>” are comments. Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section the default units will be used, i.e. degrees for the equilibrium angle and kcal mol-1 for the energy constants.
3.6.7.10.3. OPLS torsion
An OPLS torsion energy term section is:
where I, J, K and L are the atoms that are bonded to each other (in that succeeding order) forming the torsion angle \(\phi\), \(\phi\) is the instantaneous torsion angle, \(\phi_{1}\), \(\phi_{2}\), \(\phi_{3}\) and \(\phi_{4}\) are the reference torsion angles for this quadruplet of atoms and \(V_{1}\), \(V_{2}\), \(V_{3}\) and \(V_{4}\) are the energy barrier heights.
The format of the OPLS torsion energy term section is:
#torsion\_opls oplsaa
> E = SUM(n=1,4) { [V(n)/2] \* [ 1 - ((-1)^n)cos(n\*Phi + Phi0(n)) ] }
> with '1-4' interactions scaled by 0.5
@units V kcal/mol
@units Phi degree
!Ver Ref I J K L V1 Phi0 V2 Phi0 V3 Phi0 V4 Phi0
!---- --- --- --- --- --- ------ ----- ------- ----- ------- ----- ------- -----
1.1 6 O C CA CA 0.0000 0.0 2.1000 0.0 0.0000 0.0 0.0000 0.0
The third, fourth, fifth and sixth columns of the OPLS torsion section refer to atom types I, J, K and L, which correspond to the successively bonded atoms. Columns 7, 9, 11 and 13 (V1, V2, V3 and V4) refer to the energy constants, while columns 8, 10, 12 and 14 (Phi1, Phi2, Phi3 and Phi4) refer to the reference torsion angles.
All lines starting with “>” are comments. Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section the default units will be used, i.e. degrees for the equilibrium torsion angle and kcal mol-1 for the energy constants.
3.6.7.10.4. TraPPE Torsion
A TraPPE torsion energy term section is:
where I, J, K and L are the atoms that are bonded to each other (in that succeeding order) forming the torsion angle \(\phi\), \(\phi\) is the instantaneous torsion angle and \(C_{0}\), \(C_{1}\), \(C_{2}\) and \(C_{3}\) are the energy barrier heights.
The format of the TraPPE torsion energy term section is:
#torsion\_trappe trappeUA-flexible
> E = C0 + C1 \* [1 + cos(phi)] + C2 \* [1 - cos(2\*phi)] + C3 \* [1 + cos(3\*phi)]
> with '1-4' electrostatic interactions scaled by 0.5 and van der Waals ignored
@units C0 K
@units C1 K
@units C2 K
@units C3 K
@units Phi degree
!Ver Ref I J K L C0 C1 C2 C3
!--- --- ------------- ------------ -------------- -------------- ----- ----- ----- ------
1.1 7 CH-TraPPE-UA C-TraPPE-UA CH2-TraPPE-UA CH2-TraPPE-UA 0.00 0.00 0.00 461.29
The third, fourth, fifth and sixth columns of the TraPPE torsion section refer to atom types I, J, K and L, which correspond to the successively bonded atoms. Columns 7, 8, 9 and 10 (C1, C2, C3 and C4) refer to the energy constants.
All lines starting with > are comments.
Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section the default units will be used, i.e. kcal mol-1 for the energy constants.
3.6.7.10.5. AUA torsion
An AUA torsion energy term section is:
where I, J, K and L are the atoms that are bonded to each other (in that succeeding order) forming the torsion angle \(\phi\), \(\phi\) is the instantaneous torsion angle and \(A_{0}\), \(A_{1}\), \(A_{2}\), \(A_{3}\), \(A_{4}\), \(A_{5}\), \(A_{6}\), \(A_{7}\), and \(A_{8}\) are the energy barrier heights.
The format of the AUA torsion energy term section is:
#torsion\_aua AUA
> E = A(n) \* (cos(phi))^n
@lj14 0.0
@el14 0.0
@units A0 K
@units A1 K
@units A2 K
@units A3 K
@units A4 K
@units A5 K
@units A6 K
@units A7 K
@units A8 K
!Ver Ref I J K L A0 A1 A2 A3 A4 A5 A6 A7 A8
!---- --- ------------ ------------- --------- -------- ------ ------ ----- ------- --- --- --- --- ---
1.0 1 CH-aliph-AUA CHx-aliph-AUA O-OH-AUA H-OH-AUA 339.41 353.97 58.34 -751.72 0.0 0.0 0.0 0.0 0.0
The third, fourth, fifth and sixth columns of the AUA torsion section refer to atom types I, J, K and L, which correspond to the successively bonded atoms. Columns 7, 8, 9 and 10 (A0, A1, A2, A3, A4, A5, A6, A7 and A8) refer to the energy constants.
All lines starting with > are comments.
Any change in the functional form will not be reflected in the calculations in the code.
Tip
If the units are not defined in the beginning of this section the default units will be used, i.e. K for the energy constants.
3.6.7.10.6. Rigid Torsion
A rigid torsion angle is an angle for which only the bonded atoms and the equilibrium torsion angle need to be specified. In the case of a rigid torsion angle, there is no fluctuating of the torsion angle.
The format of the rigid torsion angle term section is:
#rigid\_angle AUA
!Ver Ref I J K Theta0
!---- --- ---- ---- ---- --------
1.0 1 CH-arom-AUA CH-arom-AUA CH-arom-AUA 120.00
Tip
If the units are not defined in the beginning of this section the default units will be used, i.e. degrees for the equilibrium torsion angle.
3.6.7.11. Cross Interactions
3.6.7.11.1. Quadratic Bond-bond interaction
A quadratic bond-bond interaction potential, or cross-term, section is:
where \(r\) and \(r_{0}\) are the instantaneous bond length and the equilibrium bond length of the first bond respectively, \(r'\) and \(r'_{0}\) are the instantaneous bond length and the equilibrium bond length of the second bond respectively and \(k_{bb'}\) is the interaction force constant.
The format of the quadratic bond-bond interaction energy term section is:
#bond-bond pcff+ 200
> E = K(b,b') * (R - R0) * (R' - R0')
!Ver Ref I J K K(b,b')
!--- --- ----- ----- ----- --------
1.0 6 o1= c2= o1= 275.4350
The third, fourth and fifth columns of the quadratic bond-bond interaction energy term section refer to atom types I, J and K, with J being the central atom of the angle formed from the two bonds. Atoms I and J are the atom types of the first bond and atoms J and K are the atoms of the second bond. The equilibrium bond lengths \(r_{0}\) and \(r'_{0}\) are taken from the appropriate bond section of the parameter file.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. \({\mathring{\mathrm{A}}}\) for the bonds and kcal mol-1 \({\mathring{\mathrm{A}}}\)-2 for the \(k_{bb'}\) constant.
3.6.7.11.2. Quadratic bond-angle interaction
A quadratic bond-angle interaction potential, or cross-term, section is:
where \(r\) and \(r_{0}\) are the instantaneous bond length and equilibrium bond length respectively, \(\theta\) and \(\theta_{0}\) are the instantaneous bond angle and equilibrium bond angle respectively, and \(k_{b\theta}\) is the quadratic force constant.
The format of the quadratic bond-angle interaction energy term section is:
#bond-angle pcff
> E = K * (R - R0) * (Theta - Theta0)
!Ver Ref I J K K(b,theta) K(b',theta)
!--- --- ----- ----- ----- ---------- -----------
3.0 10 oah az oah 21.5223
3.0 10 oah az oas 5.5077 45.3577
The third, fourth and fifth columns of the quadratic bond-angle interaction energy term section refer to atom types I, J and K, with J being the central atom of the angle. If the atoms I and K are the same, only one force constant is provided in the sixth column. If the atoms I and K are different, then an additional force constant is provided, in the seventh column. The values of the equilibrium bond length \(r_{0}\) and equilibrium angle \(\theta_{0}\) are taken from the appropriate bond and angle term sections of the file.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. \({\mathring{\mathrm{A}}}\) for the bonds, degrees for the angles and kcal mol-1 \({\mathring{\mathrm{A}}}\)-1 rad-1 for the quadratic force constant.
3.6.7.11.3. Angle-angle interaction
An angle-angle interaction potential, i.e. the interaction between atoms that form two angles which share a bond, is:
where \(V\) is the force constant, \(\theta\) and \(\theta_{0}\) are the instantaneous and equilibrium values for the first bond angle and \(\theta'\) and \(\theta'_{0}\) are the instantaneous and equilibrium values for the second bond angle.
The format of the angle-angle interaction energy term section is:
#angle-angle pcff
> E = K * (Theta - Theta0) * (Theta' - Theta0')
!Ver Ref I J K L K(theta,theta')
!--- --- ----- ----- ----- ----- ---------------
3.0 10 oah az oah oah 11.3873
where I, J, K and L are the atom types of the atoms involved in the two angles (\(\theta\) and \(\theta'\)). Atoms I, J and K define the first angle (\(\theta\)) and atoms K, J and L define the second angle (\(\theta'\)), with the bond between atoms J and K being the common bond of the two angles. The values of the equilibrium angles (\(\theta_{0}\) and \(\theta'_{0}\)) are obtained from the appropriate angle sections.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 rad-2 for the quadratic force constant.
3.6.7.11.4. Bond-Torsion Interaction
3.6.7.11.4.1. End Bond-Torsion
An interaction potential between a three term (cosine) torsion and the end bonds of the torsion is:
where \(r\) is the instantaneous bond length, \(r_{0}\) is the equilibrium bond length, \(V_{1}\)/\(V_{2}\)/\(V_{3}\) are the force constants and \(\phi\) is the instantaneous torsion angle.
The format of the end bond-torsion interaction energy term is:
#end_bond-torsion_3 pcff
> E = (R - R0) * SUM { V(n) * cos[n*phi] }
! LEFT RIGHT
! ----------------------------- -----------------------------
!Ver Ref I J K L F(1) F(2) F(3) F(1) F(2) F(3)
!--- --- ----- ----- ----- ----- ------- ------- ------- ------- ------- -------
1.0 1 c c c c -0.0732 0.0000 0.0000
1.0 1 c c c c= -0.6028 0.0000 0.7675 1.0356 0.0000 0.0506
where I, J, K and L are the atom types involved in the torsion. \(V_{1}\), \(V_{2}\) and \(V_{3}\) (noted in the table as LEFT \(F(1)\), \(F(2)\) and \(F(3)\)) are the force constants for the interaction of the bond between atoms I and J with the torsion between atoms I, J, K and L, \(V'_{1}\), \(V'_{2}\) and \(V'_{3}\) (noted in the table as RIGHT \(F(1)\), \(F(2)\) and \(F(3)\)) are the force constants for the interaction of the bond between atoms I and J with the torsion between atoms I, J, K and L. If the bonds between atoms I and J and between K and L are identical, i.e. atom types I and L are the same and atom types J and K are also the same, then only one set of force constants should be provided, as shown in the example above. The equilibrium bond lengths are obtained from the appropriate bond section.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 \({\mathring{\mathrm{A}}}\)-1 for the force constants, \({\mathring{\mathrm{A}}}\) for the bond length and degrees for the torsion angle.
3.6.7.11.4.2. Middle Bond-Torsion
An interaction potential between a three term (cosine) torsion and the end bonds of the torsion is:
where \(r\) is the instantaneous bond length, \(r_{0}\) is the equilibrium bond length, \(V_{1}\)/\(V_{2}\)/\(V_{3}\) are the force constants and \(\phi\) is the instantaneous torsion angle.
The format of the middle bond-torsion interaction energy term is:
#middle_bond-torsion_3 pcff
> E = (R - R0) *
> { F(1) * cos(phi) + F(2) * cos(2 * phi) + F(3) * cos(3 * phi) }
!Ver Ref I J K L F(1) F(2) F(3)
!--- --- ----- ----- ----- ----- ------- ------- -------
1.0 1 c c c c -17.7870 -7.1877 0.0000
where I, J, K and L are the atom types involved in the torsion, \(V_{1}\)/\(V_{2}\) and \(V_{3}\) (noted in the table as \(F(1)\), \(F(2)\) and \(F(3)\)) are the force constants for the interaction of the bond between atoms J and K with the torsion between atoms I, J, K and L. The equilibrium bond lengths are obtained from the appropriate bond section.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 \({\mathring{\mathrm{A}}}\)-1 for the force constants, \({\mathring{\mathrm{A}}}\) for the bond length and degrees for the torsion angle.
3.6.7.11.4.3. Angle-torsion interaction
An interaction potential between an angle and a torsion defined by the same atoms is:
where \(\theta\) is the instantaneous bond angle and \(\theta_{0}\) is the equilibrium angle, \(V_{1}\)/\(V_{2}\)/\(V_{3}\) are the force constants and \(\phi\) is the instantaneous torsion angle.
The format of the angle-torsion interaction energy term is:
#angle-torsion_3 pcff
> E = (Theta - Theta0) *
> { F(1) * cos(phi) + F(2) * cos(2 * phi) + F(3) * cos(3 * phi) }
! LEFT RIGHT
! ----------------------------- -----------------------------
!Ver Ref I J K L F(1) F(2) F(3) F(1) F(2) F(3)
!--- --- ----- ----- ----- ----- ------- ------- ------- ------- ------- -------
1.0 1 c c c c 0.3886 -0.3139 0.1389
1.0 1 c c c c- 16.6010 0.1267 3.1777 -0.7732 2.4204 -1.5184
where I, J, K and L are the atom types involved in the torsion, \(V_{1}\)/\(V_{2}\) and \(V_{3}\) (noted in the table as LEFT \(F(1)\), \(F(2)\) and \(F(3)\)) are the force constants for the interaction of the bond between atoms I, J and K with the torsion between atoms I, J, K and L, math:V’_{1}/\(V'_{2}\) and \(V'_{3}\) (noted in the table as RIGHT \(F(1)\), \(F(2)\) and \(F(3)\)) are the force constants for the interaction of the bond between atoms J, K and L with the torsion between atoms I, J, K and L. The equilibrium angles are obtained from the appropriate angle sections.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 rad-1 for the energy constants and degrees for the angles.
3.6.7.12. Out-of-plane, Improper Torsion Interaction
3.6.7.12.1. Wilson out-of-plane interaction
This section supplies the parameters used for the potential when the out-of-plane coordinate is defined according to the angle between one bond from the central atom and the plane defined by the other two bonds. An out-of-plane potential is usually applied to planar groups containing an sp2 central atom bonded to three other atoms. Examples are amide nitrogens, amide carbons, and the carbon atoms in a benzene ring. The out-of-plane potential acts to keep the central atom in the plane defined by the other three atoms. The functional form is:
where \(V\) is the force constant and \(\chi\) is the instantaneous Wilson out-of-plane angle.
The format of the out-of-plane potential (Wilson definition) section is:
#wilson_out_of_plane pcff
> E = K * (Chi - Chi0)^2
!Ver Ref I J K L KChi Chi0
!--- --- ----- ----- ----- ----- -------- ---------
1.0 1 nr c+ nr nr 54.4060 0.0000
where I, J, K and L are the atom types involved in the out-of-plane term, J being the central atom, and Chi0 the instantaneous out-of-plane angle. This term is asymmetric with respect to the outer atoms (I and L) but is made symmetric by summing over the three different out-of-planes defined by a trigonal center.
Tip
If the units are not defined in the beginning of this section, the default units will be used, i.e. kcal mol-1 rad-2 for the force constant and degrees for the angle.
3.6.7.13. Embedded Atom Method
3.6.7.13.1. Introduction and functional form
As its name implies, the embedded atom method accounts for the behavior of an atom placed in a defined electron density. The method therefore captures a significant portion of the physical reality of metallic bonding. Related to the effective medium theory of Norskov and Lang [1], the embedded atom method (EAM) was developed by Daw and Baskes [2]. This approach represents the total energy of the system as two additive terms, a pairwise sum of interactions between atoms, and a term representing the electron density of each atomic site, as shown in Equation (35) below.
E is the total energy of the system, i and j indicate the unique pairs of atoms within the N atoms of the system, \(r_{ij}\) is their interatomic separation, \(V_{ij}(r_{ij})\) is a pairwise potential, and \(F_{i}(\rho_{i})\) is the embedding function for atom i which depends on the electron density, \(\rho_{i}\), experienced by that atom:
To evaluate a given atom’s embedding function, one needs to compute the electron density at the position of atom i. This is obtained by a superposition of “atomic densities”, which are described by a density function, \(\phi_{j}(r)\), as shown in Equation (36).
A typical EAM embedding function, illustrated using the Zr EAM forcefield of Mendelev and Ackland [3].
The embedding function, \(F_{i}(\rho_{i})\), provides an essential degree of freedom in the description of metallic bonding. If this term were linear with respect to varying density, the overall energetic description would be equivalent to a standard two body representation. However, the curvature of the embedding term with varying electron density provides an account of the effects of many body interactions. A common form of the embedding function for an EAM forcefield is shown in Figure 4. Here increasing electron density yields progressively more negative embedding energies, until a minimum value is attained beyond which increasing electron density yields less favorable system energies. Figures 5 and 6 provide views of the density function, employed to compute the electron density at a given site (Figure 5), and the interaction function (Figure 6) which is reminiscent of a typical two-body interaction function.
Computing energies and forces based on Equations (35) and (36) can be achieved rapidly as each of the terms are functions of interatomic separation and such separations and their derivatives with respect to atomic coordinates can be rapidly evaluated. In practice, to avoid restricting the form of the functions employed, and to promote calculation efficiency, numerically splined tabulated look-up tables are employed in most EAM calculations for the necessary functions. The resulting forcefield files are therefore large numerical tables. For individual elements three such tables are required, representing the pairwise function, the embedding function, and the density function.
Handling alloy systems requires provision for interaction functions describing the pairwise interaction of each element, in addition to embedding, and density functions. For an n-component alloy there will be n(n+1)/2 pairwise interaction functions, n-embedding functions, and, associated density functions. The determination of these functions is challenging, and is complicated by the fact that the creation of a description suitable for a single element provides little information for the behavior of that element in an alloy or compound. Consequently, EAM forcefields are typically developed for specific systems and the description of a given element cannot trivially be combined directly with the description for another element, as assumptions about the two density functions, for example, may not be compatible.
Despite such specificity, the merit of EAM forcefields is their ability to rapidly and accurately describe the bonding of metallic systems. EAM forcefields allow the simulation of:
- Structures - for example atomic configurations in the vicinity of grain boundaries
- Energies - for example relative polymorph energies and defect energies
- Diffusivity - for example through the use of mean squared displacements of sets of atoms in molecular dynamics trajectories
- Thermal expansivity - for example employing constant pressure simulations as a function of temperature to predict the response of a lattice to a temperature ramp
- Melting of metals and thermodynamic properties of the liquid state
3.6.7.13.2. EAM forms included in the MedeA Environment
The MedeA environment supports standard ‘Finnis-Sinclair’ format EAM forcefield files, with extensions to permit detailed referencing of the source of the particular EAM description and atom type assignment. Such a file contains named sections expressing atom types, any atom equivalences, the standard Finnis-Sinclair format EAM function tables, information for partial charge assignment, and template information to assign forcefield atom types based on rules concerning topology and element type. This overall format is the standard employed by all MedeA environment forcefields.
In addition to standard Finnis-Sinclair EAM forcefields, the environment also supports the EAM parameterization described by Zhou and co-workers [4]. Here, mixing rules have been implicitly included in the design of the forcefield, and any combination of the elements: Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe, Mo, Ta, W, Mg, Co, Ti, or Zr may be handled. It is likely that this generality results in diminished accuracy in some circumstances. However, when the effects of alloy formation are of interest, in the creation of layered metallic structures, for example, this description is highly effective.
Finally, the environment also supports an EAM parameterization from Bonny et al. [5] where the pair potential is defined in terms of analytical spline functions. In this specific functional form, the density function is defined as the product of a Thomas-Fermi screening function and a cutoff:
with:
where S, \(\beta\), \(r_{c}\), and \(h\) are parameters. The format of the density function term section is:
#Bonny\_atomic\_density eam
> Phi(r) = S\*(exp(-beta\*r)/r)\*x\*\*4/(1+x\*\*4)
> x(r) = (r-rc)/h
@type S
@units S eV
@type rc
@units rc Ang
!Ver Ref I S beta rc h
!---- --- ---- ---------- ----- --- ---
1.0 1 Fe 27.8689586 2.0 4.0 0.25
All lines starting with “>” are comments and are only considered as such by . Any change in the functional form will not be reflected in the calculations in the code. The embedding function is parameterized as:
The format of the embedding function term section is:
#Bonny\_embedding\_function eam
> F(rho) = A\*sqrt(rho) + B\*rho + C\*rho\*\*2 + D\*rho\*\*4
!Ver Ref I A B C D
!--- --- ---- ----------- ----------- ---------- ------------
1.0 1 Fe -8.66624513 9.41375492 -3.23721354 0.348448677
Finally, the interaction function or pair potential is defined as a linear combination of piecewise cubic splines and of the density function:
Here, \(r_{k}\) are the knots and \(a_{k}\) the fitting parameters of the \(N_{p}\) cubic spline functions. There can be a maximum of 12 splines for each i-j pair, hence the parameters range from \(\left\{ r_{k},\ a_{k} \right\}_{k = 1}\) to \(\left\{ r_{k},\ a_{k} \right\}_{k = 12}\) . \(\Theta\left( r_{k} - r \right)\) is the Heaviside unit step function, and \(\delta_{\text{xy}}\) is the Kronecker delta. K is a fitting parameter. The format of the interaction function term section is:
#Bonny\_eam\_pair eam
> V\_XY(r) = sum\_k (ak\*(rk-r)\*\*3) \* Theta(rk-r) - K\*Phi\_Y(r)\*delta\_XY
> where
> Theta(rk-r) is the Heaviside unit step function
> Phi\_Y(r) is the electron density function of atom type Y
> delta\_XY is the Kronecker delta
@units r1 Ang
@units r2 Ang
@units r3 Ang
@units r4 Ang
@units r5 Ang
@units r6 Ang
@units r7 Ang
@units r8 Ang
@units r9 Ang
@units r10 Ang
@units r11 Ang
@units r12 Ang
@units a1 eV/Ang^3
@units a2 eV/Ang^3
@units a3 eV/Ang^3
@units a4 eV/Ang^3
@units a5 eV/Ang^3
@units a6 eV/Ang^3
@units a7 eV/Ang^3
@units a8 eV/Ang^3
@units a9 eV/Ang^3
@units a10 eV/Ang^3
@units a11 eV/Ang^3
@units a12 eV/Ang^3
!Ver Ref I J r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12 a1 a2 a3 a4 a5 a6 a7
a8 a9 a10 a11 a12 K
!--- --- --- --- ------------ ------------ ------------ ------------
------------ ------------ ------------ ------------ ------------
------------ ------------ ------------ ------------ ------------
------------ ------------ ------------ ------------ ------------
------------ ------------ ------------ ------------ ------------
-----------
1.0 1 Fe Fe 2.5 2.8 3.1 3.4 3.7 4.0 2.4 2.0 0.0 0.0 0.0 0.0 2.660526
3.56676207 -1.2729965 2.10027686 -0.901005963 0.393902864 12.0 100.0 0.0
0.0 0.0 0.0 18.8275098
In this case, line wrapping makes it difficult to identify the columns; using a text editor which does not wrap lines is a good solution if one needs to edit the frc file.
The parameterization of Bonny et al. [6] for Fe-Ni-Cr is readily available in the environment.
3.6.8. Forcefields with MedeA-LAMMPS and MedeA-GIBBS
| Interaction | MedeA LAMMPS | MedeA GIBBS | |
|---|---|---|---|
| Nonbond | LJ6-9 |  |
|
| LJ6-12 | |
|
|
| Buckingham | |
||
| SDK | |
|
|
| Bond-increments | |
|
|
| Charges | |
|
|
| mie | |
|
|
| Bond | Quadratic | |
|
| Quartic | |
|
|
| Morse | |
||
| Rigid | |
||
| Angle | Quadratic | |
|
| Quadratic cosine | |
|
|
| Quartic | |
|
|
| Rigid | |
||
| Torsion | One-term cosine | |
|
| Three-term cosine | |
|
|
| OPLS | |
|
|
| TraPPE | |
|
|
| AUA | |
||
| Rigid | |
||
| Mesoscale | Martini | |
|
| SPICA | |
|
|
| EAM | Finnis-Sinclair | |
|
| Zhou | |
|
|
| Bonny | |
|
|
| MEAM | |
||
| 3-body | Tersoff | |
|
| Stillinger-Weber | |
||
| Reactive | Streitz-Mintmire | |
|
| ReaxFF | |
||
| COMB3 | |
||
| Machine Learned | SNAP | |
|
| NNP | |
||
| [1] | Norskov, J. K., & Lang, N. D. (1980). Effective-medium theory of chemical binding: Application to chemisorption. Physical Review B, 21 (6), 2131. |
| [2] | Daw, M. S., & Baskes, M. I. (1984). Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Physical Review B, 29 (12), 6443. |
| [3] | Mendelev, M. I., & Ackland, G. J. (2007). Development of an interatomic potential for the simulation of phase transformations in zirconium. Philosophical Magazine Letters, 87 (5), 349-359. |
| [4] | Zhou, X. W., Johnson, R. A., & Wadley, H. N. G. (2004). Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers. Physical Review B, 69 (14), 144113. |
| [5] | Bonny, G., Terentyev, D., Pasianot, R. C., Ponc\({\acute{e}}\) , S., & Bakaev, A. (2011). Interatomic potential to study plasticity in stainless steels: the FeNiCr model alloy. Modelling and simulation in materials science and engineering, 19 (8), 085008. |
| [6] | Bonny, G., Terentyev, D., Pasianot, R. C., Ponc\({\acute{e}}\) , S., & Bakaev, A. (2011). Interatomic potential to study plasticity in stainless steels: the FeNiCr model alloy. Modelling and simulation in materials science and engineering, 19 (8), 085008. |
| download: | pdf |
|---|