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

Merge pull request #3480 from aliehlen/dielectic-doc-updates

Browse Source 
Documentation updates for DIELECTRIC package
This commit is contained in:
Axel Kohlmeyer
2022-10-10 19:38:16 -04:00
committed by GitHub
parent 26f05b336f 1bc16d8a78
commit 18f91bc565
8 changed files with 183 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
doc/src/atom_style.rst
View File

@ -91,7 +91,7 @@ quantities.
+--------------+-----------------------------------------------------+--------------------------------------+
| *charge* | charge | atomic system with charges |
+--------------+-----------------------------------------------------+--------------------------------------+
| *dielectric* | dipole, area, curvature | system with surface polarization |
| *dielectric* | normx normy normz area/patch ed em epsilon curv | system with surface polarization |
+--------------+-----------------------------------------------------+--------------------------------------+
| *dipole* | charge and dipole moment | system with dipolar particles |
+--------------+-----------------------------------------------------+--------------------------------------+
@ -180,16 +180,21 @@ vector with the 3 diameters of the ellipsoid and a quaternion 4-vector
with its orientation.
For the *dielectric* style, each particle can be either a physical
particle (e.g. an ion), or an interface particle representing a
boundary element. For physical particles, the per-particle properties
are the same as atom_style full. For interface particles, in addition
to these properties, each particle also has an area, a normal unit
vector, a mean local curvature, the mean and difference of the
dielectric constants of two sides of the interface, and the local
dielectric constant at the boundary element. The distinction between
the physical and interface particles is only meaningful when :doc:`fix
polarize <fix_polarize>` commands are applied to the interface
particles.
particle (e.g. an ion), or an interface particle representing a boundary
element between two regions of different dielectric constant. For
interface particles, in addition to the properties associated with
atom_style full, each particle also should be assigned a normal unit
vector (defined by normx, normy, normz), an area (area/patch), the
difference and mean of the dielectric constants of two sides of the
interface along the direction of the normal vector (ed and em), the
local dielectric constant at the boundary element (epsilon), and a mean
local curvature (curv). Physical particles must be assigned these
values, as well, but only their local dielectric constants will be used;
see documentation for associated :doc:`pair styles <pair_dielectric>`
and :doc:`fixes <fix_polarize>`. The distinction between the physical
and interface particles is only meaningful when :doc:`fix polarize
<fix_polarize>` commands are applied to the interface particles. This
style is part of the DIELECTRIC package.
For the *dipole* style, a point dipole is defined for each point
particle. Note that if you wish the particles to be finite-size

183
doc/src/fix_polarize.rst
View File

