Skip to content

Latest commit

 

History

History
445 lines (308 loc) · 21.6 KB

4.2ADMPPmeForce.md

File metadata and controls

445 lines (308 loc) · 21.6 KB
Raw

ADMPPmeForce

1. Theory

ADMPmeForce provides a support to multipolar polarizable electrostatic energy calculation.

1.1 Multipole Expansion

The electrostatic interaction between two atoms can be described using multipole expansion, in which the electron cloud of an atom can be expanded as a series of multipole moments including charges, dipoles, quadrupoles, and octupoles etc. If only the charges (zero-moment) are considered, it is reduced to the point charge model in classical force fields:

$$ V=\sum_{ij} \frac{q_i q_j}{r_{ij}} $$

where $q_i$ is the charge of atom $i$.

More complex (and supposedly more accurate) force field can be obtained by including more multipoles with higher orders. Some force fields, such as MPID, goes as high as octupoles. Currently in DMFF, we support up to quadrupoles:

$$ V=\sum_{tu} Q_t^A T^{AB}_{tu} Q_u^B $$

where $Q_t^A$ represents the $t$-component of the multipole moment of atom $A$. Note there are two (equivalent) ways to define multipole moments: cartesian and spherical harmonics. Cartesian representation is over-complete but with a simpler definition, while spherical harmonics are easier to use in real calculations. In the user API (the xml file), we use cartesian representation, in consistent with the AMOEBA and the MPID plugins in OpenMM. However, spherical harmonics are always used in the backend computation kernel, and we assume all components are arranged in the following order:

$$0, 10, 1c, 1s, 20, 21c, 21s, 22c, 22s, ...$$

The $T_{tu}^{AB}$ represents the interaction tensor between multipoles. The mathematical expression of these tensors can be found in the appendix F of Ref 1. The user can also find the conversion rule between different representations in Ref 1.

1.2 Coordinate System for Multipoles

Different to charges, the definition of multipole moments depends on local coordinate system. The exact value of the moment tensor will be rotated in accord to different coordinate systems. There are three types of frames involved in DMFF, each used in a different scenario:

  • Global frame: coordinate system binds to the simulation box. It is same for all the atoms. We use this frame to calculate the charge density structure factor $S(\vec{k})$ in reciprocal space.
  • Local frame: this frame is defined differently on each atom, determined by the positions of its peripheral atoms. Normally, atomic multipole moments are most stable in the local frame, so it is the most suitable frame for force field input. In DMFF API, the local frames are defined in the same way as the AMOEBA plugin in OpenMM. The details can found in the following references:
  • Quasi internal frame, aka. QI frame: this frame is defined for each pair of interaction sites, in which the z-axis is pointing from one site to another. In this frame, the real-space interaction tensor ($T_{tu}^{AB}$) can be greatly simplified due to symmetry. We thus use this frame in the real space calculation of PME.

1.3 Polarization

ADMPPmeForce supports polarizable force fields, in which the dipole moment of the atom can respond to the change of the external electric field. In practice, each atom has not only permanent multipoles $Q_t$, but also induced dipoles $U_{ind}$. The induced dipole-induced dipole and induced dipole-permanent multipole interactions needs to be damped at short-range to avoid polarization catastrophe. In DMFF, we use the Thole damping scheme identical to MPID (ref 5), which introduces a damping width ($a_i$) for each atom $i$. The damping function is then computed and applied to the corresponding interaction tensor. Taking $U_{ind}$-permanent charge interaction as an example, the definition of damping function is:

$$ \displaylines{ 1-\left(1+a u+\frac{1}{2} a^{2} u^{2}\right) e^{-a u} \\ a=a_i + a_j \\ u=r_{ij}/\left(\alpha_i \alpha_j\right)^{1/6} } $$

Other damping functions between multipole moments can be found in Ref 5, table I.

It is noted that the atomic damping parameter $a=a_i+a_j$ is only effective on topological neighboring pairs (with $pscale = 0$), while a default value of $a_{default}$ is set for all other pairs. In DMFF, the atomic $a_i$ is specified via the xml API, while $a_{default}$ is controlled by the DEFAULT_THOLE_WIDTH variable, which is set to 5.0 by default. This variable can be changed in the following way:

