Skip to content

4.1classical.md

Latest commit

 

History

History
233 lines (159 loc) · 11.7 KB

4.1classical.md

File metadata and controls

233 lines (159 loc) · 11.7 KB

HarmonicBondJaxForce

1. Theory

The chemical bond is represented by a harmonic potential:

$$ E = \frac{1}{2}k(b-b_0)^2 $$

where $k$ is the force constant, $b$ is the distance betweeen two particles that forming a bond and $b_0$ is the equilibrium bond length. Note that in some other MD softwares, the potential form adopts a slight different form: $E=k(b-b_0)^2$. Users should check which form to use and multiply (or divide) the force constant by 2.

2. Frontend

The way to specify a harmonic bond in DMFF is the same as the way doing it in OpenMM:

<HarmonicBondForce>
 <Bond class1="C" class2="C" length="0.1525" k="259408.0" mask="true"/>
 <Bond class1="C" class2="CA" length="0.1409" k="392459.2" mask="true"/>
 <Bond class1="C" class2="CB" length="0.1419" k="374049.6" mask="true"/>
 ...
</HarmonicBondForce>

Every <Bond> tag defines a rule for creating harmonic bond interactions between atoms. Each tag may identify the atoms either by type (using the attributes type1 and type2) or by class (using the attributes class1 and class2). length is the equilibrium bond length in $\mathrm{nm}$, and k is the force constant in $\mathrm{kJ/mol/nm^2}$.

When the tag has an attribute named mask and it's value set to true, this means the parameter is not trainable. Such information will be passed to ParamSet.mask (the corresponding mask value will be 0.0 if not trainable).

HarmonicAngleJaxForce

1. Theory

The angle is represented by a harmonic potential:

$$ E = \frac{1}{2}k(\theta - \theta_0)^2 $$

where $k$ is the force constant, $\theta$ is the angle between three particles and $\theta_0$ is the equilibrium angle value. Similiar to HarmonicBondJaxForce, the parameters in some other MD softwares are defined to follow the potential form: $E=k(\theta-\theta_0)^2$. Be aware to adjust the parameters properly when applying the DMFF parameters to other softwares.

2. Frontend

The way to specify a harmonic angle in DMFF is the same as the way doing it in OpenMM:

<HarmonicAngleForce>
 <Angle class1="C" class2="C" class3="O" angle="2.094" k="669.44"/>
 <Angle class1="C" class2="C" class3="OH" angle="2.094" k="669.44"/>
 <Angle class1="CA" class2="C" class3="CA" angle="2.094" k="527.184" mask="true"/>
 ...
</HarmonicAngleForce>

Every <Angle> tag defines a rule for creating harmonic angle interactions between triplets of atoms. Each tag may identify the atoms either by type (using the attributes type1, type2, type3) or by class (using the attributes class1, class2, class3). The force field identifies every set of three atoms in the system where the first is bonded to the second, and the second to the third. angle is the equilibrium angle in radians, and k is the spring constant in $\mathrm{kJ/mol/radian^2}$.

When the tag has an attribute named mask and it's value set to true, this means the parameter is not trainable. Such information will be passed to ParamSet.mask (the corresponding mask value will be 0.0 if not trainable).

PeriodicTorsionJaxForce

1. Theory

The torsion is represented by a truncated periodic Fourier series:

$$ E = \sum_{n=0}^{6} k_n(1 + \cos(n\phi-\phi_{0n})) $$

where $\phi$ is the dihedral angle formed by four particles, $n$ is the periodicity, $\phi_{0n}$ is the phase offset $k_{n}$ is the force constant. To perserve the symmetry, $\phi_{0n}$ usually adopts a value of $0$ (for $n=1,3,5$) or $\pi$ ($n=2,4,6$), and it is recommened to follow these definitions and not to optimize them in force field development.

2. Frontend

The way to specify a periodic torsion in DMFF is the same as the way doing it in OpenMM:

