ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2824 from ndtrung81/dielectric-updates

Browse Source 
Updates to the DIELECTRIC package
This commit is contained in:
Axel Kohlmeyer
2021-07-06 20:58:23 -04:00
committed by GitHub
parent 18f5ef1a76 91790e8586
commit c9ddee5e38
15 changed files with 7565 additions and 6847 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_kspace.rst
View File

@ -33,8 +33,10 @@ OPT.
* :doc:`pppm/cg (o) <kspace_style>`
* :doc:`pppm/dipole <kspace_style>`
* :doc:`pppm/dipole/spin <kspace_style>`
* :doc:`pppm/dielectric <kspace_style>`
* :doc:`pppm/disp (io) <kspace_style>`
* :doc:`pppm/disp/tip4p (o) <kspace_style>`
* :doc:`pppm/disp/dielectric <kspace_style>`
* :doc:`pppm/stagger <kspace_style>`
* :doc:`pppm/tip4p (o) <kspace_style>`
* :doc:`pppm/dielectric <kspace_style>`

7
doc/src/compute_efield_atom.rst
View File

@ -21,6 +21,13 @@ Examples
compute 1 all efield/atom
compute 1 all efield/atom pair yes kspace no
Used in input scripts:
.. parsed-literal::
examples/PACKAGES/dielectric/in.confined
examples/PACKAGES/dielectric/in.nopbc
Description
"""""""""""

12
doc/src/fix_polarize.rst
View File

@ -29,9 +29,17 @@ Examples
.. code-block:: LAMMPS
fix 2 all polarize/bem/gmres 5 0.0001
fix 2 interface polarize/bem/gmres 5 0.0001
fix 1 interface polarize/bem/icc 1 0.0001
fix 3 all polarize/functional 1 0.001
fix 3 interface polarize/functional 1 0.001
Used in input scripts:
.. parsed-literal::
examples/PACKAGES/dielectric/in.confined
examples/PACKAGES/dielectric/in.nopbc
Description
"""""""""""

13
doc/src/kspace_style.rst
View File

