convert pair styles dpd to exp6

This commit is contained in:
Axel Kohlmeyer
2020-02-24 15:29:57 -05:00
parent 9955c8f94b
commit 02e287bf51
60 changed files with 184 additions and 376 deletions

Binary file not shown.

Before

Width:  |  Height:  |  Size: 28 KiB

View File

@ -1,15 +0,0 @@
\documentclass[12pt]{article}
\usepackage{amsmath}
\begin{document}
\begin{align*}
E =& E_2 \sum_{i,j}e^{-k_2 r_{ij}} + E_A \sum_{\substack{i,j,k,\ell \\\in \textrm{type A}}} f(r_{ij})f(r_{k\ell}) + E_B \sum_{\substack{i,j,k,\ell \\\in \textrm{type B}}} f(r_{ij})f(r_{k\ell}) + E_C \sum_{\substack{i,j,k,\ell \\\in \textrm{type C}}} f(r_{ij})f(r_{k\ell}) \\
f(r) =& e^{-k_3 r}s(r) \\
s(r) =& \begin{cases}
1 & r<R_s \\
\displaystyle\frac{(R_f-r)^2(R_f-3R_s+2r)}{(R_f-R_s)^3} & R_s\leq r\leq R_f \\
0 & r>R_f\\
\end{cases}
\end{align*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 5.0 KiB

View File

@ -1,9 +0,0 @@
\documentstyle[12pt]{article}
\begin{document}
$$
E_{Pauli(ECP_s)}=p_1\exp\left(-\frac{p_2r^2}{p_3+s^2} \right)
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.5 KiB

View File

@ -1,8 +0,0 @@
\documentstyle[12pt]{article}
\begin{document}
$$
E_{Pauli(ECP_p)}=p_1\left( \frac{2}{p_2/s+s/p_2} \right)\left( r-p_3s\right)^2\exp \left[ -\frac{p_4\left( r-p_3s \right)^2}{p_5+s^2} \right]
$$

Binary file not shown.

Before

Width:  |  Height:  |  Size: 3.3 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
E_{KE} = \frac{\hbar^2 }{{m_{e} }}\sum\limits_i {\frac{3}{{2s_i^2 }}}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 4.1 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
E_{NN} = \frac{1}{{4\pi \varepsilon _0 }}\sum\limits_{i < j} {\frac{{Z_i Z_j }}{{R_{ij} }}}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 6.3 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
E_{Ne} = - \frac{1}{{4\pi \varepsilon _0 }}\sum\limits_{i,j} {\frac{{Z_i }}{{R_{ij} }}Erf\left( {\frac{{\sqrt 2 R_{ij} }}{{s_j }}} \right)}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 5.5 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
E_{Pauli} = \sum\limits_{\sigma _i = \sigma _j } {E\left( { \uparrow \uparrow } \right)_{ij}} + \sum\limits_{\sigma _i \ne \sigma _j } {E\left( { \uparrow \downarrow } \right)_{ij}}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 6.7 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
E_{ee} = \frac{1}{{4\pi \varepsilon _0 }}\sum\limits_{i < j} {\frac{1}{{r_{ij} }}Erf\left( {\frac{{\sqrt 2 r_{ij} }}{{\sqrt {s_i^2 + s_j^2 } }}} \right)}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.3 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
U\left(R,r,s\right) = E_{NN} \left( R \right) + E_{Ne} \left( {R,r,s} \right) + E_{ee} \left( {r,s} \right) + E_{KE} \left( {r,s} \right) + E_{PR} \left( { \uparrow \downarrow ,S} \right)
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 12 KiB

View File

@ -1,13 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
\begin{eqnarray*}
\vec{f} & = & (F^C + F^D + F^R) \hat{r_{ij}} \qquad \qquad r < r_c \\
F^C & = & A w(r) \\
F^D & = & - \gamma w^2(r) (\hat{r_{ij}} \bullet \vec{v_{ij}}) \\
F^R & = & \sigma w(r) \alpha (\Delta t)^{-1/2} \\
w(r) & = & 1 - r/r_c
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 25 KiB

View File

@ -1,12 +0,0 @@
\documentclass[12pt]{article}
\pagestyle{empty}
\begin{document}
\begin{eqnarray*}
du_{i}^{cond} & = & \kappa_{ij}(\frac{1}{\theta_{i}}-\frac{1}{\theta_{j}})\omega_{ij}^{2} + \alpha_{ij}\omega_{ij}\zeta_{ij}^{q}(\Delta{t})^{-1/2} \\
du_{i}^{mech} & = & -\frac{1}{2}\gamma_{ij}\omega_{ij}^{2}(\frac{\vec{r_{ij}}}{r_{ij}}\bullet\vec{v_{ij}})^{2} -
\frac{\sigma^{2}_{ij}}{4}(\frac{1}{m_{i}}+\frac{1}{m_{j}})\omega_{ij}^{2} -
\frac{1}{2}\sigma_{ij}\omega_{ij}(\frac{\vec{r_{ij}}}{r_{ij}}\bullet\vec{v_{ij}})\zeta_{ij}(\Delta{t})^{-1/2} \\
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 8.7 KiB

View File

@ -1,11 +0,0 @@
\documentclass[12pt]{article}
\pagestyle{empty}
\begin{document}
\begin{eqnarray*}
\alpha_{ij}^{2} & = & 2k_{B}\kappa_{ij} \\
\sigma^{2}_{ij} & = & 2\gamma_{ij}k_{B}\Theta_{ij} \\
\Theta_{ij}^{-1} & = & \frac{1}{2}(\frac{1}{\theta_{i}}+\frac{1}{\theta_{j}}) \\
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 1.6 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\pagestyle{empty}
\begin{document}
$$
\omega_{ij} = 1 - \frac{r_{ij}}{r_{c}}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 58 KiB

View File

@ -1,14 +0,0 @@
\documentclass[12pt]{article}
\usepackage{amsmath}
\usepackage{bm}
\begin{document}
\begin{eqnarray*}
E &=& \frac{1}{2} \sum_{i} \sum_{j\notin\text{layer}\,i} \phi_{ij} \\\phi_{ij} &=& f_\text{c}(x_r) \left[ e^{-\lambda(r_{ij} - z_0 )} \left[C+f(\rho_{ij})+ g(\rho_{ij}, \{\alpha_{ij}^{(m)}\}) \right]- A\left (\frac{z_0}{r_{ij}} \right)^6 \right] \\
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 6.5 KiB

View File

@ -1,10 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
E_i = F_\alpha \left(\sum_{j \neq i}\ \rho_\beta (r_{ij})\right) +
\frac{1}{2} \sum_{j \neq i} \phi_{\alpha\beta} (r_{ij})
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 6.4 KiB

View File

@ -1,11 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
$$
E_i = F_\alpha \left(\sum_{j \neq i}\
\rho_{\alpha\beta} (r_{ij})\right) +
\frac{1}{2} \sum_{j \neq i} \phi_{\alpha\beta} (r_{ij})
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 42 KiB

View File

@ -1,22 +0,0 @@
\documentclass[12pt]{article}
\usepackage{amssymb,amsmath}
\begin{document}
\begin{eqnarray*}
E & = & \sum_{j \ne i} \phi_{2}(R_{ij}, Z_{i}) + \sum_{j \ne i} \sum_{k \ne i,k > j} \phi_{3}(R_{ij}, R_{ik}, Z_{i}) \\
\phi_{2}(r, Z) & = & A\left[\left(\frac{B}{r}\right)^{\rho} - e^{-\beta Z^2}\right]exp{\left(\frac{\sigma}{r-a}\right)} \\
\phi_{3}(R_{ij}, R_{ik}, Z_i) & = & exp{\left(\frac{\gamma}{R_{ij}-a}\right)}exp{\left(\frac{\gamma}{R_{ik}-a}\right)}h(cos\theta_{ijk},Z_i) \\
Z_i & = & \sum_{m \ne i} f(R_{im}) \qquad
f(r) = \begin{cases}
1 & \quad r<c \\
\exp\left(\frac{\alpha}{1-x^{-3}}\right) & \quad c<r<a \\
0 & \quad r>a
\end{cases} \\
h(l,Z) & = & \lambda [(1-e^{-Q(Z)(l+\tau(Z))^2}) + \eta Q(Z)(l+\tau(Z))^2 ] \\
Q(Z) & = & Q_0 e^{-\mu Z} \qquad \tau(Z) = u_1 + u_2 (u_3 e^{-u_4 Z} - e^{-2u_4 Z})
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 6.7 KiB

View File

@ -1,9 +0,0 @@
\documentstyle[12pt]{article}
\begin{document}
$$
E = \frac{1}{2} \sum_{i=1}^{N} \sum_{j=i_1}^{i_N} \phi_{ij} \left(r_{ij}\right) + \sum_{i=1}^{N}E_i\left(q_i,\sigma_i\right)
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 8.3 KiB

View File

@ -1,11 +0,0 @@
\documentstyle[12pt]{article}
\begin{document}
\begin{eqnarray*}
q_i & = & \sum_{j=i_1}^{i_N} \eta_{ji}\left(r_{ij}\right) \\
\sigma_i & = & \sum_{j=i_1}^{i_N} q_j \cdot \psi_{ij} \left(r_{ij}\right) \\
E_i\left(q_i,\sigma_i\right) & = & \frac{1}{2} \cdot q_i \cdot \sigma_i
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 40 KiB

View File

@ -1,25 +0,0 @@
\documentstyle[12pt]{article}
\begin{document}
\begin{eqnarray*}
\phi_{ij}\left(r\right) = \left\{ \begin{array}{lr}
\left[\frac{E_{b,ij}\beta_{ij}}{\beta_{ij}-\alpha_{ij}}\exp\left(-\alpha_{ij} \frac{r-r_{e,ij}}{r_{e,ij}}\right)-\frac{E_{b,ij}\alpha_{ij}}{\beta_{ij}-\alpha_{ij}}\exp\left(-\beta_{ij} \frac{r-r_{e,ij}}{r_{e,ij}}\right)\right]f_c\left(r,r_{e,ij},r_{c,\phi,ij}\right),& p_{ij}=1 \\
\left[\frac{E_{b,ij}\beta_{ij}}{\beta_{ij}-\alpha_{ij}} \left(\frac{r_{e,ij}}{r}\right)^{\alpha_{ij}} -\frac{E_{b,ij}\alpha_{ij}}{\beta_{ij}-\alpha_{ij}} \left(\frac{r_{e,ij}}{r}\right)^{\beta_{ij}}\right]f_c\left(r,r_{e,ij},r_{c,\phi,ij}\right),& p_{ij}=2
\end{array}
\right.
\end{eqnarray*}
$$
\eta_{ji} = A_{\eta,ij}\left(\chi_j-\chi_i\right)f_c\left(r,r_{s,\eta,ij},r_{c,\eta,ij}\right)
$$
$$
\psi_{ij}\left(r\right) = A_{\psi,ij}\exp\left(-\zeta_{ij}r\right)f_c\left(r,r_{s,\psi,ij},r_{c,\psi,ij}\right)
$$
$$
f_{c}\left(r,r_p,r_c\right) = 0.510204 erfc\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.2 KiB

View File

@ -1,9 +0,0 @@
\documentclass[12pt]{article}
\pagestyle{empty}
\begin{document}
$$
U_{ij}(r) = \frac{\epsilon}{\alpha-6}\{6exp[\alpha(1-\frac{r_{ij}}{R_{m}})]-\alpha(\frac{R_{m}}{r_{ij}})^6\}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 16 KiB

View File

@ -1,11 +0,0 @@
\documentstyle[12pt]{article}
\pagestyle{empty}
\begin{document}
\begin{eqnarray*}
R_{m}^{3} &=& \displaystyle\sum_{a}\displaystyle\sum_{b} x_{a}x_{b}R_{m,ab}^{3} \\
\epsilon &=& \frac{1}{R_{m}^{3}}\displaystyle\sum_{a}\displaystyle\sum_{b} x_{a}x_{b}\epsilon_{ab}R_{m,ab}^{3} \\
\alpha &=& \frac{1}{\epsilon R_{m}^{3}}\displaystyle\sum_{a}\displaystyle\sum_{b} x_{a}x_{b}\alpha_{ab}\epsilon_{ab}R_{m,ab}^{3} \\
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 7.5 KiB

View File

@ -1,11 +0,0 @@
\documentstyle[12pt]{article}
\pagestyle{empty}
\begin{document}
\begin{eqnarray*}
\epsilon_{ab} &=& \sqrt{\epsilon_{a}\epsilon_{b}} \\
R_{m,ab} &=& \frac{R_{m,a}+R_{m,b}}{2} \\
\alpha_{ab} &=& \sqrt{\alpha_{a}\alpha_{b}} \\
\end{eqnarray*}
\end{document}

View File

@ -67,17 +67,25 @@ pair interaction and the thermostat for each pair of particles.
For style *dpd*\ , the force on atom I due to atom J is given as a sum
of 3 terms
.. image:: Eqs/pair_dpd.jpg
:align: center
.. math::
where Fc is a conservative force, Fd is a dissipative force, and Fr is
a random force. Rij is a unit vector in the direction Ri - Rj, Vij is
the vector difference in velocities of the two atoms = Vi - Vj, alpha
is a Gaussian random number with zero mean and unit variance, dt is
the timestep size, and w(r) is a weighting factor that varies between
0 and 1. Rc is the cutoff. Sigma is set equal to sqrt(2 Kb T gamma),
where Kb is the Boltzmann constant and T is the temperature parameter
in the pair\_style command.
\vec{f} = & (F^C + F^D + F^R) \hat{r_{ij}} \qquad \qquad r < r_c \\
F^C = & A w(r) \\
F^D = & - \gamma w^2(r) (\hat{r_{ij}} \bullet \vec{v_{ij}}) \\
F^R = & \sigma w(r) \alpha (\Delta t)^{-1/2} \\
w(r) = & 1 - r/r_c
where :math:`F^C` is a conservative force, :math:`F^D` is a dissipative
force, and :math:`F^R` is a random force. :math:`r_{ij}` is a unit
vector in the direction :math:`r_i - r_j`, :math:`V_{ij} is the vector
difference in velocities of the two atoms :math:`= \vec{v}_i -
\vec{v}_j, :math:`\alpha` is a Gaussian random number with zero mean and
unit variance, dt is the timestep size, and w(r) is a weighting factor
that varies between 0 and 1. :math:`r_c` is the cutoff. :math:`\sigma`
is set equal to :math:`\sqrt{2 k_B T \gamma}`, where :math:`k_B` is the
Boltzmann constant and T is the temperature parameter in the pair\_style
command.
For style *dpd/tstat*\ , the force on atom I due to atom J is the same
as the above equation, except that the conservative Fc term is
@ -97,7 +105,7 @@ the examples above, or in the data file or restart files read by the
commands:
* A (force units)
* gamma (force/velocity units)
* :math:`\gamma` (force/velocity units)
* cutoff (distance units)
The last coefficient is optional. If not specified, the global DPD
@ -125,7 +133,6 @@ the work of :ref:`(Afshar) <Afshar>` and :ref:`(Phillips) <Phillips>`.
includes all the components of force listed above, including the
random force.
----------

View File

@ -54,25 +54,34 @@ under isoenergetic and isoenthalpic conditions (see :ref:`(Lisal) <Lisal3>`).
For DPD simulations in general, the force on atom I due to atom J is
given as a sum of 3 terms
.. image:: Eqs/pair_dpd.jpg
:align: center
.. math::
where Fc is a conservative force, Fd is a dissipative force, and Fr is
a random force. Rij is a unit vector in the direction Ri - Rj, Vij is
the vector difference in velocities of the two atoms = Vi - Vj, alpha
is a Gaussian random number with zero mean and unit variance, dt is
the timestep size, and w(r) is a weighting factor that varies between
0 and 1. Rc is the cutoff. The weighting factor, omega\_ij, varies
between 0 and 1, and is chosen to have the following functional form:
\vec{f} = & (F^C + F^D + F^R) \hat{r_{ij}} \qquad \qquad r < r_c \\
F^C = & A w(r) \\
F^D = & - \gamma w^2(r) (\hat{r_{ij}} \bullet \vec{v_{ij}}) \\
F^R = & \sigma w(r) \alpha (\Delta t)^{-1/2} \\
w(r) = & 1 - r/r_c
.. image:: Eqs/pair_dpd_omega.jpg
:align: center
where :math:`F^C` is a conservative force, :math:`F^D` is a dissipative
force, and :math:`F^R` is a random force. :math:`r_{ij}` is a unit
vector in the direction :math:`r_i - r_j`, :math:`V_{ij} is the vector
difference in velocities of the two atoms :math:`= \vec{v}_i -
\vec{v}_j, :math:`\alpha` is a Gaussian random number with zero mean and
unit variance, dt is the timestep size, and w(r) is a weighting factor
that varies between 0 and 1. Rc is the cutoff. The weighting factor,
:math:`\omega_{ij}`, varies between 0 and 1, and is chosen to have the
following functional form:
.. math::
\omega_{ij} = 1 - \frac{r_{ij}}{r_{c}}
Note that alternative definitions of the weighting function exist, but
would have to be implemented as a separate pair style command.
For style *dpd/fdt*\ , the fluctuation-dissipation theorem defines gamma
to be set equal to sigma\*sigma/(2 T), where T is the set point
For style *dpd/fdt*\ , the fluctuation-dissipation theorem defines :math:`\gamma`
to be set equal to :math:`\sigma^2/(2 T)`, where T is the set point
temperature specified as a pair style parameter in the above examples.
The following coefficients must be defined for each pair of atoms types
via the :doc:`pair_coeff <pair_coeff>` command as in the examples above,
@ -80,33 +89,42 @@ or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>` commands:
* A (force units)
* sigma (force\*time\^(1/2) units)
* :math:`\sigma` (force\*time\^(1/2) units)
* cutoff (distance units)
The last coefficient is optional. If not specified, the global DPD
cutoff is used.
Style *dpd/fdt/energy* is used to perform DPD simulations
under isoenergetic and isoenthalpic conditions. The fluctuation-dissipation
theorem defines gamma to be set equal to sigma\*sigma/(2 dpdTheta), where
dpdTheta is the average internal temperature for the pair. The particle
internal temperature is related to the particle internal energy through
a mesoparticle equation of state (see :doc:`fix eos <fix>`). The
differential internal conductive and mechanical energies are computed
within style *dpd/fdt/energy* as:
Style *dpd/fdt/energy* is used to perform DPD simulations under
isoenergetic and isoenthalpic conditions. The fluctuation-dissipation
theorem defines :math:`\gamma` to be set equal to :math:`sigma^2/(2
\theta)`, where :math:theta` is the average internal temperature for the
pair. The particle internal temperature is related to the particle
internal energy through a mesoparticle equation of state (see :doc:`fix
eos <fix>`). The differential internal conductive and mechanical
energies are computed within style *dpd/fdt/energy* as:
.. math::
du_{i}^{cond} = & \kappa_{ij}(\frac{1}{\theta_{i}}-\frac{1}{\theta_{j}})\omega_{ij}^{2} + \alpha_{ij}\omega_{ij}\zeta_{ij}^{q}(\Delta{t})^{-1/2} \\
du_{i}^{mech} = & -\frac{1}{2}\gamma_{ij}\omega_{ij}^{2}(\frac{\vec{r_{ij}}}{r_{ij}}\bullet\vec{v_{ij}})^{2} -
\frac{\sigma^{2}_{ij}}{4}(\frac{1}{m_{i}}+\frac{1}{m_{j}})\omega_{ij}^{2} -
\frac{1}{2}\sigma_{ij}\omega_{ij}(\frac{\vec{r_{ij}}}{r_{ij}}\bullet\vec{v_{ij}})\zeta_{ij}(\Delta{t})^{-1/2}
.. image:: Eqs/pair_dpd_energy.jpg
:align: center
where
.. image:: Eqs/pair_dpd_energy_terms.jpg
:align: center
.. math::
Zeta\_ij\^q is a second Gaussian random number with zero mean and unit
\alpha_{ij}^{2} = & 2k_{B}\kappa_{ij} \\
\sigma^{2}_{ij} = & 2\gamma_{ij}k_{B}\Theta_{ij} \\
\Theta_{ij}^{-1} = & \frac{1}{2}(\frac{1}{\theta_{i}}+\frac{1}{\theta_{j}})
:math:`\zeta_ij^q` is a second Gaussian random number with zero mean and unit
variance that is used to compute the internal conductive energy. The
fluctuation-dissipation theorem defines alpha\*alpha to be set
equal to 2\*kB\*kappa, where kappa is the mesoparticle thermal
fluctuation-dissipation theorem defines :math:`alpha^2` to be set
equal to :math:2k_B\kappa`, where :math:`\kappa` is the mesoparticle thermal
conductivity parameter. The following coefficients must be defined for
each pair of atoms types via the :doc:`pair_coeff <pair_coeff>`
command as in the examples above, or in the data file or restart files
@ -114,8 +132,8 @@ read by the :doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands:
* A (force units)
* sigma (force\*time\^(1/2) units)
* kappa (energy\*temperature/time units)
* :math:`\sigma` (force\*time\^(1/2) units)
* :math:`\kappa` (energy\*temperature/time units)
* cutoff (distance units)
The last coefficient is optional. If not specified, the global DPD

View File

@ -40,8 +40,11 @@ in :ref:`(Wen) <Wen2018>`, which is based on the :ref:`(Kolmogorov) <Kolmogorov2
potential and provides an improved prediction for forces.
The total potential energy of a system is
.. image:: Eqs/pair_drip.jpg
:align: center
.. math::
E = & \frac{1}{2} \sum_{i} \sum_{j\notin\text{layer}\,i} \phi_{ij} \\
\phi_{ij} = &f_\text{c}(x_r) \left[ e^{-\lambda(r_{ij} - z_0 )} \left[C+f(\rho_{ij})+ g(\rho_{ij}, \{\alpha_{ij}^{(m)}\}) \right]- A\left (\frac{z_0}{r_{ij}} \right)^6 \right]
where the *r\^-6* term models the attractive London dispersion,
the exponential term is designed to capture the registry effect due to

View File

@ -57,8 +57,19 @@ Description
The *e3b* style computes an \"explicit three-body\" (E3B) potential for water :ref:`(Kumar 2008) <Kumar>`.
.. image:: Eqs/e3b.jpg
:align: center
.. math::
E =& E_2 \sum_{i,j}e^{-k_2 r_{ij}} + E_A \sum_{\substack{i,j,k,\ell \\
\in \textrm{type A}}} f(r_{ij})f(r_{k\ell}) + E_B \sum_{\substack{i,j,k,\ell \\
\in \textrm{type B}}} f(r_{ij})f(r_{k\ell}) + E_C \sum_{\substack{i,j,k,\ell \\
\in \textrm{type C}}} f(r_{ij})f(r_{k\ell}) \\
f(r) =& e^{-k_3 r}s(r) \\
s(r) =& \begin{cases}
1 & r<R_s \\
\displaystyle\frac{(R_f-r)^2(R_f-3R_s+2r)}{(R_f-R_s)^3} & R_s\leq r\leq R_f \\
0 & r>R_f\\
\end{cases}
This potential was developed as a water model that includes the three-body cooperativity of hydrogen bonding explicitly.
To use it in this way, it must be applied in conjunction with a conventional two-body water model, through *pair\_style hybrid/overlay*.
@ -103,7 +114,7 @@ If the neigh setting is too large, the pair style will use more memory than nece
This pair style tallies a breakdown of the total E3B potential energy into sub-categories, which can be accessed via the :doc:`compute pair <compute_pair>` command as a vector of values of length 4.
The 4 values correspond to the terms in the first equation above: the E2 term, the Ea term, the Eb term, and the Ec term.
See the examples/USER/e3b directory for a complete example script.
See the examples/USER/misc/e3b directory for a complete example script.
----------

View File

@ -102,8 +102,11 @@ Style *eam* computes pairwise interactions for metals and metal alloys
using embedded-atom method (EAM) potentials :ref:`(Daw) <Daw>`. The total
energy Ei of an atom I is given by
.. image:: Eqs/pair_eam.jpg
:align: center
.. math::
E_i = F_\alpha \left(\sum_{j \neq i}\ \rho_\beta (r_{ij})\right) +
\frac{1}{2} \sum_{j \neq i} \phi_{\alpha\beta} (r_{ij})
where F is the embedding energy which is a function of the atomic
electron density rho, phi is a pair potential interaction, and alpha
@ -371,8 +374,12 @@ alloys using a generalized form of EAM potentials due to Finnis and
Sinclair :ref:`(Finnis) <Finnis1>`. The total energy Ei of an atom I is
given by
.. image:: Eqs/pair_eam_fs.jpg
:align: center
.. math::
E_i = F_\alpha \left(\sum_{j \neq i}\
\rho_{\alpha\beta} (r_{ij})\right) +
\frac{1}{2} \sum_{j \neq i} \phi_{\alpha\beta} (r_{ij})
This has the same form as the EAM formula above, except that rho is
now a functional specific to the atomic types of both atoms I and J,

View File

@ -37,15 +37,27 @@ potentials, while *edip/multi* supports multi-element EDIP runs.
In EDIP, the energy E of a system of atoms is
.. image:: Eqs/pair_edip.jpg
:align: center
.. math::
where phi2 is a two-body term and phi3 is a three-body term. The
summations in the formula are over all neighbors J and K of atom I
within a cutoff distance = a.
Both terms depend on the local environment of atom I through its
effective coordination number defined by Z, which is unity for a
cutoff distance < c and gently goes to 0 at distance = a.
E = & \sum_{j \ne i} \phi_{2}(R_{ij}, Z_{i}) + \sum_{j \ne i} \sum_{k \ne i,k > j} \phi_{3}(R_{ij}, R_{ik}, Z_{i}) \\
\phi_{2}(r, Z) = & A\left[\left(\frac{B}{r}\right)^{\rho} - e^{-\beta Z^2}\right]exp{\left(\frac{\sigma}{r-a}\right)} \\
\phi_{3}(R_{ij}, R_{ik}, Z_i) = & exp{\left(\frac{\gamma}{R_{ij}-a}\right)}exp{\left(\frac{\gamma}{R_{ik}-a}\right)}h(cos\theta_{ijk},Z_i) \\
Z_i = & \sum_{m \ne i} f(R_{im}) \qquad
f(r) = \begin{cases}
1 & \quad r<c \\
\exp\left(\frac{\alpha}{1-x^{-3}}\right) & \quad c<r<a \\
0 & \quad r>a
\end{cases} \\
h(l,Z) = & \lambda [(1-e^{-Q(Z)(l+\tau(Z))^2}) + \eta Q(Z)(l+\tau(Z))^2 ] \\
Q(Z) = & Q_0 e^{-\mu Z} \qquad \tau(Z) = u_1 + u_2 (u_3 e^{-u_4 Z} - e^{-2u_4 Z})
where :math:`\phi_2` is a two-body term and :math:`\phi_3` is a
three-body term. The summations in the formula are over all neighbors J
and K of atom I within a cutoff distance = a. Both terms depend on the
local environment of atom I through its effective coordination number
defined by Z, which is unity for a cutoff distance < c and gently goes
to 0 at distance = a.
Only a single pair\_coeff command is used with the *edip* style which
specifies a EDIP potential file with parameters for all

View File

@ -98,25 +98,19 @@ and the quantum-derived Pauli (E\_PR) and Kinetic energy interactions
potentials between electrons (E\_KE) for a total energy expression
given as,
.. image:: Eqs/eff_energy_expression.jpg
:align: center
.. math::
U\left(R,r,s\right) = E_{NN} \left( R \right) + E_{Ne} \left( {R,r,s} \right) + E_{ee} \left( {r,s} \right) + E_{KE} \left( {r,s} \right) + E_{PR} \left( { \uparrow \downarrow ,S} \right)
The individual terms are defined as follows:
.. image:: Eqs/eff_KE.jpg
:align: center
.. math::
.. image:: Eqs/eff_NN.jpg
:align: center
.. image:: Eqs/eff_Ne.jpg
:align: center
.. image:: Eqs/eff_ee.jpg
:align: center
.. image:: Eqs/eff_Pauli.jpg
:align: center
E_{KE} = & \frac{\hbar^2 }{{m_{e} }}\sum\limits_i {\frac{3}{{2s_i^2 }}} \\
E_{NN} = & \frac{1}{{4\pi \varepsilon _0 }}\sum\limits_{i < j} {\frac{{Z_i Z_j }}{{R_{ij} }}} \\
E_{Ne} = & - \frac{1}{{4\pi \varepsilon _0 }}\sum\limits_{i,j} {\frac{{Z_i }}{{R_{ij} }}Erf\left( {\frac{{\sqrt 2 R_{ij} }}{{s_j }}} \right)} \\
E_{ee} = & \frac{1}{{4\pi \varepsilon _0 }}\sum\limits_{i < j} {\frac{1}{{r_{ij} }}Erf\left( {\frac{{\sqrt 2 r_{ij} }}{{\sqrt {s_i^2 + s_j^2 } }}} \right)} \\
E_{Pauli} = & \sum\limits_{\sigma _i = \sigma _j } {E\left( { \uparrow \uparrow } \right)_{ij}} + \sum\limits_{\sigma _i \ne \sigma _j } {E\left( { \uparrow \downarrow } \right)_{ij}} \\
where, s\_i correspond to the electron sizes, the sigmas i's to the
fixed spins of the electrons, Z\_i to the charges on the nuclei, R\_ij
@ -229,11 +223,9 @@ representations, after the "ecp" keyword.
Si. The ECP captures the orbital overlap between the core and valence
electrons (i.e. Pauli repulsion) with one of the functional forms:
.. image:: Eqs/eff_ECP1.jpg
:align: center
.. image:: Eqs/eff_ECP2.jpg
:align: center
.. math::
E_{Pauli(ECP_s)} = & p_1\exp\left(-\frac{p_2r^2}{p_3+s^2} \right) \\
E_{Pauli(ECP_p)} = & p_1\left( \frac{2}{p_2/s+s/p_2} \right)\left( r-p_3s\right)^2\exp \left[ -\frac{p_4\left( r-p_3s \right)^2}{p_5+s^2} \right]
Where the 1st form correspond to core interactions with s-type valence
electrons and the 2nd to core interactions with p-type valence

View File

@ -34,30 +34,42 @@ Style *eim* computes pairwise interactions for ionic compounds
using embedded-ion method (EIM) potentials :ref:`(Zhou) <Zhou2>`. The
energy of the system E is given by
.. image:: Eqs/pair_eim1.jpg
:align: center
.. math::
E = \frac{1}{2} \sum_{i=1}^{N} \sum_{j=i_1}^{i_N} \phi_{ij} \left(r_{ij}\right) + \sum_{i=1}^{N}E_i\left(q_i,\sigma_i\right)
The first term is a double pairwise sum over the J neighbors of all I
atoms, where phi\_ij is a pair potential. The second term sums over
atoms, where :math:`\phi_{ij}` is a pair potential. The second term sums over
the embedding energy E\_i of atom I, which is a function of its charge
q\_i and the electrical potential sigma\_i at its location. E\_i, q\_i,
and sigma\_i are calculated as
q\_i and the electrical potential :math:`\sigma_i` at its location. E\_i, q\_i,
and :math:`sigma_i` are calculated as
.. image:: Eqs/pair_eim2.jpg
:align: center
.. math::
where eta\_ji is a pairwise function describing electron flow from atom
I to atom J, and psi\_ij is another pairwise function. The multi-body
q_i = & \sum_{j=i_1}^{i_N} \eta_{ji}\left(r_{ij}\right) \\
\sigma_i = & \sum_{j=i_1}^{i_N} q_j \cdot \psi_{ij} \left(r_{ij}\right) \\
E_i\left(q_i,\sigma_i\right) = & \frac{1}{2} \cdot q_i \cdot \sigma_i
where :math:`\eta_{ji} is a pairwise function describing electron flow from atom
I to atom J, and :math:`\psi_{ij}` is another pairwise function. The multi-body
nature of the EIM potential is a result of the embedding energy term.
A complete list of all the pair functions used in EIM is summarized
below
.. image:: Eqs/pair_eim3.jpg
:align: center
.. math::
Here E\_b, r\_e, r\_(c,phi), alpha, beta, A\_(psi), zeta, r\_(s,psi),
r\_(c,psi), A\_(eta), r\_(s,eta), r\_(c,eta), chi, and pair function type
p are parameters, with subscripts ij indicating the two species of
\phi_{ij}\left(r\right) = & \left\{ \begin{array}{lr}
\left[\frac{E_{b,ij}\beta_{ij}}{\beta_{ij}-\alpha_{ij}}\exp\left(-\alpha_{ij} \frac{r-r_{e,ij}}{r_{e,ij}}\right)-\frac{E_{b,ij}\alpha_{ij}}{\beta_{ij}-\alpha_{ij}}\exp\left(-\beta_{ij} \frac{r-r_{e,ij}}{r_{e,ij}}\right)\right]f_c\left(r,r_{e,ij},r_{c,\phi,ij}\right),& p_{ij}=1 \\
\left[\frac{E_{b,ij}\beta_{ij}}{\beta_{ij}-\alpha_{ij}} \left(\frac{r_{e,ij}}{r}\right)^{\alpha_{ij}} -\frac{E_{b,ij}\alpha_{ij}}{\beta_{ij}-\alpha_{ij}} \left(\frac{r_{e,ij}}{r}\right)^{\beta_{ij}}\right]f_c\left(r,r_{e,ij},r_{c,\phi,ij}\right),& p_{ij}=2
\end{array}
\right.\\
\eta_{ji} = & A_{\eta,ij}\left(\chi_j-\chi_i\right)f_c\left(r,r_{s,\eta,ij},r_{c,\eta,ij}\right) \\
\psi_{ij}\left(r\right) = & A_{\psi,ij}\exp\left(-\zeta_{ij}r\right)f_c\left(r,r_{s,\psi,ij},r_{c,\psi,ij}\right) \\
f_{c}\left(r,r_p,r_c\right) = & 0.510204 \mathrm{erfc}\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204
Here :math:`E_b, r_e, r_(c,\phi), \alpha, \beta, A_(\psi), \zeta, r_(s,\psi),
r_(c,\psi), A_(\eta), r_(s,\eta), r_(c,\eta), \chi,` and pair function type
*p* are parameters, with subscripts *ij* indicating the two species of
atoms in the atomic pair.
.. note::

View File

@ -43,17 +43,20 @@ one CG particle can interact with a species in a neighboring CG
particle through a site-site interaction potential model. The
*exp6/rx* style computes an exponential-6 potential given by
.. image:: Eqs/pair_exp6_rx.jpg
:align: center
.. math::
where the *epsilon* parameter determines the depth of the potential
minimum located at *Rm*\ , and *alpha* determines the softness of the repulsion.
U_{ij}(r) = \frac{\epsilon}{\alpha-6}\{6\exp[\alpha(1-\frac{r_{ij}}{R_{m}})]-\alpha(\frac{R_{m}}{r_{ij}})^6\}
where the :math:`\epsilon` parameter determines the depth of the
potential minimum located at :math:`R_m`, and :math:`\alpha` determines
the softness of the repulsion.
The coefficients must be defined for each species in a given particle
type via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above, where the first argument is the filename that includes the
exponential-6 parameters for each species. The file includes the
species tag followed by the *alpha*\ , *epsilon* and *Rm*
species tag followed by the :math:`\alpha, \epsilon` and :math:`R_m`
parameters. The format of the file is described below.
The second and third arguments specify the site-site interaction
@ -74,22 +77,22 @@ to scale the EXP-6 parameters as reactions occur. Currently, there
are three scaling options: *exponent*\ , *polynomial* and *none*\ .
Exponent scaling requires two additional arguments for scaling
the *Rm* and *epsilon* parameters, respectively. The scaling factor
the :math:`R_m` and :math:`\epsilon` parameters, respectively. The scaling factor
is computed by phi\^exponent, where phi is the number of molecules
represented by the coarse-grain particle and exponent is specified
as a pair coefficient argument for *Rm* and *epsilon*\ , respectively.
The *Rm* and *epsilon* parameters are multiplied by the scaling
as a pair coefficient argument for :math:`R_m` and :math:`\epsilon`, respectively.
The :math:`R_m` and :math:`\epsilon` parameters are multiplied by the scaling
factor to give the scaled interaction parameters for the CG particle.
Polynomial scaling requires a filename to be specified as a pair
coeff argument. The file contains the coefficients to a fifth order
polynomial for the *alpha*\ , *epsilon* and *Rm* parameters that depend
polynomial for the :math:`\alpha`, :math:`\epsilon` and :math:`R_m` parameters that depend
upon phi (the number of molecules represented by the CG particle).
The format of a polynomial file is provided below.
The *none* option to the scaling does not have any additional pair coeff
arguments. This is equivalent to specifying the *exponent* option with
*Rm* and *epsilon* exponents of 0.0 and 0.0, respectively.
:math:`R_m` and :math:`\epsilon` exponents of 0.0 and 0.0, respectively.
The final argument specifies the interaction cutoff (optional).
@ -133,23 +136,30 @@ between sections.
Following a blank line, the next N lines list the species and their
corresponding parameters. The first argument is the species tag, the
second argument is the exp6 tag, the 3rd argument is the *alpha*
parameter (energy units), the 4th argument is the *epsilon* parameter
(energy-distance\^6 units), and the 5th argument is the *Rm* parameter
second argument is the exp6 tag, the 3rd argument is the :math:`\alpha`
parameter (energy units), the 4th argument is the :math:`\epsilon` parameter
(energy-distance\^6 units), and the 5th argument is the :math:`R_m` parameter
(distance units). If a species tag of "1fluid" is listed as a pair
coefficient, a one-fluid approximation is specified where a
concentration-dependent combination of the parameters is computed
through the following equations:
.. image:: Eqs/pair_exp6_rx_oneFluid.jpg
:align: center
.. math::
R_{m}^{3} = & \sum_{a}\sum_{b} x_{a}x_{b}R_{m,ab}^{3} \\
\epsilon = & \frac{1}{R_{m}^{3}}\sum_{a}\sum_{b} x_{a}x_{b}\epsilon_{ab}R_{m,ab}^{3} \\
\alpha = & \frac{1}{\epsilon R_{m}^{3}}\sum_{a}\sum_{b} x_{a}x_{b}\alpha_{ab}\epsilon_{ab}R_{m,ab}^{3}
where
.. image:: Eqs/pair_exp6_rx_oneFluid2.jpg
:align: center
.. math::
and xa and xb are the mole fractions of a and b, respectively, which
\epsilon_{ab} = & \sqrt{\epsilon_{a}\epsilon_{b}} \\
R_{m,ab} = & \frac{R_{m,a}+R_{m,b}}{2} \\
\alpha_{ab} = & \sqrt{\alpha_{a}\alpha_{b}}
and :math:`x_a` and :math:`x_b` are the mole fractions of a and b, respectively, which
comprise the gas mixture.