lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2000 from akohlmey/polymorphic-update 

Update for pair style polymorphic from Xiaowang Zhou
This commit is contained in:
Axel Kohlmeyer 2020-04-20 17:14:47 -04:00 committed by GitHub
parent 0766350838 dcc074215d
commit 0711232e5b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 75608 additions and 242 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
doc/src/pair_momb.rst
View File

@ -31,7 +31,7 @@ Style *momb* computes pairwise van der Waals (vdW) and short-range
interactions using the Morse potential and :ref:`(Grimme) <Grimme>` method interactions using the Morse potential and :ref:`(Grimme) <Grimme>` method
implemented in the Many-Body Metal-Organic (MOMB) force field implemented in the Many-Body Metal-Organic (MOMB) force field
described comprehensively in :ref:`(Fichthorn) <Fichthorn>` and described comprehensively in :ref:`(Fichthorn) <Fichthorn>` and
:ref:`(Zhou) <Zhou4>`. Grimme's method is widely used to correct for :ref:`(Zhou) <Zhou5>`. Grimme's method is widely used to correct for
dispersion in density functional theory calculations. dispersion in density functional theory calculations.
.. math:: .. math::
@ -76,6 +76,6 @@ Related commands
**(Fichthorn)** Fichthorn, Balankura, Qi, CrystEngComm, 18(29), 5410-5417 (2016). **(Fichthorn)** Fichthorn, Balankura, Qi, CrystEngComm, 18(29), 5410-5417 (2016).
.. _Zhou4: .. _Zhou5:
**(Zhou)** Zhou, Saidi, Fichthorn, J Phys Chem C, 118(6), 3366-3374 (2014). **(Zhou)** Zhou, Saidi, Fichthorn, J Phys Chem C, 118(6), 3366-3374 (2014).

376
doc/src/pair_polymorphic.rst
View File