from dmff.admp import pme
pme.DEFAULT_THOLE_WIDTH = 1.0

We solve $U_{ind}$ by minimizing the electrostatic energy:

$$ V=V_{perm-perm}+V_{perm-ind}+V_{ind-ind} $$

The last two terms are related to $U_{ind}$. Without introducing the nonlinear polarization terms (e.g., some force fields introduce $U^4$ to avoid polarization catastrophe), the last two terms are quadratic to $U_{ind}$:

$$ V_{perm-ind}+V_{ind-ind}=U^TKU-FU $$

where the off-diagonal term of $K$ matrix is induced-induced dipole interaction, the diagonal term is formation energy of the induced dipoles ($\sum_i \frac{U_i^2}{2\alpha_i}$); the $F$ matrix represents permanent multipole - induced dipole interaction. We use the gradient descent method to optimize energy to get $U_{ind}$.

In the current version, we temporarily assume that the polarizability is spherically symmetric, thus the polarizability $\alpha_i$ is a scalar, not a tensor. Thus the inputs (polarizabilityXX, polarizabilityYY, polarizabilityZZ) in the xml API is averaged internally. In future, it is relatively simple to relax this restriction: simply change the reciprocal of the polarizability to the inverse of the matrix when calculating the diagonal terms of the $K$ matrix.

1.4 Dispersion Interaction

In ADMP, we assume that the following expansion is used for the long-range dispersion interaction:

$$ V_{disp}=\sum_{ij}-\frac{C_{ij}^6}{r_{ij}^6}-\frac{C_{ij}^8}{r_{ij}^8}-\frac{C_{ij}^{10}}{r_{ij}^{10}}-... $$

where the dispersion coefficients are determined by the following combination rule:

$$ C^n_{ij}=\sqrt{C_i^n C_j^n} $$

Note that the dispersion terms should be consecutive even powers according to the perturbation theory, so the odd dispersion terms are not supported in ADMP.

In ADMP, this long-range dispersion is computed using PME (vida infra), just as electrostatic terms.

In the classical module, dispersions are treated as short-range interactions using standard cutoff scheme.

1.5 Long-Range Interaction with PME

Taking charge-charge interaction as example, the interaction decays in the form of $O(\frac{1}{r})$, and its energy does not converge with the increase of cutoff distance. The multipole electrostatics and dispersion interactions also converge slow with respect to cutoff distance. We therefore use Particle Meshed Ewald(PME) method to calculate these interactions.

In PME, the interaction tensor is splitted into the short-range part and the long-range part, which are tackled in real space and reciprocal space, respectively. For example, the Coulomb interaction is decomposed as:

$$ \frac{1}{r}=\frac{erfc(\kappa r)}{r}+\frac{erf(\kappa r)}{r} $$

The first term is a short-range term, which can be calculated directly by using a simple distance cutoff in real space. The second term is a long-range term, which needs to be calculated in reciprocal space by fast Fourier transform(FFT). The total energy of charge-charge interaction is computed as:

$$ \displaylines{ E_{real} = \sum_{ij}\frac{erfc(\kappa r_{ij})}{r_{ij}} \\ E_{recip} = \sum_{\vec{k}\neq 0} {\frac{2\pi}{Vk^2}\exp\left[\frac{k^2}{4\kappa^2}\right]\left|S(\vec{k})\right|^2}\frac{1}{\left|\theta(\vec{k})\right|^2} \\ E_{self} = -\frac{\kappa}{\sqrt{\pi}}\sum_i {q_i^2} \\ E = E_{real}+E_{recip}+E_{self} } $$

As for multipolar PME and dispersion PME, the users and developers are referred to Ref 2, 3 for mathematical details.

The key parameters in PME include:

  • $\kappa$: controls the separation of the long-range and the short-range. The larger $\kappa$ is, the faster the real space energy decays, the smaller the cutoff distance can be used in the real space, and more difficult it is to converge the reciprocal energy and the larger $K_{max}$ it needs;

  • $r_{c}$: cutoff distance in real space;

  • $K_{max}$: controls the number of maximum k-points in all three dimensions

In DMFF, we determine these parameters in the same way as in OpenMM:

$$ \displaylines{ \kappa=\sqrt{-\log (2 \delta)} / r_{c} \\ K_{max}=\frac{2 \kappa d}{3 d^{1 / 5}} } $$

where the user needs to specify the cutoff distance $r_c$ when building the neighbor list, the width of the box in each dimension $d$ (determined from the input box matrix), and the energy accuracy $\delta$.

In the current version, the dispersion PME calculator uses the same parameters as in electrostatic PME.

References

  1. Anthony's book
  2. The Multipolar Ewald paper in JCTC: J. Chem. Theory Comput. 2015, 11, 2, 436–450
  3. The dispersion Ewald/PME
  4. Frenkel & Smit book
  5. MPID Reference

2. Frontend

A typical xml file to define an ADMPPmeForce frontend is:

 <ADMPPmeForce lmax="2"
     mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="1.00" mScale16="1.00"
     pScale12="0.00" pScale13="0.00" pScale14="0.00" pScale15="1.00" pScale16="1.00"
     dScale12="1.00" dScale13="1.00" dScale14="1.00" dScale15="1.00" dScale16="1.00">
   <Atom type="380" kz="381" kx="-381"
             c0="-0.803721"
             dX="0.0" dY="0.0"  dZ="-0.00784325"
             qXX="0.000366476" qXY="0.0" qYY="-0.000381799" qXZ="0.0" qYZ="0.0" qZZ="1.53231e-05"
             />
   <Atom type="381" kz="380" kx="381"
             c0="0.401876"
             dX="-0.00121713" dY="0.0"  dZ="-0.00095895"
             qXX="6.7161e-06" qXY="0.0" qYY="-3.37874e-05" qXZ="1.25905e-05" qYZ="0.0" qZZ="2.70713e-05"
             />
   <Polarize type="380" polarizabilityXX="1.1249e-03" polarizabilityYY="1.1249e-03" polarizabilityZZ="1.1249e-03" thole="0.33"/>
   <Polarize type="381" polarizabilityXX="2.6906e-04" polarizabilityYY="2.6906e-04" polarizabilityZZ="2.6906e-04" thole="0.33"/>
 </ADMPPmeForce>

The important tags in the frontend includes:

  • lmax: the maximal rank of the multipoles defined in the force. "2" represents a truncation at quadrupole level.
  • mScales1n: the scaling factors for the nonbonding interactions between bonded atoms. mScales12 is used for 1-2 interactions, mScales13 is used for 1-3 interactions etc. These scaling factors are applied to scale the permanent multipole - permanent multipole interactions.
  • pScales1n: the scaling factors for the nonbonding interactions between bonded atoms. pScales12 is used for 1-2 interactions, pScales13 is used for 1-3 interactions etc. These scaling factors are applied to scale the permanent multipole - induced dipole interactions.
  • dScales1n: the scaling factors for the nonbonding interactions between bonded atoms. dScales12 is used for 1-2 interactions, dScales13 is used for 1-3 interactions etc. These scaling factors are applied to scale the induced dipole - induced dipole interactions.

For each atom type, there is an <Atom/> node, with the following tags:

  • type: the label for the atomtype.
  • kx, ky, kz: these are the labels for local coordinate definitions, please refer to OpenMM AMOEBA plugin for detailed explanations.
  • c0: atomic charges
  • dX, dY, dZ: atomic dipoles, defined in Cartesian, defined in e*nm.
  • qXX, qXY, ..: atomic quadrupole moment, defined in Cartesian, defined in e*nm^2

If the force field is polarizable, then there would be a <Polarize/> node for each atomtype, too, with the following tags:

  • PolarizabilityXX, PolarizabilityYY, PolarizabilityZZ: Polarizabilities of the atomtype, currently we only support isotropic polarization, so the three values are averaged internally. To avoid confusion, we suggest you to use same values for all three tags. The unit should be in nm^3.
  • thole: the thole damping parameters used in short range damping.

To create a ADMP PME potential function, one needs to do the following in python:

H = Hamiltonian('forcefield.xml')
pots = H.createPotential(pdb.topology, nonbondedCutoff=rc*unit.nanometer, nonbondedMethod=app.PME, ethresh=5e-4, step_pol=5)

Then the potential function, the parameters, and the generator can be accessed as:

efunc_pme = pots.dmff_potentials["ADMPPmeForce"]
params_pme = H.getParameters()["ADMPPmeForce"]
generator_pme = H.getGenerators()[0] # if ADMPPmeForce is the first force node in xml

A few keywords in createPotential would impact the behavior of the function:

  • steps_pol: number of iterations to converge the induced dipoles. Only when this tag is set, the resulting potential function can be jitted. Otherwise, the convergence threshold is controlled by the variable dmff.admp.settings.POL_CONV (default value 1.0). In this case, the number of iterations would not be fixed, thus the generated potential function cannot be jitted. The POL_CONV variable can be accessed as:

    from dmff.admp import settings
settings.POL_COV = 2.0

  • ethresh: the accuracy of PME, used to setup the size of the meshgrid of PME.

A typical xml file to define an ADMPDispPmeForce frontend is:

   <ADMPDispPmeForce 
       mScale12="0.00" mScale13="0.00" mScale14="0.00" mScale15="0.00" mScale16="0.00" >
       <Atom type="1" C6="1.507393e-03" C8="1.241793e-04" C10="4.890285e-06" />
       <Atom type="2" C6="1.212045e-04" C8="4.577425e-06" C10="8.708729e-08" />
       <Atom type="3" C6="7.159800e-04" C8="5.846871e-05" C10="2.282115e-06" />
       <Atom type="4" C6="1.523005e-03" C8="1.127912e-04" C10="4.005600e-06" />
       <Atom type="5" C6="1.136931e-04" C8="4.123377e-06" C10="7.495037e-08" />
   </ADMPDispPmeForce>

The important tags in the frontend includes:

  • mScales1n: the scaling factors for the nonbonding interactions between bonded atoms. mScales12 is used for 1-2 interactions, mScales13 is used for 1-3 interactions etc. These scaling factors are applied to scale the dispersion interactions.

For each atom type, there is an <Atom/> node, with the following tags:

  • type: the label for the atomtype.
  • C6, C8, C10: dispersion terms, defined in kj/mol/nm^n.

To create a DispADMP PME potential function, one needs to do the following in python:

H = Hamiltonian('forcefield.xml')
pots = H.createPotential(pdb.topology, nonbondedCutoff=rc*unit.nanometer, nonbondedMethod=app.PME, ethresh=5e-4)

Then the potential function, the parameters, and the generator can be accessed as:

efunc_pme = pots.dmff_potentials["ADMPDispPmeForce"]
params_pme = H.getParameters()["ADMPDispPmeForce"]
generator_pme = H.getGenerators()[0] # if ADMPDispPmeForce is the First force node in xml

The keyword ethresh in createPotential would impact the behavior of the function:

  • ethresh: the accuracy of PME, used to setup the size of the meshgrid of PME.

3. ADMPPmeGenerator Doc

ATTRIBUTES:

The following attributes are accessible via the ADMPPmeGenerator:

  • mScales, dScales, pScales, lmax, ethresh, step_pol, lmax: all these variables introduced above are accessible through generator

  • lpol: bool type, is the force polarizable?

  • map_atomtype, map_poltype: integer array, the typification map for atom types and polarization types. Can be used to expand the force field parameters to atomic parameters. For example:

    Q_atoms = params['ADMPPmeForce']['Q_local'][pme_generator.map_atomtype]

  • ffinfo: the force field parameters registered in the generator.

  • pme_force: the backend ADMPPmeForce object

METHODS

The following methods are accessible via ADMPPmeGenerator:

  • overwrite(paramset): overwrite the parameters registered in the generator (i.e., the ffinfo object) using the values in paramset.

ATTRIBUTES:

The following attributes are accessible via the ADMPDispPmeGenerator:

  • mScales, ethresh: all these variables introduced above are accessible through generator

  • map_atomtype: integer array, the typification map for atom types types. Can be used to expand the force field parameters to atomic parameters. For example:

    C6_atoms = params['ADMPDispPmeForce']['C6'][pme_generator.map_atomtype]

  • ffinfo: the force field parameters registered in the generator.

  • disp_pme_force: the backend ADMPDispPmeForce object

METHODS