@ -16,11 +16,11 @@ Syntax
.. parsed-literal::
fix ID group-ID style nevery tolerance ...
fix ID group-ID style nevery tolerance
* ID, group-ID are documented in :doc:`fix <fix>` command
* style = *polarize/bem/gmres* or *polarize/bem/icc* or *polarize/functional*
* Nevery = this fixed is invoked every this many timesteps
* nevery = this fixed is invoked every this many timesteps
* tolerance = the relative tolerance for the iterative solver to stop
@ -46,44 +46,53 @@ Description
These fixes compute induced charges at the interface between two
impermeable media with different dielectric constants. The interfaces
need to be discretized into vertices, each representing a boundary element.
The vertices are treated as if they were regular atoms or particles.
:doc:`atom_style dielectric <atom_style>` should be used since it defines
the additional properties of each interface particle such as
interface normal vectors, element areas, and local dielectric mismatch.
These fixes also require the use of :doc:`pair_style <pair_style>` and
:doc:`kspace_style <kspace_style>` with the *dielectric* suffix.
At every time step, given a configuration of the physical charges in the system
(such as atoms and charged particles) these fixes compute and update
the charge of the interface particles. The interfaces are allowed to move
during the simulation with appropriate time integrators (for example,
with :doc:`fix_rigid <fix_rigid>`).
need to be discretized into vertices, each representing a boundary
element. The vertices are treated as if they were regular atoms or
particles. :doc:`atom_style dielectric <atom_style>` should be used
since it defines the additional properties of each interface particle
such as interface normal vectors, element areas, and local dielectric
mismatch. These fixes also require the use of :doc:`pair_style
<pair_style>` and :doc:`kspace_style <kspace_style>` with the
*dielectric* suffix. At every time step, given a configuration of the
physical charges in the system (such as atoms and charged particles)
these fixes compute and update the charge of the interface
particles. The interfaces are allowed to move during the simulation if
the appropriate time integrators are also set (for example, with
:doc:`fix_rigid <fix_rigid>`).
Consider an interface between two media: one with dielectric constant
of 78 (water), the other of 4 (silica). The interface is discretized
into 2000 boundary elements, each represented by an interface particle. Suppose that
each interface particle has a normal unit vector pointing from the silica medium to water.
The dielectric difference along the normal vector is then 78 - 4 = 74,
the mean dielectric value is (78 + 4) / 2 = 41. Each boundary element
also has its area and the local mean curvature (which is used by these fixes
for computing a correction term in the local electric field).
To model charged interfaces, the interface particle will have a non-zero charge value,
Consider an interface between two media: one with dielectric constant of
78 (water), the other of 4 (silica). The interface is discretized into
2000 boundary elements, each represented by an interface
particle. Suppose that each interface particle has a normal unit vector
pointing from the silica medium to water. The dielectric difference
along the normal vector is then 78 - 4 = 74, the mean dielectric value
is (78 + 4) / 2 = 41. Each boundary element also has its area and the
local mean curvature, which is used by these fixes for computing a
correction term in the local electric field. To model charged
interfaces, the interface particle will have a non-zero charge value,
coming from its area and surface charge density.
For non-interface particles such as atoms and charged particles,
the interface normal vectors, element area, and dielectric mismatch are
irrelevant. Their local dielectric value is used to rescale their actual charge
when computing the Coulombic interactions. For instance, for a cation carrying
a charge of +2 (in charge unit) in an implicit solvent with dielectric constant of 40
would have actual charge of +2, and a local dielectric constant value of 40.
It is assumed that the particles cannot pass through the interface during the simulation
so that its local dielectric constant value does not change.
For non-interface particles such as atoms and charged particles, the
interface normal vectors, element area, and dielectric mismatch are
irrelevant and unused. Their local dielectric value is used internally
to rescale their given charge when computing the Coulombic
interactions. For instance, to simulate a cation carrying a charge of +2
(in simulation charge units) in an implicit solvent with a dielectric
constant of 40, the cation's charge should be set to +2 and its local
dielectric constant property (defined in the :doc:`atom_style dielectric
<atom_style>`) should be set to 40; there is no need to manually rescale
charge. This will produce the proper force for any :doc:`pair_style
<pair_style>` with the dielectric suffix. It is assumed that the
particles cannot pass through the interface during the simulation
because the value of the local dielectric constant property does not
change.
There are some example scripts for using these fixes
with LAMMPS in the ``examples/PACKAGES/dielectric`` directory. The README file
therein contains specific details on the system setup. Note that the example data files
show the additional fields (columns) needed for :doc:`atom_style dielectric <atom_style>`
beyond the conventional fields *id*, *mol*, *type*, *q*, *x*, *y*, and *z*.
There are some example scripts for using these fixes with LAMMPS in the
``examples/PACKAGES/dielectric`` directory. The README file therein
contains specific details on the system setup. Note that the example
data files show the additional fields (columns) needed for
:doc:`atom_style dielectric <atom_style>` beyond the conventional fields
*id*, *mol*, *type*, *q*, *x*, *y*, and *z*.
----------
@ -104,22 +113,24 @@ the interface, are computed using the equation:
* :math:`\mathbf{E}(\mathbf{s})` is the electrical field at the vertex
* :math:`\mathbf{n}(\mathbf{s})` is the unit normal vector at the vertex pointing from medium with :math:`\epsilon_2` to that with :math:`\epsilon_1`
Fix *polarize/bem/gmres* employs the Generalized Minimum Residual (GMRES)
as described in :ref:`(Barros) <Barros>` to solve :math:`\sigma_b`.
Fix *polarize/bem/gmres* employs the Generalized Minimum Residual
(GMRES) as described in :ref:`(Barros) <Barros>` to solve
:math:`\sigma_b`.
Fix *polarize/bem/icc* employs the successive over-relaxation algorithm
as described in :ref:`(Tyagi) <Tyagi>` to solve :math:`\sigma_b`.
The iterative solvers would terminate either when the maximum relative change
in the induced charges in consecutive iterations is below the set tolerance,
or when the number of iterations reaches *iter_max* (see below).
The iterative solvers would terminate either when the maximum relative
change in the induced charges in consecutive iterations is below the set
tolerance, or when the number of iterations reaches *iter_max* (see
below).
Fix *polarize/functional* employs the energy functional variation approach
as described in :ref:`(Jadhao) <Jadhao>` to solve :math:`\sigma_b`.
Fix *polarize/functional* employs the energy functional variation
approach as described in :ref:`(Jadhao) <Jadhao>` to solve
:math:`\sigma_b`.
More details on the implementation of these fixes and their recommended use
are described in :ref:`(NguyenTD) <NguyenTD>`.
More details on the implementation of these fixes and their recommended
use are described in :ref:`(NguyenTD) <NguyenTD>`.
Restart, fix_modify, output, run start/stop, minimize info
@ -127,35 +138,78 @@ Restart, fix_modify, output, run start/stop, minimize info
No information about this fix is written to :doc:`binary restart files <restart>`.
The :doc:`fix_modify <fix_modify>` command provides certain options to
control the induced charge solver and the initial values of the interface elements:
The :doc:`fix_modify <fix_modify>` command provides the ability to modify certain
settings:
.. parsed-literal::
*itr_max* arg
arg = maximum number of iterations for convergence
*dielectrics* ediff emean epsilon area charge
ediff = dielectric difference
emean = dielectric mean
epsilon = local dielectric value
aree = element area
charge = real interface charge
ediff = dielectric difference or NULL
emean = dielectric mean or NULL
epsilon = local dielectric value or NULL
area = element area or NULL
charge = real interface charge or NULL
*kspace* arg = yes or no
*rand* max seed
max = range of random induced charges to be generated
seed = random number seed to use when generating random charge
*mr* arg
arg = maximum number of q-vectors to use when solving (GMRES only)
*omega* arg
arg = relaxation parameter to use when iterating (ICC only)
*polarize/bem/gmres* or *polarize/bem/icc* compute a global 2-element vector
which can be accessed by various :doc:`output commands <Howto_output>`.
The first element is the number of iterations when the solver terminates
(of which the upper bound is set by *iter_max*). The second element is the RMS error.
The *itr_max* keyword sets the max number of iterations to be used for
solving each step.
The *dielectrics* keyword allows properties of the atoms in group
*group-ID* to be modified. Values passed to any of the arguments
(*ediff*, *emean*, *epsilon*, *area*, *charge*) will override existing
values for all atoms in the group *group-ID*. Passing NULL to any of
these arguments will preserve the existing value. Note that setting the
properties of the interface this way will change the properties of all
atoms associated with the fix (all atoms in *group-ID*), so multiple fix
and fix_modify commands would be needed to change the properties of two
different interfaces to different values (one fix and fix_modify for
each interface group).
The *kspace* keyword turns on long range interactions.
If the arguments of the *rand* keyword are set, then the atoms subject
to this fix will be assigned a random initial charge in a uniform
distribution from -*max*/2 to *max*/2, using random number seed *seed*.
The *mr* keyword only applies to *style* = *polarize/bem/gmres*. It is
the maximum number of q-vectors to use when solving for the surface
charge.
The *omega* keyword only applies when using *style* =
*polarize/bem/icc*. It is a relaxation parameter defined in
:ref:`(Tyagi) <Tyagi>` that should generally be set between 0 and 2.
Note that the local dielectric constant (epsilon) can also be set
independently using the :doc:`set <set>` command.
----------
*polarize/bem/gmres* or *polarize/bem/icc* compute a global 2-element
vector which can be accessed by various :doc:`output commands
<Howto_output>`. The first element is the number of iterations when the
solver terminates (of which the upper bound is set by *iter_max*). The
second element is the RMS error.
Restrictions
""""""""""""
These fixes are part of the DIELECTRIC package. It is only enabled
These fixes are part of the DIELECTRIC package. They are only enabled
if LAMMPS was built with that package, which requires that also the
KSPACE package is installed. See the :doc:`Build package
<Build_package>` page for more info.
Note that the *polarize/bem/gmres* and *polarize/bem/icc* fixes only support
:doc:`units <units>` *lj*, *real*, *metal*, *si* and *nano* at the moment.
Note that the *polarize/bem/gmres* and *polarize/bem/icc* fixes only
support :doc:`units <units>` *lj*, *real*, *metal*, *si* and *nano* at
the moment.
Related commands
@ -171,6 +225,15 @@ Default
*iter_max* = 20
*kspace* = yes
*omega* = 0.7 (ICC only)
*mr* = \# atoms in group *group-ID* minus 1 (GMRES only)
No random charge initialization happens by default.
----------
.. _Barros:

34
doc/src/pair_dielectric.rst
View File

@ -76,16 +76,19 @@ Description
"""""""""""
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.
without the *dielectric* suffix. In addition to computing atom forces
and energies, these pair styles compute the electric field vector at
each atom, which are intended to be used by the :doc:`fix polarize
<fix_polarize>` commands to compute induced charges at interfaces
between two regions of different dielectric constant.
These pair styles should be used with :doc:`atom_style dielectric <atom_style>`,
which uses atom charges rescaled by their local dielectric constant.
These pair styles should be used with :doc:`atom_style dielectric
<atom_style>`.
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.
lj/long/coul/long/dielectric should be used with their kspace style
counterparts, namely, pppm/dielectric, pppm/disp/dielectric, and
msm/dielectric, respectively.
----------
@ -97,24 +100,27 @@ Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
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.
and cutoff distances for this pair style can be mixed. The default mix
algorithm is *geometric*\ . See the :doc:`pair_modify <pair_modify>`"
command for details.
The :doc:`pair_modify <pair_modify>` table option is not relevant
for this pair style.
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.
These pair styles write 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.
This pair style can only be used via the *pair* keyword of the
These pair styles can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the
*inner*, *middle*, *outer* keywords.
Restrictions
""""""""""""
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>` page for more info.
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>` page for more info.
Related commands
""""""""""""""""

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