@ -18,21 +18,22 @@ Examples
.. code-block:: LAMMPS .. code-block:: LAMMPS
pair_style polymorphic pair_style polymorphic
pair_coeff * * TlBr_msw.polymorphic Tl Br pair_coeff * * FeCH_BOP_I.poly Fe C H
pair_coeff * * AlCu_eam.polymorphic Al Cu pair_coeff * * TlBr_msw.poly Tl Br
pair_coeff * * GaN_tersoff.polymorphic Ga N pair_coeff * * CuTa_eam.poly Cu Ta
pair_coeff * * GaN_sw.polymorphic GaN pair_coeff * * GaN_tersoff.poly Ga N
pair_coeff * * GaN_sw.poly Ga N
Description Description
""""""""""" """""""""""
The *polymorphic* pair style computes a 3-body free-form potential The *polymorphic* pair style computes a 3-body free-form potential
(:ref:`Zhou <Zhou3>`) for the energy E of a system of atoms as (:ref:`Zhou3 <Zhou3>`) for the energy E of a system of atoms as
.. math:: .. math::
E & = \frac{1}{2}\sum_{i=1}^{i=N}\sum_{j=1}^{j=N}\left[\left(1-\delta_{ij}\right)\cdot U_{IJ}\left(r_{ij}\right)-\left(1-\eta_{ij}\right)\cdot F_{IJ}\left(r_{ij}\right)\cdot V_{IJ}\left(r_{ij}\right)\right] \\ E & = \frac{1}{2}\sum_{i=1}^{i=N}\sum_{j=1}^{j=N}\left[\left(1-\delta_{ij}\right)\cdot U_{IJ}\left(r_{ij}\right)-\left(1-\eta_{ij}\right)\cdot F_{IJ}\left(X_{ij}\right)\cdot V_{IJ}\left(r_{ij}\right)\right] \\
X_{ij} & = \sum_{k=i_1,k\neq i,j}^{i_N}W_{IK}\left(r_{ik}\right)\cdot G_{JIK}\left(\theta_{jik}\right)\cdot P_{IK}\left(\Delta r_{jik}\right) \\ X_{ij} & = \sum_{k=i_1,k\neq j}^{i_N}W_{IK}\left(r_{ik}\right)\cdot G_{JIK}\left(\cos\theta_{jik}\right)\cdot P_{JIK}\left(\Delta r_{jik}\right) \\
\Delta r_{jik} & = r_{ij}-\xi_{IJ}\cdot r_{ik} \Delta r_{jik} & = r_{ij}-\xi_{IJ}\cdot r_{ik}
where I, J, K represent species of atoms i, j, and k, :math:`i_1, ..., where I, J, K represent species of atoms i, j, and k, :math:`i_1, ...,
@ -42,230 +43,267 @@ Dirac constant (i.e., :math:`\delta_{ij} = 1` when :math:`i = j`, and
constant that can be set either to :math:`\eta_{ij} = \delta_{ij}` or constant that can be set either to :math:`\eta_{ij} = \delta_{ij}` or
:math:`\eta_{ij} = 1 - \delta_{ij}` depending on the potential type, :math:`\eta_{ij} = 1 - \delta_{ij}` depending on the potential type,
:math:`U_{IJ}(r_{ij})`, :math:`V_{IJ}(r_{ij})`, :math:`W_{IK}(r_{ik})` :math:`U_{IJ}(r_{ij})`, :math:`V_{IJ}(r_{ij})`, :math:`W_{IK}(r_{ik})`
are pair functions, :math:`G_{JIK}(\cos(\theta))` is an angular are pair functions, :math:`G_{JIK}(\cos\theta_{jik})` is an angular
function, :math:`P_{IK}(\Delta r_{jik})` is a function of atomic spacing function, :math:`P_{JIK}(\Delta r_{jik})` is a function of atomic
differential :math:`\Delta r_{jik} = r_{ij} - \xi_{IJ} \cdot r_{ik}` spacing differential :math:`\Delta r_{jik} = r_{ij} - \xi_{IJ} \cdot
with :math:`\xi_{IJ}` being a pair-dependent parameter, and r_{ik}` with :math:`\xi_{IJ}` being a pair-dependent parameter, and
:math:`F_{IJ}(X_{ij})` is a function of the local environment variable :math:`F_{IJ}(X_{ij})` is a function of the local environment variable
:math:`X_{ij}`. This generic potential is fully defined once the :math:`X_{ij}`. This generic potential is fully defined once the
constants :math:`\eta_{ij}` and :math:`\xi_{IJ}`, and the six functions constants :math:`\eta_{ij}` and :math:`\xi_{IJ}`, and the six functions
:math:`U_{IJ}(r_{ij})`, :math:`V_{IJ}(r_{ij})`, :math:`W_{IK}(r_{ik})`, :math:`U_{IJ}(r_{ij})`, :math:`V_{IJ}(r_{ij})`, :math:`W_{IK}(r_{ik})`,
:math:`G_{JIK}(\cos(\theta))`, :math:`P_{IK}(\Delta r_{jik})`, and :math:`G_{JIK}(\cos\theta_{jik})`, :math:`P_{JIK}(\Delta r_{jik})`, and
:math:`F_{IJ}(X_{ij})` are given. Note that these six functions are all :math:`F_{IJ}(X_{ij})` are given. Here LAMMPS uses a global parameter
one dimensional, and hence can be provided in an analytic or tabular :math:`\eta` to represent :math:`\eta_{ij}`. When :math:`\eta = 1`,
form. This allows users to design different potentials solely based on a :math:`\eta_{ij} = 1 - \delta_{ij}`, otherwise :math:`\eta_{ij} =
manipulation of these functions. For instance, the potential reduces to \delta_{ij}`. Additionally, :math:`\eta = 3` indicates that the function
Stillinger-Weber potential (:ref:`SW <SW>`) if we set :math:`P_{JIK}(\Delta r)` depends on species I, J and K, otherwise
:math:`P_{JIK}(\Delta r) = P_{IK}(\Delta r)` only depends on species I
and K. Note that these six functions are all one dimensional, and hence
can be provided in a tabular form. This allows users to design different
potentials solely based on a manipulation of these functions. For
instance, the potential reduces to a Stillinger-Weber potential
(:ref:`SW <SW>`) if we set
.. math:: .. math::
\left\{\begin{array}{l} \eta_{ij} & = \delta_{ij} (\eta = 2~or~\eta = 0),\xi_{IJ}=0 \\
\eta_{ij} = \delta_{ij},\xi_{IJ}=0 \\ U_{IJ}\left(r\right) & = A_{IJ}\cdot\epsilon_{IJ}\cdot \left(\frac{\sigma_{IJ}}{r}\right)^q\cdot \left[B_{IJ}\cdot \left(\frac{\sigma_{IJ}}{r}\right)^{p-q}-1\right]\cdot exp\left(\frac{\sigma_{IJ}}{r-a_{IJ}\cdot \sigma_{IJ}}\right) \\
U_{IJ}\left(r\right)=A_{IJ}\cdot\epsilon_{IJ}\cdot \left(\frac{\sigma_{IJ}}{r}\right)^q\cdot \left[B_{IJ}\cdot \left(\frac{\sigma_{IJ}}{r}\right)^{p-q}-1\right]\cdot exp\left(\frac{\sigma_{IJ}}{r-a_{IJ}\cdot \sigma_{IJ}}\right) \\ V_{IJ}\left(r\right) & = \sqrt{\lambda_{IJ}\cdot \epsilon_{IJ}}\cdot exp\left(\frac{\gamma_{IJ}\cdot \sigma_{IJ}}{r-a_{IJ}\cdot \sigma_{IJ}}\right) \\
V_{IJ}\left(r\right)=\sqrt{\lambda_{IJ}\cdot \epsilon_{IJ}}\cdot exp\left(\frac{\gamma_{IJ}\cdot \sigma_{IJ}}{r-a_{IJ}\cdot \sigma_{IJ}}\right) \\ F_{IJ}\left(X\right) & = -X \\
F_{IJ}\left(X\right)=-X \\ P_{JIK}\left(\Delta r\right) & = P_{IK}\left(\Delta r\right) = 1 \\
P_{IJ}\left(\Delta r\right)=1 \\ W_{IJ}\left(r\right) & = \sqrt{\lambda_{IJ}\cdot \epsilon_{IJ}}\cdot exp\left(\frac{\gamma_{IJ}\cdot \sigma_{IJ}}{r-a_{IJ}\cdot \sigma_{IJ}}\right) \\
W_{IJ}\left(r\right)=\sqrt{\lambda_{IJ}\cdot \epsilon_{IJ}}\cdot exp\left(\frac{\gamma_{IJ}\cdot \sigma_{IJ}}{r-a_{IJ}\cdot \sigma_{IJ}}\right) \\ G_{JIK}\left(\cos\theta\right) & = \left(\cos\theta+\frac{1}{3}\right)^2
G_{JIK}\left(\theta\right)=\left(cos\theta+\frac{1}{3}\right)^2
\end{array}\right.
The potential reduces to Tersoff types of potential The potential reduces to a Tersoff potential (:ref:`Tersoff <Tersoff>`
(:ref:`Tersoff <Tersoff>` or :ref:`Albe <poly-Albe>`) if we set or :ref:`Albe <poly-Albe>`) if we set
.. math:: .. math::
\left\{\begin{array}{l} \eta_{ij} & = \delta_{ij} (\eta = 2~or~\eta = 0),\xi_{IJ}=1 \\
\eta_{ij}=\delta_{ij},\xi_{IJ}=1 \\ U_{IJ}\left(r\right) & = \frac{D_{e,IJ}}{S_{IJ}-1}\cdot exp\left[-\beta_{IJ}\sqrt{2S_{IJ}}\left(r-r_{e,IJ}\right)\right]\cdot f_{c,IJ}\left(r\right) \\
U_{IJ}\left(r\right)=\frac{D_{e,IJ}}{S_{IJ}-1}\cdot exp\left[-\beta_{IJ}\sqrt{2S_{IJ}\left(r-r_{e,IJ}\right)}\right]\cdot f_{c,IJ}\left(r\right) \\ V_{IJ}\left(r\right) & = \frac{S_{IJ}\cdot D_{e,IJ}}{S_{IJ}-1}\cdot exp\left[-\beta_{IJ}\sqrt{\frac{2}{S_{IJ}}}\left(r-r_{e,IJ}\right)\right]\cdot f_{c,IJ}\left(r\right) \\
V_{IJ}\left(r\right)=\frac{S_{IJ}\cdot D_{e,IJ}}{S_{IJ}-1}\cdot exp\left[-\beta_{IJ}\sqrt{\frac{2}{S_{IJ}}\left(r-r_{e,IJ}\right)}\right]\cdot f_{c,IJ}\left(r\right) \\ F_{IJ}\left(X\right) & = \left(1+X\right)^{-\frac{1}{2}} \\
F_{IJ}\left(X\right)=\left(1+X\right)^{-\frac{1}{2}} \\ P_{JIK}\left(\Delta r\right) & = P_{IK}\left(\Delta r\right) = exp\left(2\mu_{IK}\cdot \Delta r\right) \\
P_{IJ}\left(\Delta r\right)=exp\left(2\mu_{IK}\cdot \Delta r\right) \\ W_{IJ}\left(r\right) & = f_{c,IJ}\left(r\right) \\
W_{IJ}\left(r\right)=f_{c,IK}\left(r\right) \\ G_{JIK}\left(\cos\theta\right) & = \gamma_{IK}\left[1+\frac{c_{IK}^2}{d_{IK}^2}-\frac{c_{IK}^2}{d_{IK}^2+\left(h_{IK}+\cos\theta\right)^2}\right]
G_{JIK}\left(\theta\right)=\gamma_{IK}\left[1+\frac{c_{IK}^2}{d_{IK}^2}-\frac{c_{IK}^2}{d_{IK}^2+\left(h_{IK}+cos\theta\right)^2}\right]
\end{array}\right. where
.. math:: .. math::
f_{c,IJ}=\left\{\begin{array}{lr} f_{c,IJ}\left(r\right)=\left\{\begin{array}{l}
1, & r\leq r_{s,IJ} \\ 1, r\leq R_{IJ}-D_{IJ} \\
\frac{1}{2}+\frac{1}{2} cos \left[\frac{\pi \left(r-r_{s,IJ}\right)}{r_{c,IJ}-r_{s,IJ}}\right], & r_{s,IJ}<r<r_{c,IJ} \\ \frac{1}{2}+\frac{1}{2}cos\left[\frac{\pi\left(r+D_{IJ}-R_{IJ}\right)}{2D_{IJ}}\right], R_{IJ}-D_{IJ} < r < R_{IJ}+D_{IJ} \\
0, & r \geq r_{c,IJ} \\ 0, r \geq R_{IJ}+D_{IJ}
\end{array}\right. \end{array}\right.
The potential reduces to Rockett-Tersoff (:ref:`Wang <Wang3>`) type if we set The potential reduces to a modified Stillinger-Weber potential
(:ref:`Zhou3 <Zhou3>`) if we set
.. math:: .. math::
\left\{\begin{array}{l} \eta_{ij} & = \delta_{ij} (\eta = 2~or~\eta = 0),\xi_{IJ}=0 \\
\eta_{ij}=\delta_{ij},\xi_{IJ}=1 \\ U_{IJ}\left(r\right) & = \varphi_{R,IJ}\left(r\right)-\varphi_{A,IJ}\left(r\right) \\
U_{IJ}\left(r\right)=\left\{\begin{array}{lr} V_{IJ}\left(r\right) & = u_{IJ}\left(r\right) \\
A_{IJ}\cdot exp\left(-\lambda_{1,IJ}\cdot r\right)\cdot f_{c,IJ}\left(r\right), & r\leq r_{s,1,IJ} \\ F_{IJ}\left(X\right) & = -X \\
A_{IJ}\cdot exp\left(-\lambda_{1,IJ}\cdot r\right)\cdot f_{c,IJ}\left(r\right)\cdot f_{c,1,IJ}\left(r\right), & r_{s,1,IJ}<r<r_{c,1,IJ} \\ P_{JIK}\left(\Delta r\right) & = P_{IK}\left(\Delta r\right) = 1 \\
0, & r\ge r_{c,1,IJ} W_{IJ}\left(r\right) & = u_{IJ}\left(r\right) \\
\end{array}\right. \\ G_{JIK}\left(\cos\theta\right) & = g_{JIK}\left(\cos\theta\right)
V_{IJ}\left(r\right)=\left\{\begin{array}{lr}
B_{IJ} \cdot exp\left(-\lambda_{2,IJ}\cdot r\right)\cdot f_{c,IJ}\left(r\right), & r\le r_{s,1,IJ} \\ The potential reduces to a Rockett-Tersoff potential (:ref:`Wang3
B_{IJ} \cdot exp\left(-\lambda_{2,IJ}\cdot r\right)\cdot f_{c,IJ}\left(r\right)+A_{IJ}\cdot exp\left(-\lambda_{1,IJ}\cdot r\right)\cdot & \\ ~~~~~~ f_{c,IJ}\left(r\right)\cdot \left[1-f_{c,1,IJ}\left(r\right)\right], & r_{s,1,IJ}<r<r_{c,1,IJ} \\ <Wang3>`) if we set
B_{IJ} \cdot exp\left(-\lambda_{2,IJ}\cdot r\right)\cdot f_{c,IJ}\left(r\right)+A_{IJ}\cdot exp\left(-\lambda_{1,IJ}\cdot r\right)\cdot & \\ ~~~~~~ f_{c,IJ}\left(r\right) & r \ge r_{c,1,IJ}
\end{array}\right. \\
F_{IJ}\left(X\right)=\left[1+\left(\beta_{IJ}\cdot X\right)^{n_{IJ}}\right]^{-\frac{1}{2n_{IJ}}} \\
P_{IJ}\left(\Delta r\right)=exp\left(\lambda_{3,IK}\cdot \Delta r^3\right) \\
W_{IJ}\left(r\right)=f_{c,IK}\left(r\right) \\
G_{JIK}\left(\theta\right)=1+\frac{c_{IK}^2}{d_{IK}^2}-\frac{c_{IK}^2}{d_{IK}^2+\left(h_{IK}+cos\theta\right)^2}
\end{array}\right.
.. math:: .. math::
f_{c,IJ}=\left\{\begin{array}{lr} \eta_{ij} & = \delta_{ij} (\eta = 2~or~\eta = 0),\xi_{IJ}=1 \\
1, & r\leq r_{s,IJ} \\ U_{IJ}\left(r\right) & = A_{IJ}exp\left(-\lambda_{1,IJ}\cdot r\right)f_{c,IJ}\left(r\right)f_{ca,IJ}\left(r\right) \\
\frac{1}{2}+\frac{1}{2} cos \left[\frac{\pi \left(r-r_{s,IJ}\right)}{r_{c,IJ}-r_{s,IJ}}\right], & r_{s,IJ}<r<r_{c,IJ} \\ V_{IJ}\left(r\right) & = \left\{\begin{array}{l}B_{IJ}exp\left(-\lambda_{2,IJ}\cdot r\right)f_{c,IJ}\left(r\right)+ \\ A_{IJ}exp\left(-\lambda_{1,IJ}\cdot r\right)f_{c,IJ}\left(r\right) \left[1-f_{ca,IJ}\left(r\right)\right]\end{array} \right\} \\
0, & r \geq r_{c,IJ} \\ F_{IJ}\left(X\right) & = \left[1+\left(\beta_{IJ}X\right)^{n_{IJ}}\right]^{-\frac{1}{2n_{IJ}}} \\
\end{array}\right. P_{JIK}\left(\Delta r\right) & = P_{IK}\left(\Delta r\right) = exp\left(\lambda_{3,IK}\cdot \Delta r^3\right) \\
W_{IJ}\left(r\right) & = f_{c,IJ}\left(r\right) \\
G_{JIK}\left(\cos\theta\right) & = 1+\frac{c_{IK}^2}{d_{IK}^2}-\frac{c_{IK}^2}{d_{IK}^2+\left(h_{IK}+\cos\theta\right)^2}
where :math:`f_{ca,IJ}(r)` is similar to the :math:`f_{c,IJ}(r)` defined
above:
.. math:: .. math::
f_{c,1,IJ}=\left\{\begin{array}{lr} f_{ca,IJ}\left(r\right)=\left\{\begin{array}{l}
1, & r\leq r_{s,1,IJ} \\ 1, r\leq R_{a,IJ}-D_{a,IJ} \\
\frac{1}{2}+\frac{1}{2} cos \left[\frac{\pi \left(r-r_{s,1,IJ}\right)}{r_{c,1,IJ}-r_{s,1,IJ}}\right], & r_{s,1,IJ}<r<r_{c,1,IJ} \\ \frac{1}{2}+\frac{1}{2}cos\left[\frac{\pi\left(r+D_{a,IJ}-R_{a,IJ}\right)}{2D_{a,IJ}}\right], R_{a,IJ}-D_{a,IJ} < r < R_{a,IJ}+D_{a,IJ} \\
0, & r \geq r_{c,1,IJ} \\ 0, r \geq R_{a,IJ}+D_{a,IJ}
\end{array}\right. \end{array}\right.
The potential becomes embedded atom method (:ref:`Daw <poly-Daw>`) if we set The potential becomes the embedded atom method (:ref:`Daw <poly-Daw>`)
if we set
.. math:: .. math::
\left\{\begin{array}{l} \eta_{ij} & = 1-\delta_{ij} (\eta = 1),\xi_{IJ}=0 \\
\eta_{ij}=1-\delta_{ij},\xi_{IJ}=0 \\ U_{IJ}\left(r\right) & = \phi_{IJ}\left(r\right) \\
U_{IJ}\left(r\right)=\phi_{IJ}\left(r\right) \\ V_{IJ}\left(r\right) & = 1 \\
V_{IJ}\left(r\right)=1 \\ F_{II}\left(X\right) & = -2F_I\left(X\right) \\
F_{II}\left(X\right)=-2F_I\left(X\right) \\ P_{JIK}\left(\Delta r\right) & = P_{IK}\left(\Delta r\right) = 1 \\
P_{IJ}\left(\Delta r\right)=1 \\ W_{IJ}\left(r\right) & = f_{J}\left(r\right) \\
W_{IJ}\left(r\right)=f_{K}\left(r\right) \\ G_{JIK}\left(\cos\theta\right) & = 1
G_{JIK}\left(\theta\right)=1
\end{array}\right.
In the embedded atom method case, :math:`\phi_{IJ}(r_{ij})` is the pair In the embedded atom method case, :math:`\phi_{IJ}(r)` is the pair
energy, :math:`F_I(X)` is the embedding energy, *X* is the local energy, :math:`F_I(X)` is the embedding energy, *X* is the local
electron density, and :math:`f_K(r)` is the atomic electron density function. electron density, and :math:`f_J(r)` is the atomic electron density
function.
If the tabulated functions are created using the parameters of sw, The potential reduces to another type of Tersoff potential (:ref:`Zhou4
tersoff, and eam potentials, the polymorphic pair style will produce <Zhou4>`) if we set
the same global properties (energies and stresses) and the same forces
as the sw, tersoff, and eam pair styles. The polymorphic pair style
also produces the same atom properties (energies and stresses) as the
corresponding tersoff and eam pair styles. However, due to a different
partition of global properties to atom properties, the polymorphic
pair style will produce different atom properties (energies and
stresses) as the sw pair style. This does not mean that polymorphic
pair style is different from the sw pair style in this case. It just
means that the definitions of the atom energies and atom stresses are
different.
Only a single pair_coeff command is used with the polymorphic style .. math::
which specifies an potential file for all needed elements. These are
\eta_{ij} & = \delta_{ij} (\eta = 3),\xi_{IJ}=1 \\
U_{IJ}\left(r\right) & = \frac{D_{e,IJ}}{S_{IJ}-1}\cdot exp\left[-\beta_{IJ}\sqrt{2S_{IJ}}\left(r-r_{e,IJ}\right)\right]\cdot f_{c,IJ}\left(r\right) \cdot T_{IJ}\left(r\right)+V_{ZBL,IJ}\left(r\right)\left[1-T_{IJ}\left(r\right)\right] \\
V_{IJ}\left(r\right) & = \frac{S_{IJ}\cdot D_{e,IJ}}{S_{IJ}-1}\cdot exp\left[-\beta_{IJ}\sqrt{\frac{2}{S_{IJ}}}\left(r-r_{e,IJ}\right)\right]\cdot f_{c,IJ}\left(r\right) \cdot T_{IJ}\left(r\right) \\
F_{IJ}\left(X\right) & = \left(1+X\right)^{-\frac{1}{2}} \\
P_{JIK}\left(\Delta r\right) & = \omega_{JIK} \cdot exp\left(\alpha_{JIK}\cdot \Delta r\right) \\
W_{IJ}\left(r\right) & = f_{c,IJ}\left(r\right) \\
G_{JIK}\left(\cos\theta\right) & = \gamma_{JIK}\left[1+\frac{c_{JIK}^2}{d_{JIK}^2}-\frac{c_{JIK}^2}{d_{JIK}^2+\left(h_{JIK}+\cos\theta\right)^2}\right] \\
T_{IJ}\left(r\right) & = \frac{1}{1+exp\left[-b_{f,IJ}\left(r-r_{f,IJ}\right)\right]} \\
V_{ZBL,IJ}\left(r\right) & = 14.4 \cdot \frac{Z_I \cdot Z_J}{r}\sum_{k=1}^{4}\mu_k \cdot exp\left[-\nu_k \left(Z_I^{0.23}+Z_J^{0.23}\right) r\right]
where :math:`f_{c,IJ}(r)` is the same as defined above. This Tersoff
potential differs from the one above because the :math:`P_{JIK}(\Delta
r)` function is now dependent on all three species I, J, and K.
If the tabulated functions are created using the parameters of
Stillinger-Weber, Tersoff, and EAM potentials, the polymorphic pair
style will produce the same global properties (energies and stresses)
and the same forces as the :doc:`sw <pair_sw>`, :doc:`tersoff
<pair_tersoff>`, and :doc:`eam <pair_eam>` pair styles. The polymorphic
pair style also produces the same per-atom properties (energies and
stresses) as the corresponding :doc:`tersoff <pair_tersoff>` and
:doc:`eam <pair_eam>` pair styles. However, due to a different
partitioning of global properties to per-atom properties, the
polymorphic pair style will produce different per-atom properties
(energies and stresses) as the :doc:`sw <pair_sw>` pair style. This does
not mean that polymorphic pair style is different from the sw pair
style. It just means that the definitions of the atom energies and atom
stresses are different.
Only a single pair_coeff command is used with the polymorphic pair style
which specifies a potential file for all needed elements. These are
mapped to LAMMPS atom types by specifying N additional arguments after mapped to LAMMPS atom types by specifying N additional arguments after
the filename in the pair_coeff command, where N is the number of the filename in the pair_coeff command, where N is the number of LAMMPS
LAMMPS atom types: atom types:
* filename * filename
* N element names = mapping of Tersoff elements to atom types * N element names = mapping of polymorphic potential elements to atom types
See the pair_coeff doc page for alternate ways to specify the path for See the pair_coeff doc page for alternate ways to specify the path for
the potential file. Several files for polymorphic potentials are the potential file. Several files for polymorphic potentials are
included in the potentials directory of the LAMMPS distribution. They included in the potentials directory of the LAMMPS distribution. They
have a "poly" suffix. have a "poly" suffix.
As an example, imagine the SiC_tersoff.poly file has tabulated As an example, imagine the GaN_tersoff.poly file has tabulated functions
functions for Si-C tersoff potential. If your LAMMPS simulation has 4 for Ga-N tersoff potential. If your LAMMPS simulation has 4 atoms types
atoms types and you want the 1st 3 to be Si, and the 4th to be C, you and you want the 1st 3 to be Ga, and the 4th to be N, you would use the
would use the following pair_coeff command: following pair_coeff command:
.. code-block:: LAMMPS .. code-block:: LAMMPS
pair_coeff * * SiC_tersoff.poly Si Si Si C pair_coeff * * GaN_tersoff.poly Ga Ga Ga N
The 1st 2 arguments must be \* \* so as to span all LAMMPS atom The first two arguments must be \* \* to span all pairs of LAMMPS atom
types. The first three Si arguments map LAMMPS atom types 1,2,3 to the types. The first three Ga arguments map LAMMPS atom types 1,2,3 to the
Si element in the polymorphic file. The final C argument maps LAMMPS Ga element in the polymorphic file. The final N argument maps LAMMPS
atom type 4 to the C element in the polymorphic file. If a mapping atom type 4 to the N element in the polymorphic file. If a mapping value
value is specified as NULL, the mapping is not performed. This can be is specified as NULL, the mapping is not performed. This can be used
used when an polymorphic potential is used as part of the hybrid pair when an polymorphic potential is used as part of the hybrid pair
style. The NULL values are placeholders for atom types that will be style. The NULL values are placeholders for atom types that will be used
used with other potentials. with other potentials.
Potential files in the potentials directory of the LAMMPS distribution Potential files in the potentials directory of the LAMMPS distribution
have a ".poly" suffix. At the beginning of the files, an unlimited have a ".poly" suffix. At the beginning of the files, an unlimited
number of lines starting with '#' are used to describe the potential number of lines starting with '#' are used to describe the potential and
and are ignored by LAMMPS. The next line lists two numbers: are ignored by LAMMPS. The next line lists two numbers:
.. parsed-literal:: .. parsed-literal::
ntypes :math:`\eta` ntypes eta
Here ntypes represent total number of species defined in the potential Here *ntypes* represent total number of species defined in the potential
file, and :math:`\eta = 0` or 1. The number ntypes must equal the total file, :math:`\eta = 1` reduces to embedded atom method, :math:`\eta = 3`
number of different species defined in the pair_coeff command. When assumes a three species dependent :math:`P_{JIK}(\Delta r)` function,
:math:`\eta = 1`, :math:\eta_{ij}` defined in the potential functions and all other values of :math:`\eta` assume a two species dependent
above is set to :math:`1 - \delta_{ij}`, otherwise :math:`\eta_{ij}` is :math:`P_{JK}(\Delta r)` function. The value of *ntypes* must equal the
set to :math:`\delta_{ij}`. The next ntypes lines each lists two numbers total number of different species defined in the pair_coeff command. The
and a character string representing atomic number, atomic mass, and name next *ntypes* lines each lists two numbers and a character string
of the species of the ntypes elements: representing atomic number, atomic mass, and name of the species of the
ntypes elements:
.. parsed-literal:: .. parsed-literal::
atomic_number atomic-mass element (1) atomic-number atomic-mass element-name(1)
atomic_number atomic-mass element (2) atomic-number atomic-mass element-name(2)
... ...
atomic_number atomic-mass element (ntypes) atomic-number atomic-mass element-name(ntypes)
The next line contains four numbers:
.. parsed-literal::
nr ntheta nx xmax
Here nr is total number of tabular points for radial functions U, V, W,
P, ntheta is total number of tabular points for the angular function G,
nx is total number of tabular points for the function F, xmax is a
maximum value of the argument of function F. Note that the pair
functions :math:`U_{IJ}(r)`, :math:`V_{IJ}(r)`, :math:`W_{IJ}(r)` are
uniformly tabulated between 0 and cutoff distance of the IJ pair,
:math:`G_{JIK}(\cos\theta)` is uniformly tabulated between -1 and 1,
:math:`P_{JIK}(\Delta r)` is uniformly tabulated between -rcmax and
rcmax where rcmax is the maximum cutoff distance of all pairs, and
:math:`F_{IJ}(X)` is uniformly tabulated between 0 and xmax. Linear
extrapolation is assumed if actual simulations exceed these ranges.
The next ntypes\*(ntypes+1)/2 lines contain two numbers: The next ntypes\*(ntypes+1)/2 lines contain two numbers:
.. parsed-literal:: .. parsed-literal::
cut :math:`xi` (1) cut xi(1)
cut :math:`xi` (2) cut xi(2)
... ...
cut :math:`xi` (ntypes\*(ntypes+1)/2) cut xi(ntypes\*(ntypes+1)/2)
Here cut means the cutoff distance of the pair functions, :math:`\xi` is Here cut means the cutoff distance of the pair functions, "xi" is
the same as defined in the potential functions above. The :math:`\xi` as defined in the potential functions above. The
ntypes\*(ntypes+1)/2 lines are related to the pairs according to the ntypes\*(ntypes+1)/2 lines are related to the pairs according to the
sequence of first ii (self) pairs, i = 1, 2, ..., ntypes, and then then sequence of first ii (self) pairs, i = 1, 2, ..., ntypes, and then ij
ij (cross) pairs, i = 1, 2, ..., ntypes-1, and j = i+1, i+2, ..., ntypes (cross) pairs, i = 1, 2, ..., ntypes-1, and j = i+1, i+2, ..., ntypes
(i.e., the sequence of the ij pairs follows 11, 22, ..., 12, 13, 14, (i.e., the sequence of the ij pairs follows 11, 22, ..., 12, 13, 14,
..., 23, 24, ...). ..., 23, 24, ...).
The final blocks of the potential file are the U, V, W, P, G, and F In the final blocks of the potential file, U, V, W, P, G, and F
functions are listed sequentially. First, U functions are given for functions are listed sequentially. First, U functions are given for each
each of the ntypes\*(ntypes+1)/2 pairs according to the sequence of the ntypes\*(ntypes+1)/2 pairs according to the sequence described
described above. For each of the pairs, nr values are listed. Next, above. For each of the pairs, nr values are listed. Next, similar arrays
similar arrays are given for V, W, and P functions. Then G functions are given for V and W functions. If P functions depend only on pair
are given for all the ntypes\*ntypes\*ntypes ijk triplets in a natural species, i.e., :math:`\eta \neq 3`, then P functions are also listed the
sequence i from 1 to ntypes, j from 1 to ntypes, and k from 1 to same way the next. If P functions depend on three species, i.e.,
ntypes (i.e., ijk = 111, 112, 113, ..., 121, 122, 123 ..., 211, 212, :math:`\eta = 3`, then P functions are listed for all the
...). Each of the ijk functions contains ng values. Finally, the F ntypes*ntypes*ntypes IJK triplets in a natural sequence I from 1 to
functions are listed for all ntypes\*(ntypes+1)/2 pairs, each ntypes, J from 1 to ntypes, and K from 1 to ntypes (i.e., IJK = 111,
containing nx values. Either analytic or tabulated functions can be 112, 113, ..., 121, 122, 123 ..., 211, 212, ...). Next, G functions are
specified. Currently, constant, exponential, sine and cosine analytic listed for all the ntypes*ntypes*ntypes IJK triplets similarly. For each
functions are available which are specified with: constant c1 , where of the G functions, ntheta values are listed. Finally, F functions are
f(x) = c1 exponential c1 c2 , where f(x) = c1 exp(c2\*x) sine c1 c2 , listed for all the ntypes*(ntypes+1)/2 pairs in the same sequence as
where f(x) = c1 sin(c2\*x) cos c1 c2 , where f(x) = c1 cos(c2\*x) described above. For each of the F functions, nx values are listed.
Tabulated functions are specified by spline n x1 x2, where n=number of
point, (x1,x2)=range and then followed by n values evaluated uniformly
over these argument ranges. The valid argument ranges of the
functions are between 0 <= r <= cut for the U(r), V(r), W(r)
functions, -cutmax <= delta_r <= cutmax for the P(delta_r) functions,
-1 <= :math:`\cos\theta` <= 1 for the G(:math:`\cos\theta`) functions,
and 0 <= X <= maxX for the F(X) functions.
**Mixing, shift, table tail correction, restart**\ : **Mixing, shift, table tail correction, restart**\ :
This pair styles does not support the :doc:`pair_modify <pair_modify>` This pair styles does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options. shift, table, and tail options.
This pair style does not write their information to :doc:`binary restart files <restart>`, since it is stored in potential files. Thus, you This pair style does not write their information to :doc:`binary restart
need to re-specify the pair_style and pair_coeff commands in an input files <restart>`, since it is stored in potential files. Thus, you need
script that reads a restart file. to re-specify the pair_style and pair_coeff commands in an input script
that reads a restart file.
---------- ----------
@ -277,15 +315,15 @@ input script. If using read_data, atomic masses must be defined in the
atomic structure data file. atomic structure data file.
This pair style is part of the MANYBODY package. It is only enabled if This pair style is part of the MANYBODY package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info. LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` doc page for more info.
This pair potential requires the :doc:`newtion <newton>` setting to be This pair potential requires the :doc:`newtion <newton>` setting to be
"on" for pair interactions. "on" for pair interactions.
The potential files provided with LAMMPS (see the potentials The potential files provided with LAMMPS (see the potentials directory)
directory) are parameterized for metal :doc:`units <units>`. You can use are parameterized for metal :doc:`units <units>`. You can use any LAMMPS
any LAMMPS units, but you would need to create your own potential units, but you would need to create your own potential files.
files.
Related commands Related commands
"""""""""""""""" """"""""""""""""
@ -296,12 +334,15 @@ Related commands
.. _Zhou3: .. _Zhou3:
**(Zhou)** X. W. Zhou, M. E. Foster, R. E. Jones, P. Yang, H. Fan, and **(Zhou3)** X. W. Zhou, M. E. Foster, R. E. Jones, P. Yang, H. Fan, and F. P. Doty, J. Mater. Sci. Res., 4, 15 (2015).
F. P. Doty, J. Mater. Sci. Res., 4, 15 (2015).
.. _Zhou4:
**(Zhou4)** X. W. Zhou, M. E. Foster, J. A. Ronevich, and C. W. San Marchi, J. Comp. Chem., 41, 1299 (2020).
.. _SW: .. _SW:
**(SW)** F. H. Stillinger-Weber, and T. A. Weber, Phys. Rev. B, 31, 5262 (1985). **(SW)** F. H. Stillinger, and T. A. Weber, Phys. Rev. B, 31, 5262 (1985).
.. _Tersoff: .. _Tersoff:
@ -309,8 +350,7 @@ F. P. Doty, J. Mater. Sci. Res., 4, 15 (2015).
.. _poly-Albe: .. _poly-Albe:
**(Albe)** K. Albe, K. Nordlund, J. Nord, and A. Kuronen, Phys. Rev. B, **(Albe)** K. Albe, K. Nordlund, J. Nord, and A. Kuronen, Phys. Rev. B, 66, 035205 (2002).
66, 035205 (2002).
.. _Wang3: .. _Wang3:

3
doc/utils/sphinx-config/false_positives.txt
View File

@ -1663,6 +1663,7 @@ manybody
MANYBODY MANYBODY
Maras Maras
Marchetti Marchetti
Marchi
Mariella Mariella
Marrink Marrink
Marroquin Marroquin
@ -2442,6 +2443,7 @@ Ravelo
rc rc
Rc Rc
Rcm Rcm
rcmax
Rcmx Rcmx
Rcmy Rcmy
Rcut Rcut
@ -2560,6 +2562,7 @@ Rockett
Rodrigues Rodrigues
Rohart Rohart
Ronchetti Ronchetti
Ronevich
Rosati Rosati
Rosenberger Rosenberger
Rossky Rossky

33617
potentials/FeCH_BOP_I.poly Normal file
View File

File diff suppressed because it is too large Load Diff

33618
potentials/FeCH_BOP_II.poly Normal file
View File

File diff suppressed because it is too large Load Diff

8014
potentials/TlBr_msw.poly Normal file
View File

File diff suppressed because it is too large Load Diff

127
src/MANYBODY/pair_polymorphic.cpp
View File

@ -46,6 +46,7 @@ PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp)
nelements = 0; nelements = 0;
elements = NULL; elements = NULL;
match = NULL;
pairParameters = NULL; pairParameters = NULL;
tripletParameters = NULL; tripletParameters = NULL;
elem2param = NULL; elem2param = NULL;
@ -168,7 +169,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
firstneighW1 = new int[neighsize]; firstneighW1 = new int[neighsize];
} }
if (eta) { if (eta == 1) {
iparam_ii = elem2param[itype][itype]; iparam_ii = elem2param[itype][itype];
PairParameters & p = pairParameters[iparam_ii]; PairParameters & p = pairParameters[iparam_ii];
emb = (p.F)->get_vmax(); emb = (p.F)->get_vmax();
@ -193,7 +194,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
PairParameters & p = pairParameters[iparam_ij]; PairParameters & p = pairParameters[iparam_ij];
// do not include the neighbor if get_vmax() <= epsilon because the function is near zero // do not include the neighbor if get_vmax() <= epsilon because the function is near zero
if (eta) { if (eta == 1) {
if (emb > epsilon) { if (emb > epsilon) {
iparam_jj = elem2param[jtype][jtype]; iparam_jj = elem2param[jtype][jtype];
PairParameters & q = pairParameters[iparam_jj]; PairParameters & q = pairParameters[iparam_jj];
@ -255,7 +256,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
evdwl,0.0,fpair,delx,dely,delz); evdwl,0.0,fpair,delx,dely,delz);
} }
if (eta) { if (eta == 1) {
if (emb > epsilon) { if (emb > epsilon) {
@ -356,7 +357,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
PairParameters & q = pairParameters[iparam_ik]; PairParameters & q = pairParameters[iparam_ik];
(q.W)->value(r2,wfac,1,fpair,0); (q.W)->value(r2,wfac,1,fpair,0);
(q.P)->value(r1-(p.xi)*r2,pfac,1,fpair,0); (trip.P)->value(r1-(p.xi)*r2,pfac,1,fpair,0);
(trip.G)->value(costheta,gfac,1,fpair,0); (trip.G)->value(costheta,gfac,1,fpair,0);
zeta_ij += wfac*pfac*gfac; zeta_ij += wfac*pfac*gfac;
@ -397,7 +398,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
iparam_ik = elem2param[itype][ktype]; iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik]; PairParameters & q = pairParameters[iparam_ik];
attractive(&q,&trip,prefactor,r1,r2,delr1,delr2,fi,fj,fk); attractive(&p,&q,&trip,prefactor,r1,r2,delr1,delr2,fi,fj,fk);
f[i][0] += fi[0]; f[i][0] += fi[0];
f[i][1] += fi[1]; f[i][1] += fi[1];
@ -586,7 +587,7 @@ void PairPolymorphic::read_file(char *file)
error->all(FLERR,"Incorrect number of elements in potential file"); error->all(FLERR,"Incorrect number of elements in potential file");
match = new int[nelements]; match = new int[nelements];
ptr = strtok(NULL," \t\n\r\f"); // 1st line, 2nd token ptr = strtok(NULL," \t\n\r\f"); // 1st line, 2nd token
eta = (atoi(ptr)>0) ? true:false; eta = atoi(ptr);
// map the elements in the potential file to LAMMPS atom types // map the elements in the potential file to LAMMPS atom types
for (int i = 0; i < nelements; i++) { for (int i = 0; i < nelements; i++) {
@ -654,17 +655,9 @@ void PairPolymorphic::read_file(char *file)
p.cut = atof(ptr); p.cut = atof(ptr);
p.cutsq = p.cut*p.cut; p.cutsq = p.cut*p.cut;
ptr = strtok(NULL," \t\n\r\f"); // 2nd token ptr = strtok(NULL," \t\n\r\f"); // 2nd token
p.xi = (atoi(ptr)>0) ? true:false; p.xi = atof(ptr);
} }
// set cutmax to max of all params
cutmax = 0.0;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
if (p.cut > cutmax) cutmax = p.cut;
}
cutmaxsq = cutmax*cutmax;
// start reading tabular functions // start reading tabular functions
double * singletable = new double[nr]; double * singletable = new double[nr];
for (int i = 0; i < npair; i++) { // U for (int i = 0; i < npair; i++) { // U
@ -694,14 +687,53 @@ void PairPolymorphic::read_file(char *file)
p.W = new tabularFunction(nr,0.0,p.cut); p.W = new tabularFunction(nr,0.0,p.cut);
(p.W)->set_values(nr,0.0,p.cut,singletable,epsilon); (p.W)->set_values(nr,0.0,p.cut,singletable,epsilon);
} }
for (int i = 0; i < npair; i++) { // P
cutmax = 0.0;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i]; PairParameters & p = pairParameters[i];
if (comm->me == 0) { if (p.cut > cutmax) cutmax = p.cut;
grab(fp,nr,singletable); }
cutmaxsq = cutmax*cutmax;
if (eta != 3) {
for (int j = 0; j < nelements; j++) { // P
if (comm->me == 0) {
grab(fp,nr,singletable);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters & p = tripletParameters[i*nelements*nelements+j*nelements+j];
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
}
}
for (int j = 0; j < nelements-1; j++) { // P
for (int k = j+1; k < nelements; k++) {
if (comm->me == 0) {
grab(fp,nr,singletable);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters & p = tripletParameters[i*nelements*nelements+j*nelements+k];
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
TripletParameters & q = tripletParameters[i*nelements*nelements+k*nelements+j];
q.P = new tabularFunction(nr,-cutmax,cutmax);
(q.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
}
}
}
}
if (eta == 3) {
for (int i = 0; i < ntriple; i++) { // P
TripletParameters & p = tripletParameters[i];
if (comm->me == 0) {
grab(fp,nr,singletable);
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
} }
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
} }
delete[] singletable; delete[] singletable;
singletable = new double[ng]; singletable = new double[ng];
@ -730,6 +762,22 @@ void PairPolymorphic::read_file(char *file)
fclose(fp); fclose(fp);
} }
// recalculate cutoffs of all params
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
p.cut = (p.U)->get_xmax();
if (p.cut < (p.V)->get_xmax()) p.cut = (p.V)->get_xmax();
if (p.cut < (p.W)->get_xmax()) p.cut = (p.W)->get_xmax();
p.cutsq = p.cut*p.cut;
}
// set cutmax to max of all params
cutmax = 0.0;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
if (cutmax < p.cut) cutmax = p.cut;
}
cutmaxsq = cutmax*cutmax;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -771,15 +819,16 @@ void PairPolymorphic::setup_params()
// for debugging, call write_tables() to check the tabular functions // for debugging, call write_tables() to check the tabular functions
// if (comm->me == 0) { // if (comm->me == 0) {
// write_tables(51); // write_tables(51);
// errorX->all(FLERR,"Test potential tables");
// } // }
// error->all(FLERR,"Test potential tables");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
attractive term attractive term
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairPolymorphic::attractive(PairParameters *p, TripletParameters *trip, void PairPolymorphic::attractive(PairParameters *p, PairParameters *q,
TripletParameters *trip,
double prefactor, double rij, double rik, double prefactor, double rij, double rik,
double *delrij, double *delrik, double *delrij, double *delrik,
double *fi, double *fj, double *fk) double *fi, double *fj, double *fk)
@ -793,7 +842,7 @@ void PairPolymorphic::attractive(PairParameters *p, TripletParameters *trip,
rikinv = 1.0/rik; rikinv = 1.0/rik;
vec3_scale(rikinv,delrik,rik_hat); vec3_scale(rikinv,delrik,rik_hat);
ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,p,trip); ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,p,q,trip);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -802,46 +851,39 @@ void PairPolymorphic::ters_zetaterm_d(double prefactor,
double *rij_hat, double rij, double *rij_hat, double rij,
double *rik_hat, double rik, double *rik_hat, double rik,
double *dri, double *drj, double *drk, double *dri, double *drj, double *drk,
PairParameters *p, TripletParameters *trip) PairParameters *p, PairParameters *q,
TripletParameters *trip)
{ {
double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta; double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta;
double dcosdri[3],dcosdrj[3],dcosdrk[3]; double dcosdri[3],dcosdrj[3],dcosdrk[3];
cos_theta = vec3_dot(rij_hat,rik_hat); cos_theta = vec3_dot(rij_hat,rik_hat);
(p->W)->value(rik,fc,1,dfc,1); (q->W)->value(rik,fc,1,dfc,1);
(p->P)->value(rij-(p->xi)*rik,ex_delr,1,ex_delr_d,1); (trip->P)->value(rij-(p->xi)*rik,ex_delr,1,ex_delr_d,1);
(trip->G)->value(cos_theta,gijk,1,gijk_d,1); (trip->G)->value(cos_theta,gijk,1,gijk_d,1);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Ri // compute the derivative wrt Ri
// dri = -dfc*gijk*ex_delr*rik_hat;
// dri += fc*gijk_d*ex_delr*dcosdri;
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri); vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri); vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri); vec3_scaleadd(fc*gijk*ex_delr_d*(p->xi),rik_hat,dri,dri);
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri); vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
vec3_scale(prefactor,dri,dri); vec3_scale(prefactor,dri,dri);
// compute the derivative wrt Rj // compute the derivative wrt Rj
// drj = fc*gijk_d*ex_delr*dcosdrj;
// drj += fc*gijk*ex_delr_d*rij_hat;
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj); vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj); vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
vec3_scale(prefactor,drj,drj); vec3_scale(prefactor,drj,drj);
// compute the derivative wrt Rk // compute the derivative wrt Rk
// drk = dfc*gijk*ex_delr*rik_hat;
// drk += fc*gijk_d*ex_delr*dcosdrk;
// drk += -fc*gijk*ex_delr_d*rik_hat;
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk); vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk); vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk); vec3_scaleadd(-fc*gijk*ex_delr_d*(p->xi),rik_hat,drk,drk);
vec3_scale(prefactor,drk,drk); vec3_scale(prefactor,drk,drk);
} }
@ -919,26 +961,29 @@ void PairPolymorphic::write_tables(int npts)
} }
for (int i = 0; i < nelements; i++) { for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) { for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
strcpy(line,elements[i]); strcpy(line,elements[i]);
strcat(line,elements[j]); strcat(line,elements[j]);
strcat(line,elements[k]);
strcat(line,"_P"); strcat(line,"_P");
strcat(line,tag); strcat(line,tag);
fp = fopen(line, "w"); fp = fopen(line, "w");
int iparam_ij = elem2param[i][j]; int iparam_ij = elem3param[i][j][k];
PairParameters & pair = pairParameters[iparam_ij]; TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.P)->get_xmin(); xmin = (pair.P)->get_xmin();
xmax = (pair.P)->get_xmax(); xmax = (pair.P)->get_xmax();
double xl = xmax - xmin; double xl = xmax - xmin;
xmin = xmin - 0.5*xl; xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl; xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) { for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * k / (npts-1); x = xmin + (xmax-xmin) * n / (npts-1);
(pair.P)->value(x, pf, 1, pfp, 1); (pair.P)->value(x, pf, 1, pfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp); fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp);
} }
fclose(fp); fclose(fp);
} }
} }
}
for (int i = 0; i < nelements; i++) { for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) { for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) { for (int k = 0; k < nelements; k++) {

87
src/MANYBODY/pair_polymorphic.h
View File

@ -104,32 +104,59 @@ class PairPolymorphic : public Pair {
} }
void set_values(int n, double x1, double x2, double * values, double epsilon) void set_values(int n, double x1, double x2, double * values, double epsilon)
{ {
int i0; int shrink = 1;
i0 = n-1; int ilo,ihi;
// shrink (remove near zero points) reduces cutoff radius, and therefore computational cost double vlo,vhi;
// do not shrink when x2 < 1.1 (angular function) or x2 > 20.0 (non-radial function) ilo = 0;
if (x2 >= 1.1 && x2 <= 20.0) { ihi = n-1;
for (int i = n-1; i >= 0; i--) { for (int i = 0; i < n; i++) {
if (fabs(values[i]) > epsilon) { if (fabs(values[i]) <= epsilon) {
i0 = i; ilo = i;
break; } else {
} break;
} }
} }
// do not shrink when when list is abnormally small for (int i = n-1; i >= 0; i--) {
if (i0 < 10/n) { if (fabs(values[i]) <= epsilon) {
i0 = n-1; ihi = i;
} else if (i0 < n-1) { } else {
values[i0] = 0.0; break;
i0 = i0 + 1; }
values[i0] = 0.0;
} }
xmin = x1; if (ihi < ilo) ihi = ilo;
xmax = x1 + (x2-x1)/(n -1)*i0; vlo = values[ilo];
vhi = values[ilo];
for (int i = ilo; i <= ihi; i++) {
if (vlo > values[i]) vlo = values[i];
if (vhi < values[i]) vhi = values[i];
}
// shrink (remove near zero points) reduces cutoff radius, and therefore computational cost
// do not shrink when x2 < 1.1 (angular function) or x2 > 20.0 (non-radial function)
if (x2 < 1.1 || x2 > 20.0) {
shrink = 0;
}
// do not shrink when when list is abnormally small
if (ihi - ilo < 50) {
shrink = 0;
}
// shrink if it is a constant
if (vhi - vlo <= epsilon) {
// shrink = 1;
}
if (shrink == 0) {
ilo = 0;
ihi = n-1;
}
xmin = x1 + (x2-x1)/(n -1)*ilo;
xmax = xmin + (x2-x1)/(n -1)*(ihi-ilo);
xmaxsq = xmax*xmax; xmaxsq = xmax*xmax;
n = i0+1; n = ihi - ilo + 1;
resize(n); resize(n);
memcpy(ys,values,n*sizeof(double)); for (int i = ilo; i <= ihi; i++) {
ys[i-ilo] = values[i];
}
initialize(); initialize();
} }
void value(double x, double &y, int ny, double &y1, int ny1) void value(double x, double &y, int ny, double &y1, int ny1)
@ -228,32 +255,32 @@ class PairPolymorphic : public Pair {
struct PairParameters { struct PairParameters {
double cut; double cut;
double cutsq; double cutsq;
bool xi; // "indicator" double xi;
class tabularFunction * U; class tabularFunction * U;
class tabularFunction * V; class tabularFunction * V;
class tabularFunction * W; class tabularFunction * W;
class tabularFunction * P;
class tabularFunction * F; class tabularFunction * F;
PairParameters() { PairParameters() {
cut = 0.0; cut = 0.0;
cutsq = 0.0; cutsq = 0.0;
xi = true; xi = 1.0;
U = NULL; U = NULL;
V = NULL; V = NULL;
W = NULL; W = NULL;
P = NULL;
F = NULL; F = NULL;
}; };
}; };
struct TripletParameters { struct TripletParameters {
class tabularFunction * P;
class tabularFunction * G; class tabularFunction * G;
TripletParameters() { TripletParameters() {
P = NULL;
G = NULL; G = NULL;
}; };
}; };
double epsilon; double epsilon;
bool eta; // global indicator int eta;
int nx,nr,ng; // table sizes int nx,nr,ng; // table sizes
double maxX; double maxX;
@ -283,11 +310,13 @@ class PairPolymorphic : public Pair {
void setup_params(); void setup_params();
void write_tables(int); void write_tables(int);
void attractive(PairParameters *, TripletParameters *, double, double, void attractive(PairParameters *, PairParameters *, TripletParameters *,
double, double *, double *, double *, double *, double *); double, double, double, double *, double *, double *,
double *, double *);
void ters_zetaterm_d(double, double *, double, double *, double, double *, void ters_zetaterm_d(double, double *, double, double *, double, double *,
double *, double *, PairParameters *, TripletParameters *); double *, double *, PairParameters *, PairParameters *,
TripletParameters *);
void costheta_d(double *, double, double *, double, void costheta_d(double *, double, double *, double,
double *, double *, double *); double *, double *, double *);

4
src/USER-MISC/dihedral_nharmonic.cpp
View File

@ -261,8 +261,8 @@ void DihedralNHarmonic::allocate()
allocated = 1; allocated = 1;
int n = atom->ndihedraltypes; int n = atom->ndihedraltypes;
memory->create(nterms,n+1,"dihedral:nt"); nterms = new int[n+1];
a = new double * [n+1]; a = new double *[n+1];
memory->create(setflag,n+1,"dihedral:setflag"); memory->create(setflag,n+1,"dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0; for (int i = 1; i <= n; i++) setflag[i] = 0;