The following methods are accessible via ADMPDispPmeGenerator:

  • overwrite(paramset): overwrite the parameters registered in the generator (i.e., the ffinfo object) using the values in paramset.

4. ADMPPmeForce Doc

The backend of the ADMP PME energy is an ADMPPmeForce object. It contains the following attributes and methods:

ATTRIBUTES:

  • axis_indices: int array, neighbor atom indices used to define the local coordinate system of each atom
  • axis_type: int array, the type of local coordinate definition (e.g., z-only, bisector, z-then-x etc.)
  • ethresh: float, accuracy for PME setups
  • kappa: float, the range-separation parameter $\kappa$ used in PME calculation
  • lconverg: bool, if the last calculation converged?
  • lmax: int, maximal multipole rank in the force
  • lpme: bool, use PME or not (should be TRUE 99.9% of time)
  • lpol: bool, polarizable or not?
  • n_atoms: number of atoms
  • n_cycle: number of iterations in the last SCF calculation
  • pme_order: order of PME interpolation, now only support 6
  • rc: real space cutoffs
  • steps_pol: number of iterations (if specified, the code will do exact this number of iterations when computing induced dipole)

METHOS

(ALERT!!!: in ADMP backend, all inputs are assumed to be in Angstrom, instead of nm, different to frontend):

  • construct_local_frames:

    This function constructs the local frames for each site

Inputs:
    positions:
        N * 3: the positions matrix
    box:
Outputs:
    local_frames:
        N*3*3, the local frames, axes arranged in rows

  • get_energy:

    • if nonpolarizable:
    This is the top-level wrapper for multipole PME, returns total PME energy

Input:
    positions:
        Na * 3: positions
    box:
        3 * 3: box
    pairs:
        Np * 3: interacting pair indices and topology distance
    Q_local:
        Na * (lmax+1)^2: harmonic multipoles of each site in local frame
    mScales:
        (Nexcl,): multipole-multipole interaction exclusion scalings: 1-2, 1-3 ...
        for permanent-permanent interactions
Output:
    energy: total PME energy
    • if polarizable
    This is the top-level wrapper for multipole PME

Input:
    positions:
        Na * 3: positions
    box:
        3 * 3: box
    pairs:
        Np * 3: interacting pair indices and topology distance
    Q_local:
        Na * (lmax+1)^2: harmonic multipoles of each site in local frame
    pol:
        (Na,) float: the polarizability of each site, unit in A**3
    tholes:
        (Na,) float: the thole damping widths for each atom, it's dimensionless, default is 8 according to MPID paper
    mScales, pScale, dScale:
        (Nexcl,): multipole-multipole interaction exclusion scalings: 1-2, 1-3 ...
        for permanent-permanent, permanent-induced, induced-induced interactions
    U_init:
        (Na * 3) float: initial guesses for induced dipoles (default is None)

Output:
    energy: total PME energy

  • get_forces

    Same as get_energy, only return forces per atom, instead of total energy, equivalent to jax.grad(get_energy)

The backend of the ADMPDisp PME energy is an ADMPDispPmeForce object. It contains the following attributes and methods:

ATTRIBUTES:

  • ethresh: float, accuracy for PME setups
  • kappa: float, the range-separation parameter $\kappa$ used in PME calculation
  • lpme: bool, use PME or not (should be TRUE 99.9% of time)
  • n_atoms: number of atoms
  • pme_order: order of PME interpolation, now only support 6
  • rc: real space cutoffs

METHOS

(ALERT!!!: in ADMPDisp backend, all inputs are assumed to be in Angstrom, instead of nm, different to frontend):

  • get_energy:

    This is the top-level wrapper for DispPME, returns total DispPME energy

Input:
    positions:
        Na * 3: positions
    box:
        3 * 3: box
    pairs:
        Np * 3: interacting pair indices and topology distance
    c_list:
      Na * (pmax-4)/2: atomic dispersion coefficients 
    mScales: 
     (Nexcl,):  dispersion interaction exclusion scalings: 1-2, 1-3 ...
    pmax:
      int array: maximal exponents (p) to compute, e.g., (6, 8 ,10)
Output:
    energy: total disp PME energy

  • get_forces

    Same as get_energy, only return forces per atom, instead of total energy, equivalent to jax.grad(get_energy)