Skip to content

Latest commit

 

History

History
320 lines (238 loc) · 12.9 KB

4.3ADMPQeqForce.md

File metadata and controls

320 lines (238 loc) · 12.9 KB
Raw

ADMPQeqForce

1. Theory

ADMPQeqForce provides a support to coulombic energy calculation for constant potential model and constant charge model. Net charges on all atoms were equilibrated at setted constraint first, then charge related energys were carried out next.

1.1 General Energy Equations

First, assume the coulombic energy of the total system is:

$$E_{coul}(q_i) = E_{pc}\left(q_i\right) + E_{sr}(q_i) + E_{corr} + E_{on-site}({q_i}; \chi_i, J_i)$$

where $q_i$ is the charge of atom $i$, $\chi_i$ and $J_i$ are the electronegtivitity and hardness of atom $i$ respectively.

The total coulombic energy is composed of the followling four parts.

1.1.1 Point Charge Part

Here, the first term $E_{pc}$ is electrostatic energy of point-charge systems, which is typically done in Ewald/PME method. For isolated molecule/clusters, it can also be computed using cutoff scheme:

$$E_{pc}\left(q_i\right) = \sum_{i< j}\frac{q_i q_j}{r_{ij}}$$

Note that we currently only consider Ewad3D (so called Ewald3DC method), instead of rigorous Ewald2D.

1.1.2 Short Range Damping Part

The second term is the short-range damping term, which can be different for different models and users can choose it by a key word DampMod or edit their own damping mode.

  • For example, for EQeq model (JPCL 3, 2506), the damping kernel is ($K$ is the dielectric constant):
$$E_{sr}\left(q_i\right) = \sum_{i< j} f(r_{ij})\frac{q_i q_j}{r_{ij}} \\\ f(r) = - \exp\left[-\left(\frac{J_{ij}r}{K}\right)^2\right]\left(1-\frac{J_{ij}r}{K}+\left(\frac{J_{ij}r}{K}\right)^2\right) \\\ J_{ij} = \sqrt{J_i J_j}$$
  • Another typical example is the Tang-Toennies damping:
$$f_{TT}^1\left(r\right) = -\exp\left(-B_{ij}r\right)\left(1 + B_{ij}r\right) \\\ B_{ij} = \sqrt{B_i B_j}$$ $$f_{gg}(r) = -\text{erfc}\left(\eta_{ij}r\right) \\\ \eta_{ij} = \frac{\eta_i \eta_j}{\sqrt{\eta_i^2 + \eta_j^2}}$$

Note that one can easily recover the point charge formula by setting $\eta_i \rightarrow \infty$.

For this combination rule, it is perhaps it is more code-friendly if we define: $l_i = 1/\eta_i$, then equivalent to the above equation we have:

$$l_{ij} = \sqrt{l_i^2 + l_j^2}$$

Then $\eta_i \rightarrow \infty$ is equivalent to $l_i \rightarrow 0$.

  • The short-range damping model Prof. Zhu Tong uses (double check with him for more details): $$f(r) = \frac{r}{\left[r^3 + (1/\gamma_{ij})^3\right]^{1/3}} - 1$$

Typically, these damping functions must satisfy:

$$f(r) \rightarrow 0 \text{ when } r\rightarrow \infty \\\ f(r) \rightarrow -1 \text{ when } r\rightarrow 0$$

And are related to some nonlinear atomic parameters ($J_i, B_i, \eta_i$ etc.) subject to optimization.

1.1.3 Non-Neutral and Slab Boundary Corrections

The $E_{corr}$ is other correction term, including nonneutral background charge corrections and dipole (or more generally, slab boundary) correction, which are quite commonly used in slab electrode simulations (see JCP 147 126101). The nonneutral correction is:

$$E_{corr}^Q = - \frac{\pi Q_{tot}^2}{2\kappa^2 V} \\$$

And there are also dipole boundary correction, which, in the case of nonneutral box, is (JCP 131 094107):

$$E_{corr}^{dip} = \frac{2\pi}{V}\left(M_z^2-Q_{tot}\sum_i {q_i z_i^2} - Q_{tot}^2\frac{L_z^2}{12}\right) \\\ M_z = \sum_i {q_i z_i} \\\ Q_{tot} = \sum_i q_i$$