<PeriodicTorsionForce>
 <Proper class1="HC" class2="CT" class3="CT" class4="CT" periodicity1="3" phase1="0.0" k1="0.66944"/>
 <Proper class1="HC" class2="CT" class3="CT" class4="HC" periodicity1="3" phase1="0.0" k1="0.6276"/>
 ...
 <Improper class1="N" class2="C" class3="CT" class4="O" periodicity1="2" phase1="3.14159265359" k1="4.6024" mask="true"/>
 <Improper class1="N" class2="C" class3="CT" class4="H" periodicity1="2" phase1="3.14159265359" k1="4.6024" mask="true"/>
 ...
</PeriodicTorsionForce>

Every child tag defines a rule for creating periodic torsion interactions between sets of four atoms. Each tag may identify the atoms either by type (using the attributes type1, type2, ...) or by class (using the attributes class1, class2, ...).

The force field recognizes two different types of torsions: Proper and Improper. A proper torsion involves four atoms that are bonded in sequence: 1 to 2, 2 to 3, and 3 to 4. An improper torsion involves a central atom and three others that are bonded to it: atoms 2, 3, and 4 are all bonded to atom 1. periodicity1 is the periodicity of the torsion, phase1 is the phase offset in radians, and k1 is the force constant in kJ/mol. To add a second periodicity, just add three more attributes: periodicity2, phase2, and k2. The maxium periodicity supported in DMFF is 6, which is different from OpenMM.

You can also use wildcards when defining torsions. To do this, simply leave the type or class name for an atom empty. That will cause it to match any atom:

<Proper class1="" class2="OS" class3="P" class4="" periodicity1="3" phase1="0.0" k1="1.046"/>

When the tag has an attribute named mask and it's value set to true, this means the parameter is not trainable. Such information will be passed to ParamSet.mask (the corresponding mask value will be 0.0 if not trainable).

LennardJonesForce

1. Theory

The Lennard-Jones intearction between two particles follows the potential form:

$$ E = 4\epsilon\left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^{6}}{r^{6}}\right) $$

where $r$ is the distance between two particles, $\epsilon$ is the depth of the potential wall and $\sigma$ defines the distance where the interaction energy is zero.

The pairwise parameter $\sigma$ and $\epsilon$ are determined from the parameters of the individual particles using the Lorentz-Berthelot combining rule:

$$\sigma=\frac{\sigma_1+\sigma_2}{2}$$

$$\epsilon = \sqrt{\epsilon_1\epsilon_2} $$

2. Frontend

To specify Lennard-Jones interactions, include a tag that looks like this:

<LennardJonesForce lj14scale="0.5">
  <Atom sigma="0.247679397463798523" epsilon="0.0869048" type="gaff2-hc"/>
  <Atom sigma="0.233928874135017395" epsilon="0.0873315" type="gaff2-h1"/>
  ...
</LennardJonesForce>

The <LennardJones> tag an attribute lj14scale that specify the scale factors between pairs of atoms separated by three bonds. Each <Atom> tag specifies the nonbonded parameters for one atom type (specified with the type attribute) or atom class (specified with the class attribute). sigma is in nm, and epsilon is in $\mathrm{kJ/mol}$.

CoulombNoCutoffForce

1. Theory

The form of the Coulomb interaction between each pair of particles depends on the NonbondedMethod in use. For NoCutoff, it is given by

$$ E = \frac{1}{4\pi\epsilon_0}\frac{q_1q_2}{r} $$

where $q_1$ and $q_2$ are the charges of the two particles, and $r$ is the distance between them. $\epsilon_0$ is the dielectric constant of vacuum.

2. Frontend

To specify Coulomb interactions, include a tag that looks like this:

<CoulombForce coulomb14scale="0.83333333">
</CoulombForce>

