Skip to content

Commit

Permalink
review qeq docs and update the math display error (#139)
Browse files Browse the repository at this point in the history
Co-authored-by: WangXinyan940 <[email protected]>
  • Loading branch information
taipinghu and WangXinyan940 authored Nov 8, 2023
1 parent f2c6975 commit e241696
Showing 1 changed file with 48 additions and 74 deletions.
122 changes: 48 additions & 74 deletions docs/user_guide/4.3ADMPQeqForce.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,30 +7,27 @@ 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:

$$
```math
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.
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_{j>i}\frac{q_i q_j}{r_{ij}}\\
$$

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:
```math
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):
* For example, for EQeq model (`JPCL 3, 2506`), the damping kernel is ($`K`$ is the dielectric constant):

```math
E_{sr}\left(q_i\right) = \sum_{i<j} f(r_{ij})\frac{q_i q_j}{r_{ij}} \\
Expand All @@ -39,83 +36,69 @@ J_{ij} = \sqrt{J_i J_j} \\
```
* Another typical example is the Tang-Toennies damping:

$$
```math
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)):

$$
```math
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:
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:
```math
l_{ij} = \sqrt{l_i^2 + l_j^2}
$$

Then $\eta_i \rightarrow \infty$ is equivalent to $l_i \rightarrow 0$.
```
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):

$$
```math
f(r) = \frac{r}{\left[r^3 + (1/\gamma_{ij})^3\right]^{1/3}} - 1
$$
```




Typically, these damping functions must satisfy:

$$
```math
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.
```
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:

$$
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:
```math
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`):

$$
```math
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:

$$
```math
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]}
$$

```math
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)):

$$
```math
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
Expand All @@ -125,46 +108,37 @@ 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

$$
```math
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:

$$
```
And we need to find the proper value for both $`\psi_A`$ and $`q_i`$ to find the stationary point, such that:
```math
\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$

$$
This constraint sets the electrostatic potential of a group of atoms to a predefined value $`\psi_A`$
```math
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:

$$
```
We need to find the proper value for $`q_i`$ to find the stationary point, such that:
```math
\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.
```
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:

$$
```math
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
Expand Down

0 comments on commit e241696

Please sign in to comment.