update of DIELECTRIC package documentation and examples

doc/src/atom_style.rst

@@ -91,7 +91,7 @@ quantities.
 +--------------+-----------------------------------------------------+--------------------------------------+
 | *charge*     | charge                                              | atomic system with charges           |
 +--------------+-----------------------------------------------------+--------------------------------------+
-| *dielectric* | dipole, area, curvature                             | system with surface polarization     |
+| *dielectric* | normx normy normz area_per_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
+particle (e.g. an ion), or an interface particle representing a boundary
-boundary element. For physical particles, the per-particle properties
+element between two regions of different dielectric constant. For 
-are the same as atom_style full.  For interface particles, in addition
+interface particles, in addition to the properties associated with
-to these properties, each particle also has an area, a normal unit
+atom_style full, each particle also should be assigned a normal unit vector 
-vector, a mean local curvature, the mean and difference of the
+(defined by normx, normy, normz), an area (area_per_patch), the difference and 
-dielectric constants of two sides of the interface, and the local
+mean of the dielectric constants of two sides of the interface along the 
-dielectric constant at the boundary element.  The distinction between
+direction of the normal vector (ed and em), the local dielectric constant at the
-the physical and interface particles is only meaningful when :doc:`fix
+boundary element (epsilon), and a mean local curvature (curv). 
-polarize <fix_polarize>` commands are applied to the interface
+Physical particles must be assigned these values, as well, 
-particles.
+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

doc/src/fix_polarize.rst

@@ -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
 
 
@@ -56,7 +56,7 @@ These fixes also require the use of :doc:`pair_style <pair_style>` and
 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,
+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
@@ -65,19 +65,25 @@ into 2000 boundary elements, each represented by an interface particle. Suppose
 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
+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).
+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
+irrelevant. Their local dielectric value is used to rescale their given charge
-when computing the Coulombic interactions. For instance, for a cation carrying
+when computing the Coulombic interactions. For instance, to simulate a cation 
-a charge of +2 (in charge unit) in an implicit solvent with dielectric constant of 40
+carrying a charge of +2 (in simulation charge units) in an implicit solvent with 
-would have actual charge of +2, and a local dielectric constant value of 40.
+a dielectric constant of 40, the cation's charge should be set to +2/40 = +0.05 and 
-It is assumed that the particles cannot pass through the interface during the simulation
+its local dielectric constant property (defined in the 
-so that its local dielectric constant value does not change.
+:doc:`atom_style dielectric <atom_style>`) should be set to 40. This will produce 
+the proper force for any :doc:`pair_style <pair_style>` with the dielectric suffix. 
+As noted in :doc:`pair_style <pair_dielectric>`, the charge (+0.05) will be scaled by the 
+local dielectric constant (40) to produce the appropriate charge (+2) during 
+force computation. 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
@@ -127,18 +133,40 @@ 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
+The :doc:`fix_modify <fix_modify>` command provides the ability to modify certain 
-control the induced charge solver and the initial values of the interface elements:
+settings:
 
   .. parsed-literal::
       *itr_max* arg
          arg = maximum number of iterations for convergence
       *dielectrics* ediff emean epsilon area charge
-         ediff = dielectric difference
+         ediff = dielectric difference or NULL
-         emean = dielectric mean
+         emean = dielectric mean or NULL
-         epsilon = local dielectric value
+         epsilon = local dielectric value or NULL
-         aree = element area
+         area = element area or NULL
-         charge = real interface charge
+         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)
+
+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.
+
+The *kspace* keyword turns on long range interactions.
+
+If the argumnts 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, see :ref:`(Barros) <Barros>`.
+
+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.
+
+----------
 
 *polarize/bem/gmres* or *polarize/bem/icc* compute a global 2-element vector
 which can be accessed by various :doc:`output commands <Howto_output>`.
@@ -171,6 +199,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:

doc/src/pair_dielectric.rst

@@ -76,12 +76,16 @@ Description
 """""""""""
 
 All these pair styles are derived from the corresponding pair styles
-without the *dielectric*\ suffix. In addition to computing atom forces
+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.
+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.
+which uses atom charges rescaled by their local dielectric constant. See the note in 
+:doc:`fix polarize <fix_polarize>` about properly setting ion charge given a 
+local dielectric constant.
 
 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,

examples/PACKAGES/dielectric/README

@@ -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
@@ -40,3 +40,7 @@ where
 * 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_per_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.

examples/PACKAGES/dielectric/in.confined

@@ -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 # must be the same as the local epsilon defined in data.confined
 
 variable    data index data.confined