Update 4.3ADMPQeqForce.md

docs/user_guide/4.3ADMPQeqForce.md

@@ -7,6 +7,7 @@ ADMPQeqForce provides a support to coulombic energy calculation for constant pot
 ### 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)
 $$
@@ -18,9 +19,11 @@ 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}}
+E_{pc}\left(q_i\right) =\sum_{j>i}\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
@@ -30,16 +33,18 @@ The second term is the short-range damping term, which can be different for diff
 * 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}} \\
+E_{sr}\left(q_i\right) = \sum_{j>i} 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}
 $$
+
 * The third common option is the damping corresponding to Guassian smeared charge ([ewald.pdf (gitlab.com)](https://gitlab.com/ampere2/metalwalls/-/raw/release/doc/theory/ewald-summation/ewald.pdf)):
 
 $$
@@ -50,33 +55,40 @@ $$
 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} \\
@@ -86,19 +98,25 @@ $$
 #### 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]}
+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)](https://gitlab.com/ampere2/metalwalls/-/raw/release/doc/theory/ewald-summation/ewald.pdf)):
+
 $$
 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
@@ -108,33 +126,42 @@ Once we have the coulombic energy, we can add different constraints to construct
 #### 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}}