The attribute coulomb14scale specifies the scale factors between pairs of atoms separated by three bonds. Notice that the atomic charges are not specified in this tag because different force fields may adopt different approaches to assign atomic charges. Currently, DMFF supports two approaches to assign atomic charges: template-based and AM1-BCC.

  • In force fields like AMBER, the atomic charges are defined in a template-based manner. Such force fields are usually for protein or nucleic acids that have well-defined building blocks. Hence, the atomic charges can be specified in the <Residues> node:
<Residues>
    <Residue name="ALA">
      <Atom charge="-0.4157" name="N" type="protein-N"/>
      <Atom charge="0.2719" name="H" type="protein-H"/>
      <Atom charge="0.0337" name="CA" type="protein-CX"/>
      <Atom charge="0.0823" name="HA" type="protein-H1"/>
      <Atom charge="-0.1825" name="CB" type="protein-CT"/>
      <Atom charge="0.0603" name="HB1" type="protein-HC"/>
      <Atom charge="0.0603" name="HB2" type="protein-HC"/>
      <Atom charge="0.0603" name="HB3" type="protein-HC"/>
      <Atom charge="0.5973" name="C" type="protein-C"/>
      <Atom charge="-0.5679" name="O" type="protein-O"/>
      <Bond atomName1="N" atomName2="H"/>
      <Bond atomName1="N" atomName2="CA"/>
      <Bond atomName1="CA" atomName2="HA"/>
      <Bond atomName1="CA" atomName2="CB"/>
      <Bond atomName1="CA" atomName2="C"/>
      <Bond atomName1="CB" atomName2="HB1"/>
      <Bond atomName1="CB" atomName2="HB2"/>
      <Bond atomName1="CB" atomName2="HB3"/>
      <Bond atomName1="C" atomName2="O"/>
      <ExternalBond atomName="N"/>
      <ExternalBond atomName="C"/>
    </Residue>
    ...
</Residues>
  • For arbitrary organic molecules, AM1-BCC approach appears to be an efficient way to derive high-quality charges. The atomic charges are calculated by applying bond charge corrections (BCC) to Mulliken charges under AM1 level of theory. In DMFF, one can define BCC parameters with SMIRKS patterns and call sqm to calculate AM1 charges:
<ForceField>
<Operators>
    <AM1ChargeOperator>
        <SQM path="sqm"/>
    </AM1ChargeOperator>
</Operators>
<CoulombForce coulomb14scale="0.83333333">
    <UseBondChargeCorrection name="MOL"/>
    <BondChargeCorrection bcc="0.0000" smirks="[#6X4:1]-[#6X4:2]"/>
    <BondChargeCorrection bcc="0.0278" smirks="[#6X4:1]-[#6X3:2]"/>
    ...
</CoulombForce>

NonbondedForce

1. Theory

The NonbondedForce is a summary of CoulombForce and LennardJonesForce for consistency with OpenMM. With NonbondedForce, the force field library of OpenMM can be fluentely used in DMFF. The form of the Lennard-Jones and Coulomb interaction between each pair of particles depends on the NonbondedMethod in use.

2. Frontend

To specify NonbondedForce interactions, include a tag that looks like this:

  <NonbondedForce coulomb14scale="0.8333333333333334" lj14scale="0.5">
    <UseAttributeFromResidue name="charge"/>
    <Atom epsilon="0.359824" sigma="0.3399669508423535" type="protein-C"/>
    <Atom epsilon="0.359824" sigma="0.3399669508423535" type="protein-CA"/>
  </NonbondedForce>

The attribute coulomb14scale and lj14scale specifies the scale factors between pairs of atoms separated by three bonds. The atomic charges are defined in a template-based manner since the UseAttributeFromResidue node is added. If node UseAttributeFromResidue does not exist, the atomic charges should be specified with Atom node, such as:

  <NonbondedForce coulomb14scale="0.8333333333333334" lj14scale="0.5">
    <Atom epsilon="0.359824" charge="0.131452" sigma="0.3399669508423535" type="protein-C"/>
    <Atom epsilon="0.359824" charge="-0.051223" sigma="0.3399669508423535" type="protein-CA"/>
  </NonbondedForce>
⚡