1.1.4 On-Site Energy

The on-site energy is the self energy of each site, which is usually in the following form:

$$E_{on-site} = \sum_i {\left[\chi_i q_i + J_iq_i^2\right]}$$

Note that quite often, it can be written as:

$$E_{on-site} = \sum_i {\left[\chi_i (q_i-q_i^*) + J_i(q_i-q_i^*)^2\right]}$$

Which is equivalent to the above equation subject to a constant shift.

Also, for Gaussian smear charge, there is often a self-interaction term (ewald.pdf (gitlab.com)):

$$E_{self}^g = \sum_i {\frac{q_i^2 \eta_i}{\sqrt{2\pi}}}$$

which can be also merged into Eqn 11. Therefore there is no need to design a separate term for this (?).

1.2 Constraints

Once we have the coulombic energy, we can add different constraints to construct the optimization target function for charges. Depending on different scenarios, we usually have two types of constraints.

1.2.1 Constant Charge (ConstQ) Constraint

This constraint ensures that the total charge of a group of atoms is constant

$$E = E_{coul}(q_i) - \sum_A {\psi_A \left({\sum_{i\in A} q_i - Q_A}\right)}$$

And we need to find the proper value for both $\psi_A$ and $q_i$ to find the stationary point, such that:

$$\begin{cases} \frac{\partial E}{\partial q_i} = 0 \\\ \frac{\partial E}{\partial \psi_A} = 0 \\\ \end{cases}$$

Note that in this condition, the constraint term is always zero at stationary point, so it does not enter into force evaluations. Bottom line is, at stationary point, we can simply neglect its existence while calculating forces (NOTE that it is not true at points other than stationary points).

1.2.2 Constant Potential (ConstP) Constraint

This constraint sets the electrostatic potential of a group of atoms to a predefined value $\psi_A$

$$E = E_{coul}(q_i) - \sum_A {\psi_A \left({\sum_{i\in A} q_i }\right)}$$

We need to find the proper value for $q_i$ to find the stationary point, such that:

$$\frac{\partial E}{\partial q_i} = 0$$

Note that for this condition, the constraint term ($\psi q_i$) can be viewed as a physical energy, which represents the energy cost to remove the charge from a constant potential reservoir. And this term SHOULD enter into force evaluations.

Point is, once we find the stationary point $q^\star$, the force evaluation is quite easy in either scenario with Feynman-Hellman theorem applicable:

$$F_i = -\frac{\partial E}{\partial r_i} - \sum_j {\frac{\partial E}{\partial q_j^\star} \frac{\partial q_j^\star}{\partial r_i}} \\\ = \left.-\frac{\partial E}{\partial r_i}\right|_{q=q^{\star}}$$

References

  1. J. Phys. Chem. Lett. 2012, 3, 17, 2506–2511; https://doi.org/10.1021/jz3008485
  2. J. Chem. Phys. 147, 126101 (2017); https://doi.org/10.1063/1.4998320
  3. J. Chem. Phys. 131, 094107 (2009); https://doi.org/10.1063/1.3216473
  4. J. Chem. Phys. 153, 174704 (2020); https://doi.org/10.1063/5.0028232

2. Frontend

A typical xml file to define an ADMPQeqForce frontend is:

 <ADMPQeqForce coulomb14scale="1.0" DampMod="3" >
   <Atom type="1" chi="1.39198761e+02"  J="-6.51683561e-01" eta="0.08360402" />
   <Atom type="2" chi="1.47411181e+02"  J="-3.42450470e-02" eta="0.08360402" />
   <Atom type="3" chi="7.58342711e+01"  J="-5.62083783e-01" eta="0.08360402" />
   <Atom type="4" chi="1.40533063e+02"  J="-6.33544439e-01" eta="0.08360402" />
   <Atom type="5" chi="1.37525951e+02"  J="-6.63409032e-01" eta="0.08360402" />
 </ADMPQeqForce>