@ -600,6 +600,7 @@ Cummins
Cundall
cundall
Curk
curv
Cusentino
customIDs
cutbond

10
examples/PACKAGES/dielectric/README
View File

@ -2,7 +2,7 @@ This folder contains some example data and input scripts for the DIELECTRIC pack
Nguyen TD, Li H, Bagchi D, Solis FJ, Olvera de la Cruz, Incorporating surface polarization effects into large-scale coarse-grained molecular dynamics simulation, Computer Physics Communications 2019, 241, 80--91.
- data.confined : two point opposite charges confined between two interfaces (epsilon1=2/epsilon2=10/epsilon2=2)
- data.confined : two point opposite charges confined between two interfaces (epsilon2=2/epsilon1=10/epsilon2=2)
- data.sphere : two point opposite charges outside a spherical interface (epsilon_in=1/epsilon2=10)
- in.confined : read in data.confined
@ -10,7 +10,7 @@ Nguyen TD, Li H, Bagchi D, Solis FJ, Olvera de la Cruz, Incorporating surface po
For "atom_style dielectric" the Atoms section in the data file contains 15 following columns:
id mol type q x y z normx normy normz area_per_patch ed em epsilon curvature
id mol type q x y z normx normy normz area/patch ed em epsilon curvature
where
@ -34,9 +34,13 @@ where
For interface particles, epsilon is set to be em
(the mean dielectric value above).
* area_per_patch: the surface area of the patch (element).
* area/patch: the surface area of the patch (element).
For real charges, this value is irrelevant, can be 1.0.
* curvature: surface mean curvature at the patch.
For example, for spherical interfaces, curvature = 1/spherical radius.
For planar interfaces, curvature = 0.
Note that the properties normx, normy, normz, area/patch, ed, em, and curvature are not
used for the non-interface beads. epsilon is used to scale the charge of any non-interface
ion, see the documentation for pair styles with the dielectric suffix and fix polarize.

2
examples/PACKAGES/dielectric/in.confined
View File

@ -7,7 +7,7 @@
# Dielectric constants can be set to be different from the input data file
variable epsilon1 index 20
variable epsilon2 index 8
variable epsilon2 index 10
variable data index data.confined

23
src/DIELECTRIC/fix_polarize_bem_gmres.cpp
View File

@ -265,16 +265,19 @@ void FixPolarizeBEMGMRES::setup(int /*vflag*/)
else
error->all(FLERR, "Pair style not compatible with fix polarize/bem/gmres");
if (kspaceflag) {
if (force->kspace) {
if (strcmp(force->kspace_style, "pppm/dielectric") == 0)
efield_kspace = (dynamic_cast<PPPMDielectric *>(force->kspace))->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = (dynamic_cast<MSMDielectric *>(force->kspace))->efield;
else
error->all(FLERR, "Kspace style not compatible with fix polarize/bem/gmres");
} else
error->all(FLERR, "No Kspace style available for fix polarize/bem/gmres");
if (force->kspace) {
kspaceflag = 1;
if (strcmp(force->kspace_style, "pppm/dielectric") == 0)
efield_kspace = (dynamic_cast<PPPMDielectric *>(force->kspace))->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = (dynamic_cast<MSMDielectric *>(force->kspace))->efield;
else
error->all(FLERR, "Kspace style not compatible with fix polarize/bem/gmres");
} else {
if (kspaceflag == 1) { // users specified kspace yes but there is no kspace pair style
error->warning(FLERR, "No Kspace pair style available for fix polarize/bem/gmres");
kspaceflag = 0;
}
}
// NOTE: epsilon0e2q converts (epsilon0 * efield) to the unit of (charge unit / squared distance unit)

5
src/DIELECTRIC/fix_polarize_bem_icc.cpp
View File

@ -175,10 +175,9 @@ void FixPolarizeBEMICC::setup(int /*vflag*/)
efield_kspace = (dynamic_cast<MSMDielectric *>(force->kspace))->efield;
else
error->all(FLERR, "Kspace style not compatible with fix polarize/bem/icc");
} else {
if (kspaceflag == 1) { // users specified kspace yes
error->warning(FLERR, "No Kspace style available for fix polarize/bem/icc");
if (kspaceflag == 1) { // users specified kspace yes but there is no kspace pair style
error->warning(FLERR, "No Kspace pair style available for fix polarize/bem/icc");
kspaceflag = 0;
}
}