@ -9,6 +9,7 @@
.. index:: kspace_style pppm/gpu
.. index:: kspace_style pppm/intel
.. index:: kspace_style pppm/cg
.. index:: kspace_style pppm/dielectric
.. index:: kspace_style pppm/dipole
.. index:: kspace_style pppm/dipole/spin
.. index:: kspace_style pppm/disp
@ -16,7 +17,7 @@
.. index:: kspace_style pppm/disp/tip4p
.. index:: kspace_style pppm/disp/tip4p/omp
.. index:: kspace_style pppm/disp/intel
.. index:: kspace_style pppm/dielectric
.. index:: kspace_style pppm/disp/dielectric
.. index:: kspace_style pppm/cg/omp
.. index:: kspace_style pppm/stagger
.. index:: kspace_style pppm/tip4p
@ -38,7 +39,7 @@ Syntax
kspace_style style value
* style = *none* or *ewald* or *ewald/dipole* or *ewald/dipole/spin* or *ewald/disp* or *ewald/omp* or *pppm* or *pppm/cg* or *pppm/disp* or *pppm/tip4p* or *pppm/stagger* or *pppm/disp/tip4p* or *pppm/gpu* or *pppm/intel* or *pppm/disp/intel* or *pppm/kk* or *pppm/omp* or *pppm/cg/omp* or *pppm/disp/tip4p/omp* or *pppm/tip4p/omp* or *pppm/dielectic* or *msm* or *msm/cg* or *msm/omp* or *msm/cg/omp* or *msm/dielectric* or *scafacos*
* style = *none* or *ewald* or *ewald/dipole* or *ewald/dipole/spin* or *ewald/disp* or *ewald/omp* or *pppm* or *pppm/cg* or *pppm/disp* or *pppm/tip4p* or *pppm/stagger* or *pppm/disp/tip4p* or *pppm/gpu* or *pppm/intel* or *pppm/disp/intel* or *pppm/kk* or *pppm/omp* or *pppm/cg/omp* or *pppm/disp/tip4p/omp* or *pppm/tip4p/omp* or *pppm/dielectic* or *pppm/disp/dielectric* or *msm* or *msm/cg* or *msm/omp* or *msm/cg/omp* or *msm/dielectric* or *scafacos*
.. parsed-literal::
@ -91,6 +92,8 @@ Syntax
accuracy = desired relative error in forces
*pppm/dielectric* value = accuracy
accuracy = desired relative error in forces
*pppm/disp/dielectric* value = accuracy
accuracy = desired relative error in forces
*msm* value = accuracy
accuracy = desired relative error in forces
*msm/cg* value = accuracy (smallq)
@ -118,6 +121,12 @@ Examples
kspace style scafacos fmm 1.0e-4
kspace_style none
Used in input scripts:
.. parsed-literal::
examples/peptide/in.peptide
Description
"""""""""""

296
doc/src/pair_dielectric.rst
View File

@ -36,263 +36,53 @@ pair_style lj/long/coul/long/dielectric command
Syntax
""""""
TODO FIX the rest of the file
.. code-block:: LAMMPS
pair_style lj/cut/dipole/cut cutoff (cutoff2)
pair_style lj/sf/dipole/sf cutoff (cutoff2)
pair_style lj/cut/dipole/long cutoff (cutoff2)
pair_style lj/long/dipole/long flag_lj flag_coul cutoff (cutoff2)
pair_style style args
* cutoff = global cutoff LJ (and Coulombic if only 1 arg) (distance units)
* cutoff2 = global cutoff for Coulombic and dipole (optional) (distance units)
* flag_lj = *long* or *cut* or *off*
.. parsed-literal::
*long* = use long-range damping on dispersion 1/r\^6 term
*cut* = use a cutoff on dispersion 1/r\^6 term
*off* = omit disperion 1/r\^6 term entirely
* flag_coul = *long* or *off*
.. parsed-literal::
*long* = use long-range damping on Coulombic 1/r and point-dipole terms
*off* = omit Coulombic and point-dipole terms entirely
* style = *lj/cut/coul/cut/dielectric* or *lj/cut/coul/long/dielectric* or *lj/cut/coul/msm/dielectric* or *lj/long/coul/msm/dielectric*
* args = list of arguments for a particular style
Examples
""""""""
.. code-block:: LAMMPS
pair_style lj/cut/dipole/cut 10.0
pair_coeff * * 1.0 1.0
pair_coeff 2 3 1.0 1.0 2.5 4.0
pair_style coul/cut/dielectric 10.0
pair_coeff * *
pair_coeff 1 1 9.0
pair_style lj/sf/dipole/sf 9.0
pair_coeff * * 1.0 1.0
pair_coeff 2 3 1.0 1.0 2.5 4.0 scale 0.5
pair_coeff 2 3 1.0 1.0 2.5 4.0
pair_style lj/cut/coul/cut/dielectric 10.0
pair_style lj/cut/coul/cut/dielectric 10.0 8.0
pair_coeff * * 100.0 3.0
pair_coeff 1 1 100.0 3.5 9.0
pair_style lj/cut/dipole/long 10.0
pair_coeff * * 1.0 1.0
pair_coeff 2 3 1.0 1.0 2.5 4.0
pair_style lj/cut/coul/long/dielectric 10.0
pair_style lj/cut/coul/long/dielectric 10.0 8.0
pair_coeff * * 100.0 3.0
pair_coeff 1 1 100.0 3.5 9.0
pair_style lj/long/dipole/long long long 3.5 10.0
pair_coeff * * 1.0 1.0
pair_coeff 2 3 1.0 1.0 2.5 4.0
Used in input scripts:
.. parsed-literal::
examples/PACKAGES/dielectric/in.confined
examples/PACKAGES/dielectric/in.nopbc
Description
"""""""""""
Style *lj/cut/dipole/cut* computes interactions between pairs of particles
that each have a charge and/or a point dipole moment. In addition to
the usual Lennard-Jones interaction between the particles (Elj) the
charge-charge (Eqq), charge-dipole (Eqp), and dipole-dipole (Epp)
interactions are computed by these formulas for the energy (E), force
(F), and torque (T) between particles I and J.
All these pair styles are derived from the corresponding pair styles
without the *dielectric*\ suffix. In addition to computing atom forces
and energies, these pair styles compute the electrical field vector
at each atom, which are to be used in the :doc:`fix polarize <fix_polarize>` commands.
.. math::
These pair styles should be used with :doc:`atom_style dielectric <atom_style>`,
which uses atom charges rescaled by their local dielectric constant.
E_{LJ} = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
\left(\frac{\sigma}{r}\right)^6 \right] \\
E_{qq} = & \frac{q_i q_j}{r} \\
E_{qp} = & \frac{q}{r^3} (p \bullet \vec{r}) \\
E_{pp} = & \frac{1}{r^3} (\vec{p_i} \bullet \vec{p_j}) -
\frac{3}{r^5} (\vec{p_i} \bullet \vec{r}) (\vec{p_j} \bullet \vec{r}) \\
& \\
F_{qq} = & \frac{q_i q_j}{r^3} \vec{r} \\
F_{qp} = & -\frac{q}{r^3} \vec{p} + \frac{3q}{r^5}
(\vec{p} \bullet \vec{r}) \vec{r} \\
F_{pp} = & \frac{3}{r^5} (\vec{p_i} \bullet \vec{p_j}) \vec{r} -
\frac{15}{r^7} (\vec{p_i} \bullet \vec{r})
(\vec{p_j} \bullet \vec{r}) \vec{r} +
\frac{3}{r^5} \left[ (\vec{p_j} \bullet \vec{r}) \vec{p_i} +
(\vec{p_i} \bullet \vec{r}) \vec{p_j} \right] \\
& \\
T_{pq} = T_{ij} = & \frac{q_j}{r^3} (\vec{p_i} \times \vec{r}) \\
T_{qp} = T_{ji} = & - \frac{q_i}{r^3} (\vec{p_j} \times \vec{r}) \\
T_{pp} = T_{ij} = & -\frac{1}{r^3} (\vec{p_i} \times \vec{p_j}) +
\frac{3}{r^5} (\vec{p_j} \bullet \vec{r})
(\vec{p_i} \times \vec{r}) \\
T_{pp} = T_{ji} = & -\frac{1}{r^3} (\vec{p_j} \times \vec{p_i}) +
\frac{3}{r^5} (\vec{p_i} \bullet \vec{r})
(\vec{p_j} \times \vec{r})
where :math:`q_i` and :math:`q_j` are the charges on the two particles,
:math:`\vec{p_i}` and :math:`\vec{p_j}` are the dipole moment vectors of
the two particles, r is their separation distance, and the vector r =
Ri - Rj is the separation vector between the two particles. Note that
Eqq and Fqq are simply Coulombic energy and force, Fij = -Fji as
symmetric forces, and Tij != -Tji since the torques do not act
symmetrically. These formulas are discussed in :ref:`(Allen) <Allen2>`
and in :ref:`(Toukmaji) <Toukmaji2>`.
Also note, that in the code, all of these terms (except Elj) have a
:math:`C/\epsilon` prefactor, the same as the Coulombic term in the LJ +
Coulombic pair styles discussed :doc:`here <pair_lj>`. C is an
energy-conversion constant and epsilon is the dielectric constant
which can be set by the :doc:`dielectric <dielectric>` command. The
same is true of the equations that follow for other dipole pair
styles.
Style *lj/sf/dipole/sf* computes "shifted-force" interactions between
pairs of particles that each have a charge and/or a point dipole
moment. In general, a shifted-force potential is a (slightly) modified
potential containing extra terms that make both the energy and its
derivative go to zero at the cutoff distance; this removes
(cutoff-related) problems in energy conservation and any numerical
instability in the equations of motion :ref:`(Allen) <Allen2>`. Shifted-force
interactions for the Lennard-Jones (E_LJ), charge-charge (Eqq),
charge-dipole (Eqp), dipole-charge (Epq) and dipole-dipole (Epp)
potentials are computed by these formulas for the energy (E), force
(F), and torque (T) between particles I and J:
.. math::
E_{LJ} = & 4\epsilon \left\{ \left[ \left( \frac{\sigma}{r} \right)^{\!12} -
\left( \frac{\sigma}{r} \right)^{\!6} \right] +
\left[ 6\left( \frac{\sigma}{r_c} \right)^{\!12} -
3\left(\frac{\sigma}{r_c}\right)^{\!6}\right]\left(\frac{r}{r_c}\right)^{\!2}
- 7\left( \frac{\sigma}{r_c} \right)^{\!12} +
4\left( \frac{\sigma}{r_c} \right)^{\!6}\right\} \\
E_{qq} = & \frac{q_i q_j}{r}\left(1-\frac{r}{r_c}\right)^{\!2} \\
E_{pq} = & E_{ji} = -\frac{q}{r^3} \left[ 1 -
3\left(\frac{r}{r_c}\right)^{\!2} +
2\left(\frac{r}{r_c}\right)^{\!3}\right] (\vec{p}\bullet\vec{r}) \\
E_{qp} = & E_{ij} = \frac{q}{r^3} \left[ 1 -
3\left(\frac{r}{r_c}\right)^{\!2} +
2\left(\frac{r}{r_c}\right)^{\!3}\right] (\vec{p}\bullet\vec{r}) \\
E_{pp} = & \left[1-4\left(\frac{r}{r_c}\right)^{\!3} +
3\left(\frac{r}{r_c}\right)^{\!4}\right]\left[\frac{1}{r^3}
(\vec{p_i} \bullet \vec{p_j}) - \frac{3}{r^5}
(\vec{p_i} \bullet \vec{r}) (\vec{p_j} \bullet \vec{r})\right] \\
& \\
F_{LJ} = & \left\{\left[48\epsilon \left(\frac{\sigma}{r}\right)^{\!12} -
24\epsilon \left(\frac{\sigma}{r}\right)^{\!6} \right]\frac{1}{r^2} -
\left[48\epsilon \left(\frac{\sigma}{r_c}\right)^{\!12} - 24\epsilon
\left(\frac{\sigma}{r_c}\right)^{\!6} \right]\frac{1}{r_c^2}\right\}\vec{r}\\
F_{qq} = & \frac{q_i q_j}{r}\left(\frac{1}{r^2} -
\frac{1}{r_c^2}\right)\vec{r} \\
F_{pq} = & F_{ij } = -\frac{3q}{r^5} \left[ 1 -
\left(\frac{r}{r_c}\right)^{\!2}\right](\vec{p}\bullet\vec{r})\vec{r} +
\frac{q}{r^3}\left[1-3\left(\frac{r}{r_c}\right)^{\!2} +
2\left(\frac{r}{r_c}\right)^{\!3}\right] \vec{p} \\
F_{qp} = & F_{ij} = \frac{3q}{r^5} \left[ 1 -
\left(\frac{r}{r_c}\right)^{\!2}\right] (\vec{p}\bullet\vec{r})\vec{r} -
\frac{q}{r^3}\left[1-3\left(\frac{r}{r_c}\right)^{\!2} +
2\left(\frac{r}{r_c}\right)^{\!3}\right] \vec{p} \\
F_{pp} = &\frac{3}{r^5}\Bigg\{\left[1-\left(\frac{r}{r_c}\right)^{\!4}\right]
\left[(\vec{p_i}\bullet\vec{p_j}) - \frac{3}{r^2} (\vec{p_i}\bullet\vec{r})
(\vec{p_j} \bullet \vec{r})\right] \vec{r} + \\
& \left[1 -
4\left(\frac{r}{r_c}\right)^{\!3}+3\left(\frac{r}{r_c}\right)^{\!4}\right]
\left[ (\vec{p_j} \bullet \vec{r}) \vec{p_i} + (\vec{p_i} \bullet \vec{r})
\vec{p_j} -\frac{2}{r^2} (\vec{p_i} \bullet \vec{r})
(\vec{p_j} \bullet \vec{r})\vec{r}\right] \Bigg\}
.. math::
T_{pq} = T_{ij} = & \frac{q_j}{r^3} \left[ 1 -
3\left(\frac{r}{r_c}\right)^{\!2} +
2\left(\frac{r}{r_c}\right)^{\!3}\right] (\vec{p_i}\times\vec{r}) \\
T_{qp} = T_{ji} = & - \frac{q_i}{r^3} \left[ 1 -
3\left(\frac{r}{r_c}\right)^{\!2} +
2\left(\frac{r}{r_c}\right)^{\!3} \right] (\vec{p_j}\times\vec{r}) \\
T_{pp} = T_{ij} = & -\frac{1}{r^3}\left[1-4\left(\frac{r}{r_c}\right)^{\!3} +
e3\left(\frac{r}{r_c}\right)^{\!4}\right] (\vec{p_i} \times \vec{p_j}) + \\
& \frac{3}{r^5}\left[1-4\left(\frac{r}{r_c}\right)^{\!3} +
3\left(\frac{r}{r_c}\right)^{\!4}\right] (\vec{p_j}\bullet\vec{r})
(\vec{p_i} \times \vec{r}) \\
T_{pp} = T_{ji} = & -\frac{1}{r^3}\left[1-4\left(\frac{r}{r_c}\right)^{\!3} +
3\left(\frac{r}{r_c}\right)^{\!4}\right](\vec{p_j} \times \vec{p_i}) + \\
& \frac{3}{r^5}\left[1-4\left(\frac{r}{r_c}\right)^{\!3} +
3\left(\frac{r}{r_c}\right)^{\!4}\right] (\vec{p_i} \bullet \vec{r})
(\vec{p_j} \times \vec{r})
where :math:`\epsilon` and :math:`\sigma` are the standard LJ
parameters, :math:`r_c` is the cutoff, :math:`q_i` and :math:`q_j` are
the charges on the two particles, :math:`\vec{p_i}` and
:math:`\vec{p_j}` are the dipole moment vectors of the two particles, r
is their separation distance, and the vector r = Ri - Rj is the
separation vector between the two particles. Note that Eqq and Fqq are
simply Coulombic energy and force, Fij = -Fji as symmetric forces, and
Tij != -Tji since the torques do not act symmetrically. The
shifted-force formula for the Lennard-Jones potential is reported in
:ref:`(Stoddard) <Stoddard>`. The original (non-shifted) formulas for
the electrostatic potentials, forces and torques can be found in
:ref:`(Price) <Price2>`. The shifted-force electrostatic potentials have
been obtained by applying equation 5.13 of :ref:`(Allen) <Allen2>`. The
formulas for the corresponding forces and torques have been obtained by
applying the 'chain rule' as in appendix C.3 of :ref:`(Allen) <Allen2>`.
If one cutoff is specified in the pair_style command, it is used for
both the LJ and Coulombic (q,p) terms. If two cutoffs are specified,
they are used as cutoffs for the LJ and Coulombic (q,p) terms
respectively. This pair style also supports an optional *scale* keyword
as part of a pair_coeff statement, where the interactions can be
scaled according to this factor. This scale factor is also made available
for use with fix adapt.
Style *lj/cut/dipole/long* computes long-range point-dipole
interactions as discussed in :ref:`(Toukmaji) <Toukmaji2>`. Dipole-dipole,
dipole-charge, and charge-charge interactions are all supported, along
with the standard 12/6 Lennard-Jones interactions, which are computed
with a cutoff. A :doc:`kspace_style <kspace_style>` must be defined to
use this pair style. Currently, only :doc:`kspace_style ewald/disp <kspace_style>` support long-range point-dipole
interactions.
Style *lj/long/dipole/long* also computes point-dipole interactions as
discussed in :ref:`(Toukmaji) <Toukmaji2>`. Long-range dipole-dipole,
dipole-charge, and charge-charge interactions are all supported, along
with the standard 12/6 Lennard-Jones interactions. LJ interactions
can be cutoff or long-ranged.
For style *lj/long/dipole/long*\ , if *flag_lj* is set to *long*\ , no
cutoff is used on the LJ 1/r\^6 dispersion term. The long-range
portion is calculated by using the :doc:`kspace_style ewald_disp <kspace_style>` command. The specified LJ cutoff then
determines which portion of the LJ interactions are computed directly
by the pair potential versus which part is computed in reciprocal
space via the Kspace style. If *flag_lj* is set to *cut*\ , the LJ
interactions are simply cutoff, as with :doc:`pair_style lj/cut <pair_lj>`. If *flag_lj* is set to *off*\ , LJ interactions
are not computed at all.
If *flag_coul* is set to *long*\ , no cutoff is used on the Coulombic or
dipole interactions. The long-range portion is calculated by using
*ewald_disp* of the :doc:`kspace_style <kspace_style>` command. If
*flag_coul* is set to *off*\ , Coulombic and dipole interactions are not
computed at all.
Atoms with dipole moments should be integrated using the :doc:`fix nve/sphere update dipole <fix_nve_sphere>` or the :doc:`fix nvt/sphere update dipole <fix_nvt_sphere>` command to rotate the
dipole moments. The *omega* option on the :doc:`fix langevin <fix_langevin>` command can be used to thermostat the
rotational motion. The :doc:`compute temp/sphere <compute_temp_sphere>`
command can be used to monitor the temperature, since it includes
rotational degrees of freedom. The :doc:`atom_style hybrid dipole sphere <atom_style>` command should be used since
it defines the point dipoles and their rotational state.
The magnitude and orientation of the dipole moment for each particle
can be defined by the :doc:`set <set>` command or in the "Atoms" section
of the data file read in by the :doc:`read_data <read_data>` command.
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 read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands, or by mixing as described below:
* :math:`\epsilon` (energy units)
* :math:`\sigma` (distance units)
* cutoff1 (distance units)
* cutoff2 (distance units)
The latter 2 coefficients are optional. If not specified, the global
LJ and Coulombic cutoffs specified in the pair_style command are used.
If only one cutoff is specified, it is used as the cutoff for both LJ
and Coulombic interactions for this type pair. If both coefficients
are specified, they are used as the LJ and Coulombic cutoffs for this
type pair.
The styles lj/cut/coul/long/dielectric, lj/cut/coul/msm/dielectric, and
lj/long/coul/long/dielectric should be used with their kspace style counterparts,
namely, pppm/dielectric, pppm/disp/dielectric, and msm/dielectric, respectively.
----------
@ -307,23 +97,9 @@ For atom type pairs I,J and I != J, the epsilon and sigma coefficients
and cutoff distances for this pair style can be mixed. The default
mix value is *geometric*\ . See the "pair_modify" command for details.
For atom type pairs I,J and I != J, the A, sigma, d1, and d2
coefficients and cutoff distance for this pair style can be mixed. A
is an energy value mixed like a LJ epsilon. D1 and d2 are distance
values and are mixed like sigma. The default mix value is
*geometric*\ . See the "pair_modify" command for details.
This pair style does not support the :doc:`pair_modify <pair_modify>`
shift option for the energy of the Lennard-Jones portion of the pair
interaction; such energy goes to zero at the cutoff by construction.
The :doc:`pair_modify <pair_modify>` table option is not relevant
for this pair style.
This pair style does not support the :doc:`pair_modify <pair_modify>`
tail option for adding long-range tail corrections to energy and
pressure.
This pair style writes its information to :doc:`binary restart files <restart>`, so pair_style and pair_coeff commands do not need
to be specified in an input script that reads a restart file.
@ -334,21 +110,13 @@ This pair style can only be used via the *pair* keyword of the
Restrictions
""""""""""""
The *lj/cut/dipole/cut*\ , *lj/cut/dipole/long*\ , and
*lj/long/dipole/long* styles are part of the DIPOLE package. They are
only enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
The *lj/sf/dipole/sf* style is part of the USER-MISC package. It is
only enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Using dipole pair styles with *electron* :doc:`units <units>` is not
currently supported.
These styles are part of the DIELECTRIC package. They are only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`, :doc:`set <set>`, :doc:`read_data <read_data>`,
:doc:`fix nve/sphere <fix_nve_sphere>`, :doc:`fix nvt/sphere <fix_nvt_sphere>`
:doc:`pair_coeff <pair_coeff>`, :doc:`fix polarize <fix_polarize>`, :doc:`read_data <read_data>`
Default
"""""""

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

@ -2407,6 +2407,7 @@ orthorhombic
Ortner
oso
Otype
Ouadfel
Ouldridge
outfile
outmost
@ -2680,7 +2681,6 @@ qoverride
qqr
qqrd
qtb
Quadfel
quadratically
quadrupolar
Quant
@ -3155,6 +3155,7 @@ Suter
Sutmann
Suzen
svn
Svoboda
sw
Swegat
swiggle

4
examples/PACKAGES/dielectric/README
View File

@ -31,9 +31,7 @@ where
* epsilon = the local epsilon value at the vertex or at the ion.
For real charges, epsilon is the medium dielectric constant,
and q is the real (unscaled) charges.
For interface particles:
+ if q is zero (zero surface charges), epsilon is set to be 1.0;
+ if q is nonzero (charged surfaces), epsilon is set to be em
For interface particles, epsilon is set to be em
(the mean dielectric value above).
* area_per_patch: the surface area of the patch (element).

8000
examples/PACKAGES/dielectric/data.confined
View File

File diff suppressed because it is too large Load Diff

5124
examples/PACKAGES/dielectric/data.sphere
View File

File diff suppressed because it is too large Load Diff

2
src/DIELECTRIC/Install.sh
View File

@ -31,7 +31,7 @@ action () {
if (test $1 = 1) then
if (test ! -e ../pppm.cpp) then
echo "Must install KSPACE package with DIELECTRIC"
echo "Must install KSPACE package with DIELECTRIC package"
exit 1
fi
fi

3
src/DIELECTRIC/fix_polarize_functional.cpp
View File

@ -1117,9 +1117,6 @@ void FixPolarizeFunctional::calculate_grad_greens_ewald(double *vec, double dx,
void FixPolarizeFunctional::calculate_matrix_multiply_vector(double **matrix, double *in_vec,
double *out_vec, int M)
{
#if defined(OPENMP)
#pragma parallel omp for
#endif
for (int k = 0; k < M; ++k) {
double temp = 0.0;
for (int l = 0; l < M; ++l) { temp += matrix[k][l] * in_vec[l]; }

6
src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp
View File

@ -53,9 +53,9 @@ using namespace MathExtra;
PairLJLongCoulLongDielectric::PairLJLongCoulLongDielectric(LAMMPS *lmp) : PairLJLongCoulLong(lmp)
{
respa_enable = 0;
cut_respa = NULL;
efield = NULL;
epot = NULL;
cut_respa = nullptr;
efield = nullptr;
epot = nullptr;
nmax = 0;
}

10
src/DIELECTRIC/pppm_dielectric.cpp
View File

@ -469,21 +469,19 @@ void PPPMDielectric::slabcorr()
}
/* ----------------------------------------------------------------------
compute qsum,qsqsum,q2 and give error/warning if not charge neutral
called initially, when particle count changes, when charges are changed
compute qsum,qsqsum,q2 and ignore error/warning if not charge neutral
called whenever charges are changed
------------------------------------------------------------------------- */
void PPPMDielectric::qsum_qsq()
{
const double * const q = atom->q;
const double * const eps = atom->epsilon;
const int nlocal = atom->nlocal;
double qsum_local(0.0), qsqsum_local(0.0);
for (int i = 0; i < nlocal; i++) {
double qtmp = eps[i]*q[i];
qsum_local += qtmp;
qsqsum_local += qtmp*qtmp;
qsum_local += q[i];
qsqsum_local += q[i]*q[i];
}
MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);

868
src/DIELECTRIC/pppm_disp_dielectric.cpp Normal file
View File

@ -0,0 +1,868 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Trung Nguyen (Northwestern)
------------------------------------------------------------------------- */
#include "pppm_disp_dielectric.h"
#include "angle.h"
#include "atom.h"
#include "atom_vec_dielectric.h"
#include "bond.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fft3d_wrap.h"
#include "force.h"
#include "gridcomm.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "pair.h"
#include "remap_wrap.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXORDER 7
#define OFFSET 16384
#define SMALL 0.00001
#define LARGE 10000.0
#define EPS_HOC 1.0e-7
enum{REVERSE_RHO,REVERSE_RHO_GEOM,REVERSE_RHO_ARITH,REVERSE_RHO_NONE};
enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM,
FORWARD_IK_GEOM,FORWARD_AD_GEOM,
FORWARD_IK_PERATOM_GEOM,FORWARD_AD_PERATOM_GEOM,
FORWARD_IK_ARITH,FORWARD_AD_ARITH,
FORWARD_IK_PERATOM_ARITH,FORWARD_AD_PERATOM_ARITH,
FORWARD_IK_NONE,FORWARD_AD_NONE,FORWARD_IK_PERATOM_NONE,
FORWARD_AD_PERATOM_NONE};
#ifdef FFT_SINGLE
#define ZEROF 0.0f
#define ONEF 1.0f
#else
#define ZEROF 0.0
#define ONEF 1.0
#endif
/* ---------------------------------------------------------------------- */
PPPMDispDielectric::PPPMDispDielectric(LAMMPS *lmp) : PPPMDisp(lmp)
{
dipoleflag = 0; // turned off for now, until dipole works
group_group_enable = 0;
mu_flag = 0;
efield = nullptr;
phi = nullptr;
potflag = 0;
avec = (AtomVecDielectric *) atom->style_match("dielectric");
if (!avec) error->all(FLERR,"pppm/dielectric requires atom style dielectric");
}
/* ---------------------------------------------------------------------- */
PPPMDispDielectric::~PPPMDispDielectric()
{
memory->destroy(efield);
memory->destroy(phi);
}
/* ----------------------------------------------------------------------
compute the PPPM long-range force, energy, virial
------------------------------------------------------------------------- */
void PPPMDispDielectric::compute(int eflag, int vflag)
{
int i;
// set energy/virial flags
// invoke allocate_peratom() if needed for first time
ev_init(eflag,vflag);
if (evflag_atom && !peratom_allocate_flag) allocate_peratom();
// convert atoms from box to lamda coords
if (triclinic == 0) boxlo = domain->boxlo;
else {
boxlo = domain->boxlo_lamda;
domain->x2lamda(atom->nlocal);
}
// extend size of per-atom arrays if necessary
if (atom->nmax > nmax) {
if (function[0]) {
memory->destroy(part2grid);
memory->destroy(efield);
memory->destroy(phi);
}
if (function[1] + function[2] + function[3]) memory->destroy(part2grid_6);
nmax = atom->nmax;
if (function[0]) {
memory->create(part2grid,nmax,3,"pppm/disp:part2grid");
memory->create(efield,nmax,3,"pppm/disp:efield");
memory->create(phi,nmax,"pppm/disp:phi");
}
if (function[1] + function[2] + function[3])
memory->create(part2grid_6,nmax,3,"pppm/disp:part2grid_6");
}
energy = 0.0;
energy_1 = 0.0;
energy_6 = 0.0;
if (vflag) for (i = 0; i < 6; i++) virial_6[i] = virial_1[i] = 0.0;
// find grid points for all my particles
// distribute partcles' charges/dispersion coefficients on the grid
// communication between processors and remapping two fft
// Solution of poissons equation in k-space and backtransformation
// communication between processors
// calculation of forces
if (function[0]) {
// perform calculations for coulomb interactions only
particle_map_c(delxinv,delyinv,delzinv,shift,part2grid,nupper,nlower,
nxlo_out,nylo_out,nzlo_out,nxhi_out,nyhi_out,nzhi_out);
make_rho_c();
gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
brick2fft(nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
density_brick,density_fft,work1,remap);
if (differentiation_flag == 1) {
poisson_ad(work1,work2,density_fft,fft1,fft2,
nx_pppm,ny_pppm,nz_pppm,nfft,
nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft,
nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
energy_1,greensfn,
virial_1,vg,vg2,
u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick);
gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
fieldforce_c_ad();
if (vflag_atom)
gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
} else {
poisson_ik(work1,work2,density_fft,fft1,fft2,
nx_pppm,ny_pppm,nz_pppm,nfft,
nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft,
nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in,
energy_1,greensfn,
fkx,fky,fkz,fkx2,fky2,fkz2,
vdx_brick,vdy_brick,vdz_brick,virial_1,vg,vg2,
u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick);
gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
fieldforce_c_ik();
if (evflag_atom)
gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM,
gc_buf1,gc_buf2,MPI_FFT_SCALAR);
}
if (evflag_atom) fieldforce_c_peratom();
}
if (function[1]) {
// perform calculations for geometric mixing
particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6,
nupper_6,nlower_6,
nxlo_out_6,nylo_out_6,nzlo_out_6,
nxhi_out_6,nyhi_out_6,nzhi_out_6);
make_rho_g();
gc6->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO_GEOM,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
brick2fft(nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
density_brick_g,density_fft_g,work1_6,remap_6);
if (differentiation_flag == 1) {
poisson_ad(work1_6,work2_6,density_fft_g,fft1_6,fft2_6,
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
energy_6,greensfn_6,
virial_6,vg_6,vg2_6,
u_brick_g,v0_brick_g,v1_brick_g,v2_brick_g,
v3_brick_g,v4_brick_g,v5_brick_g);
gc6->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD_GEOM,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
fieldforce_g_ad();
if (vflag_atom)
gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_GEOM,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
} else {
poisson_ik(work1_6,work2_6,density_fft_g,fft1_6,fft2_6,
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
energy_6,greensfn_6,
fkx_6,fky_6,fkz_6,fkx2_6,fky2_6,fkz2_6,
vdx_brick_g,vdy_brick_g,vdz_brick_g,virial_6,vg_6,vg2_6,
u_brick_g,v0_brick_g,v1_brick_g,v2_brick_g,
v3_brick_g,v4_brick_g,v5_brick_g);
gc6->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK_GEOM,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
fieldforce_g_ik();
if (evflag_atom)
gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_GEOM,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
}
if (evflag_atom) fieldforce_g_peratom();
}
if (function[2]) {
// perform calculations for arithmetic mixing
particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6,
nupper_6,nlower_6,
nxlo_out_6,nylo_out_6,nzlo_out_6,
nxhi_out_6,nyhi_out_6,nzhi_out_6);
make_rho_a();
gc6->reverse_comm_kspace(this,7,sizeof(FFT_SCALAR),REVERSE_RHO_ARITH,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
brick2fft_a();
if (differentiation_flag == 1) {
poisson_ad(work1_6,work2_6,density_fft_a3,fft1_6,fft2_6,
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
energy_6,greensfn_6,
virial_6,vg_6,vg2_6,
u_brick_a3,v0_brick_a3,v1_brick_a3,v2_brick_a3,
v3_brick_a3,v4_brick_a3,v5_brick_a3);
poisson_2s_ad(density_fft_a0,density_fft_a6,
u_brick_a0,v0_brick_a0,v1_brick_a0,v2_brick_a0,
v3_brick_a0,v4_brick_a0,v5_brick_a0,
u_brick_a6,v0_brick_a6,v1_brick_a6,v2_brick_a6,
v3_brick_a6,v4_brick_a6,v5_brick_a6);
poisson_2s_ad(density_fft_a1,density_fft_a5,
u_brick_a1,v0_brick_a1,v1_brick_a1,v2_brick_a1,
v3_brick_a1,v4_brick_a1,v5_brick_a1,
u_brick_a5,v0_brick_a5,v1_brick_a5,v2_brick_a5,
v3_brick_a5,v4_brick_a5,v5_brick_a5);
poisson_2s_ad(density_fft_a2,density_fft_a4,
u_brick_a2,v0_brick_a2,v1_brick_a2,v2_brick_a2,
v3_brick_a2,v4_brick_a2,v5_brick_a2,
u_brick_a4,v0_brick_a4,v1_brick_a4,v2_brick_a4,
v3_brick_a4,v4_brick_a4,v5_brick_a4);
gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_ARITH,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
fieldforce_a_ad();
if (evflag_atom)
gc6->forward_comm_kspace(this,42,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_ARITH,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
} else {
poisson_ik(work1_6,work2_6,density_fft_a3,fft1_6,fft2_6,
nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6,
nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6,
nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6,
energy_6,greensfn_6,
fkx_6,fky_6,fkz_6,fkx2_6,fky2_6,fkz2_6,
vdx_brick_a3,vdy_brick_a3,vdz_brick_a3,virial_6,vg_6,vg2_6,
u_brick_a3,v0_brick_a3,v1_brick_a3,v2_brick_a3,
v3_brick_a3,v4_brick_a3,v5_brick_a3);
poisson_2s_ik(density_fft_a0,density_fft_a6,
vdx_brick_a0,vdy_brick_a0,vdz_brick_a0,
vdx_brick_a6,vdy_brick_a6,vdz_brick_a6,
u_brick_a0,v0_brick_a0,v1_brick_a0,v2_brick_a0,
v3_brick_a0,v4_brick_a0,v5_brick_a0,
u_brick_a6,v0_brick_a6,v1_brick_a6,v2_brick_a6,
v3_brick_a6,v4_brick_a6,v5_brick_a6);
poisson_2s_ik(density_fft_a1,density_fft_a5,
vdx_brick_a1,vdy_brick_a1,vdz_brick_a1,
vdx_brick_a5,vdy_brick_a5,vdz_brick_a5,
u_brick_a1,v0_brick_a1,v1_brick_a1,v2_brick_a1,
v3_brick_a1,v4_brick_a1,v5_brick_a1,
u_brick_a5,v0_brick_a5,v1_brick_a5,v2_brick_a5,
v3_brick_a5,v4_brick_a5,v5_brick_a5);
poisson_2s_ik(density_fft_a2,density_fft_a4,
vdx_brick_a2,vdy_brick_a2,vdz_brick_a2,
vdx_brick_a4,vdy_brick_a4,vdz_brick_a4,
u_brick_a2,v0_brick_a2,v1_brick_a2,v2_brick_a2,
v3_brick_a2,v4_brick_a2,v5_brick_a2,
u_brick_a4,v0_brick_a4,v1_brick_a4,v2_brick_a4,
v3_brick_a4,v4_brick_a4,v5_brick_a4);
gc6->forward_comm_kspace(this,21,sizeof(FFT_SCALAR),FORWARD_IK_ARITH,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
fieldforce_a_ik();
if (evflag_atom)
gc6->forward_comm_kspace(this,49,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_ARITH,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
}
if (evflag_atom) fieldforce_a_peratom();
}
if (function[3]) {
// perform calculations if no mixing rule applies
particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6,
nupper_6,nlower_6,
nxlo_out_6,nylo_out_6,nzlo_out_6,
nxhi_out_6,nyhi_out_6,nzhi_out_6);
make_rho_none();
gc6->reverse_comm_kspace(this,nsplit_alloc,sizeof(FFT_SCALAR),REVERSE_RHO_NONE,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
brick2fft_none();
if (differentiation_flag == 1) {
int n = 0;
for (int k = 0; k < nsplit_alloc/2; k++) {
poisson_none_ad(n,n+1,density_fft_none[n],density_fft_none[n+1],
u_brick_none[n],u_brick_none[n+1],
v0_brick_none,v1_brick_none,v2_brick_none,
v3_brick_none,v4_brick_none,v5_brick_none);
n += 2;
}
gc6->forward_comm_kspace(this,1*nsplit_alloc,sizeof(FFT_SCALAR),
FORWARD_AD_NONE,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
fieldforce_none_ad();
if (vflag_atom)
gc6->forward_comm_kspace(this,6*nsplit_alloc,sizeof(FFT_SCALAR),
FORWARD_AD_PERATOM_NONE,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
} else {
int n = 0;
for (int k = 0; k < nsplit_alloc/2; k++) {
poisson_none_ik(n,n+1,density_fft_none[n],density_fft_none[n+1],
vdx_brick_none[n],vdy_brick_none[n],vdz_brick_none[n],
vdx_brick_none[n+1],vdy_brick_none[n+1],vdz_brick_none[n+1],
u_brick_none,v0_brick_none,v1_brick_none,v2_brick_none,
v3_brick_none,v4_brick_none,v5_brick_none);
n += 2;
}
gc6->forward_comm_kspace(this,3*nsplit_alloc,sizeof(FFT_SCALAR),
FORWARD_IK_NONE,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
fieldforce_none_ik();
if (evflag_atom)
gc6->forward_comm_kspace(this,7*nsplit_alloc,sizeof(FFT_SCALAR),
FORWARD_IK_PERATOM_NONE,
gc6_buf1,gc6_buf2,MPI_FFT_SCALAR);
}
if (evflag_atom) fieldforce_none_peratom();
}
// update qsum and qsqsum, if atom count has changed and energy needed
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum energy across procs and add in volume-dependent term
const double qscale = force->qqrd2e * scale;
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy_1,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy_1 = energy_all;
MPI_Allreduce(&energy_6,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy_6 = energy_all;
energy_1 *= 0.5*volume;
energy_6 *= 0.5*volume;
energy_1 -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij +
1.0/12.0*pow(g_ewald_6,6)*csum;
energy_1 *= qscale;
}
// sum virial across procs
if (vflag_global) {
double virial_all[6];
MPI_Allreduce(virial_1,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i];
MPI_Allreduce(virial_6,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) virial[i] += 0.5*volume*virial_all[i];
if (function[1]+function[2]+function[3]) {
double a = MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij;
virial[0] -= a;
virial[1] -= a;
virial[2] -= a;
}
}
if (eflag_atom) {
if (function[0]) {
double *q = atom->q;
// coulomb self energy correction
for (i = 0; i < atom->nlocal; i++) {
eatom[i] -= qscale*g_ewald*q[i]*q[i]/MY_PIS +
qscale*MY_PI2*q[i]*qsum / (g_ewald*g_ewald*volume);
}
}
if (function[1] + function[2] + function[3]) {
int tmp;
for (i = 0; i < atom->nlocal; i++) {
tmp = atom->type[i];
eatom[i] += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp] +
1.0/12.0*pow(g_ewald_6,6)*cii[tmp];
}
}
}
if (vflag_atom) {
if (function[1] + function[2] + function[3]) {
int tmp;
// dispersion self virial correction
for (i = 0; i < atom->nlocal; i++) {
tmp = atom->type[i];
for (int n = 0; n < 3; n++)
vatom[i][n] -= MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp];
}
}
}
// 2d slab correction
if (slabflag) slabcorr(eflag);
if (function[0]) energy += energy_1;
if (function[1] + function[2] + function[3]) energy += energy_6;
// convert atoms back from lamda to box coords
if (triclinic) domain->lamda2x(atom->nlocal);
}
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles
for ik scheme
------------------------------------------------------------------------- */
void PPPMDispDielectric::fieldforce_c_ik()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz,x0,y0,z0;
FFT_SCALAR ekx,eky,ekz,u;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
double *q = atom->q;
double **x = atom->x;
double **f = atom->f;
double *eps = atom->epsilon;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
u = ekx = eky = ekz = ZEROF;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
z0 = rho1d[2][n];
for (m = nlower; m <= nupper; m++) {
my = m+ny;
y0 = z0*rho1d[1][m];
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
x0 = y0*rho1d[0][l];
if (potflag) u += x0*u_brick[mz][my][mx];
ekx -= x0*vdx_brick[mz][my][mx];
eky -= x0*vdy_brick[mz][my][mx];
ekz -= x0*vdz_brick[mz][my][mx];
}
}
}
// electrostatic potential
if (potflag) phi[i] = u;
// convert E-field to force
const double efactor = scale * eps[i];
efield[i][0] = efactor*ekx;
efield[i][1] = efactor*eky;
efield[i][2] = efactor*ekz;
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i];
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
if (slabflag != 2) f[i][2] += qfactor*ekz;
}
}
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles
for ad scheme
------------------------------------------------------------------------- */
void PPPMDispDielectric::fieldforce_c_ad()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz;
FFT_SCALAR ekx,eky,ekz,u;
double s1,s2,s3;
double sf = 0.0;
double *prd;
if (triclinic == 0) prd = domain->prd;
else prd = domain->prd_lamda;
double xprd = prd[0];
double yprd = prd[1];
double zprd = prd[2];
double zprd_slab = zprd*slab_volfactor;
double hx_inv = nx_pppm/xprd;
double hy_inv = ny_pppm/yprd;
double hz_inv = nz_pppm/zprd_slab;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
double *q = atom->q;
double **x = atom->x;
double **f = atom->f;
double *eps = atom->epsilon;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
compute_drho1d(dx,dy,dz, order, drho_coeff, drho1d);
u = ekx = eky = ekz = ZEROF;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
for (m = nlower; m <= nupper; m++) {
my = m+ny;
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
u += rho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
ekx += drho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
eky += rho1d[0][l]*drho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx];
ekz += rho1d[0][l]*rho1d[1][m]*drho1d[2][n]*u_brick[mz][my][mx];
}
}
}
ekx *= hx_inv;
eky *= hy_inv;
ekz *= hz_inv;
// electrical potential
if (potflag) phi[i] = u;
// convert E-field to force and substract self forces
const double qfactor = qqrd2e * scale;
double qtmp = eps[i]*q[i];
s1 = x[i][0]*hx_inv;
s2 = x[i][1]*hy_inv;
s3 = x[i][2]*hz_inv;
sf = sf_coeff[0]*sin(2*MY_PI*s1);
sf += sf_coeff[1]*sin(4*MY_PI*s1);
sf *= 2*q[i]*q[i];
f[i][0] += qfactor*(ekx*q[i] - sf);
sf = sf_coeff[2]*sin(2*MY_PI*s2);
sf += sf_coeff[3]*sin(4*MY_PI*s2);
sf *= 2*q[i]*q[i];
f[i][1] += qfactor*(eky*q[i] - sf);
sf = sf_coeff[4]*sin(2*MY_PI*s3);
sf += sf_coeff[5]*sin(4*MY_PI*s3);
sf *= 2*q[i]*q[i];
if (slabflag != 2) f[i][2] += qfactor*(ekz*q[i] - sf);
}
}
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles
------------------------------------------------------------------------- */
void PPPMDispDielectric::fieldforce_c_peratom()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz,x0,y0,z0;
FFT_SCALAR u_pa,v0,v1,v2,v3,v4,v5;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
double *q = atom->q;
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
u_pa = v0 = v1 = v2 = v3 = v4 = v5 = ZEROF;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
z0 = rho1d[2][n];
for (m = nlower; m <= nupper; m++) {
my = m+ny;
y0 = z0*rho1d[1][m];
for (l = nlower; l <= nupper; l++) {
mx = l+nx;
x0 = y0*rho1d[0][l];
if (eflag_atom) u_pa += x0*u_brick[mz][my][mx];
if (vflag_atom) {
v0 += x0*v0_brick[mz][my][mx];
v1 += x0*v1_brick[mz][my][mx];
v2 += x0*v2_brick[mz][my][mx];
v3 += x0*v3_brick[mz][my][mx];
v4 += x0*v4_brick[mz][my][mx];
v5 += x0*v5_brick[mz][my][mx];
}
}
}
}
// electrostatic potential
phi[i] = u_pa;
// convert E-field to force
const double qfactor = 0.5*force->qqrd2e * scale * q[i];
if (eflag_atom) eatom[i] += u_pa*qfactor;
if (vflag_atom) {
vatom[i][0] += v0*qfactor;
vatom[i][1] += v1*qfactor;
vatom[i][2] += v2*qfactor;
vatom[i][3] += v3*qfactor;
vatom[i][4] += v4*qfactor;
vatom[i][5] += v5*qfactor;
}
}
}
/* ----------------------------------------------------------------------
Slab-geometry correction term to dampen inter-slab interactions between
periodically repeating slabs. Yields good approximation to 2D Ewald if
adequate empty space is left between repeating slabs (J. Chem. Phys.
111, 3155). Slabs defined here to be parallel to the xy plane. Also
extended to non-neutral systems (J. Chem. Phys. 131, 094107).
------------------------------------------------------------------------- */
void PPPMDispDielectric::slabcorr(int eflag)
{
// compute local contribution to global dipole moment
double *q = atom->q;
double **x = atom->x;
double *eps = atom->epsilon;
double zprd = domain->zprd;
int nlocal = atom->nlocal;
double dipole = 0.0;
for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2];
if (mu_flag) {
double **mu = atom->mu;
for (int i = 0; i < nlocal; i++) dipole += mu[i][2];
}
// sum local contributions to get global dipole moment
double dipole_all;
MPI_Allreduce(&dipole,&dipole_all,1,MPI_DOUBLE,MPI_SUM,world);
// need to make non-neutral systems and/or
// per-atom energy translationally invariant
double dipole_r2 = 0.0;
if (eflag_atom || fabs(qsum) > SMALL) {
if (mu_flag)
error->all(FLERR,"Cannot (yet) use kspace slab correction with "
"long-range dipoles and non-neutral systems or per-atom energy");
for (int i = 0; i < nlocal; i++)
dipole_r2 += q[i]*x[i][2]*x[i][2];
// sum local contributions
double tmp;
MPI_Allreduce(&dipole_r2,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
dipole_r2 = tmp;
}
// compute corrections
const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
const double qscale = qqrd2e * scale;
if (eflag_global) energy += qscale * e_slabcorr;
// per-atom energy
if (eflag_atom) {
double efact = qscale * MY_2PI/volume;
for (int i = 0; i < nlocal; i++)
eatom[i] += efact * eps[i]*q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 +
qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0);
}
// add on force corrections
double ffact = qscale * (-4.0*MY_PI/volume);
double **f = atom->f;
for (int i = 0; i < nlocal; i++) {
f[i][2] += ffact * eps[i]*q[i]*(dipole_all - qsum*x[i][2]);
efield[i][2] += ffact * eps[i]*(dipole_all - qsum*x[i][2]);
}
// add on torque corrections
if (mu_flag && atom->torque) {
double **mu = atom->mu;
double **torque = atom->torque;
for (int i = 0; i < nlocal; i++) {
torque[i][0] += ffact * dipole_all * mu[i][1];
torque[i][1] += -ffact * dipole_all * mu[i][0];
}
}
}
/* ----------------------------------------------------------------------
memory usage of local arrays
------------------------------------------------------------------------- */
double PPPMDispDielectric::memory_usage()
{
double bytes = PPPMDisp::memory_usage();
bytes += nmax*3 * sizeof(double);
bytes += nmax * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
compute qsum,qsqsum,q2 and give error/warning if not charge neutral
called initially, when particle count changes, when charges are changed
------------------------------------------------------------------------- */
void PPPMDispDielectric::qsum_qsq()
{
const double * const q = atom->q;
const int nlocal = atom->nlocal;
double qsum_local(0.0), qsqsum_local(0.0);
for (int i = 0; i < nlocal; i++) {
qsum_local += q[i];
qsqsum_local += q[i]*q[i];
}
MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
q2 = qsqsum * force->qqrd2e;
}

62
src/DIELECTRIC/pppm_disp_dielectric.h Normal file
View File

@ -0,0 +1,62 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef KSPACE_CLASS
// clang-format off
KSpaceStyle(pppm/disp/dielectric,PPPMDispDielectric);
// clang-format on
#else
#ifndef LMP_PPPM_DISP_DIELECTRIC_H
#define LMP_PPPM_DISP_DIELECTRIC_H
#include "pppm_disp.h"
namespace LAMMPS_NS {
class PPPMDispDielectric : public PPPMDisp {
public:
PPPMDispDielectric(class LAMMPS *);
virtual ~PPPMDispDielectric();
virtual double memory_usage();
virtual void compute(int, int);
void qsum_qsq();
void slabcorr(int);
double **efield;
double *phi;
int potflag; // 1/0 if per-atom electrostatic potential phi is needed
protected:
virtual void fieldforce_c_ik();
virtual void fieldforce_c_ad();
virtual void fieldforce_c_peratom();
class AtomVecDielectric *avec;
int mu_flag;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/