The important tags in the frontend includes:

  • coulomb14scale: the scaling factors for the 1-3 coulombic interactions between bonded atoms.
  • DampMod: the tag to set Short Range Damping model.

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

  • type: the label for the atomtype.
  • chi: the electronegtivitity for each atom atomtype.
  • J: the hardness for each atom atomtype
  • eta: the length damping parameter used in calculating corresponding short range damping correction. The unit should be in nm

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

hamilt = Hamiltonian('forcefield.xml')
pot = hamilt.createPotential(dmfftop, nonbondedCutoff=rc*unit.nanometer, nonbondedMethod=app.PME,
                            ethresh=1e-3, neutral=True, slab=True, constQ=True,
                            const_list=const_list, const_vals=const_val,
                            has_aux=True
                             )

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

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

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

  • neutral: bool type, if the total net charge of each residue is not zero, we should set this tag as True to add nonneutral correction.

  • slab: bool type, if the system is with slab boundary, we should set this tag as True to add Slab Boundary Corrections.

  • constQ: bool type, if the constraint condition is constant charge, we should set this tag as True, else False for constant potential condition.

  • const_list: the atomic numbers of all residues with different constraints.

  • const_vals: the corresponding constraint valu on each residue.

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

  • has_aux: bool type, whether to intruduce auxilliary parameters, equilibrated atomic charges, etc.

3. ADMPQeqGenerator Doc

ATTRIBUTES:

The following attributes are accessible via the ADMPQeqGenerator:

  • coulomb14scale, damp_mod, ethresh, has_aux: all these variables introduced above are accessible through generator.

  • atom_keys: character strings to identify the atomic type.

  • key_type: atomic keys to set attributes, class or type.

  • coeff_method, fourier_spacing : parameters to setup ewald parameters used for calculating electrostatic energy of point-charge .

  • mscales_coul: mscale for PME.

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

METHODS

The following methods are accessible via ADMPQeqGenerator:

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

4. ADMPQeqForce Doc

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

ATTRIBUTES:

  • damp_mod: bool type, the tag to set Short Range Damping model
  • neutrial_flag: bool type, the tag to add nonneutral correction
  • slab_flag: bool type, the tag to add Slab Boundary Corrections
  • pbc_flag: bool type, the tag to choose cutoff option or nocutoff
  • constQ: bool type, the tag to choose constant charge constraint or constant potential constraint
  • const_list: int array,the atomic numbers of all residues with different constraints
  • const_value: float array, the corresponding constraint valu on each residue
  • init_q: float array, intial atomic charge
  • init_lagmt: float array, inital lagrange multipliers on each constQ constraint
  • kappa: float, the range-separation parameter $\kappa$ used in PME calculation
  • e_constraint: constraint function that made a virtual energy
  • e_sr: short range damping energy function
  • e_site: on site energy function
  • coul_energy: point charge interaction energy function

METHOS

  • E_full:

    This function constructs the charge related coulombic energy plus constraint energy.

Inputs:
    q:
        Na: the inital atomic charge 
    lagmt:
        Nr: the inital lagrange multipliers on each constQ constraint 
    chi:
        Na: the atomic electronegtivity 
    J:
        Na: the atomic hardness
    pos:
        N * 3: the positions matrix
    box:
        3 * 3: box
    pairs:
        Np: interacting pair indices within cutoff distance
    eta:
        N: the atomic length damping parameter
    ds:
        Np: topology distance of pairs
    buffer_scales:
        N * 3: buffer scales on pairs
    mscales:
        (6,): mscale for PME
Outputs:
    energy: charge related coulombic energy plus constraint energy

  • E_grads:

Same as E_full, only return the concatenation of gradients of funtion E_ful on the first two paramters.

  • get_energy:

    This contains charge equilbrant process, and retruns total energy after charge equilibration.

Input:
    positions:
        Na * 3: positions
    box:
        3 * 3: box
    pairs:
        Np * 3: interacting pair indices within cutoff topology 
    mScales:
        (6,): mscale for PME
    eta:
        Na: the atomic length damping parameter
    chi:
        Na: the atomic electronegtivity 
    J:
        Na: the atomic hardness
    aux:
        dict[str,array]: auxilliary parameters, including initial atomic charges and lagmt
Output:
    energy: total energy after charge equilibration
    aux: auxilliary parameters, including equilibrated atomic charges and lagmt