From b3fa2b7c4f97afad17f19683ba9d4a583af3cfcb Mon Sep 17 00:00:00 2001 From: ali Date: Fri, 5 Aug 2022 11:25:03 -0500 Subject: [PATCH] update of DIELECTRIC package documentation and examples --- doc/src/atom_style.rst | 27 +++++---- doc/src/fix_polarize.rst | 73 ++++++++++++++++++------ doc/src/pair_dielectric.rst | 10 +++- examples/PACKAGES/dielectric/README | 6 +- examples/PACKAGES/dielectric/in.confined | 2 +- 5 files changed, 84 insertions(+), 34 deletions(-) diff --git a/doc/src/atom_style.rst b/doc/src/atom_style.rst index b753fdbe5f..dd650c6384 100644 --- a/doc/src/atom_style.rst +++ b/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 -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 ` 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_per_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 ` and :doc:`fixes `. +The distinction between the physical and interface +particles is only meaningful when :doc:`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 diff --git a/doc/src/fix_polarize.rst b/doc/src/fix_polarize.rst index 6ed3b36c55..5217563e89 100644 --- a/doc/src/fix_polarize.rst +++ b/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 ` 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 ` 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 `). 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 -for computing a correction term in the local electric field). +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. +irrelevant. Their local dielectric value is used 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/40 = +0.05 and +its local dielectric constant property (defined in the +:doc:`atom_style dielectric `) should be set to 40. This will produce +the proper force for any :doc:`pair_style ` with the dielectric suffix. +As noted in :doc:`pair_style `, 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 `. -The :doc:`fix_modify ` command provides certain options to -control the induced charge solver and the initial values of the interface elements: +The :doc:`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) + +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) `. + +The *omega* keyword only applies when using *style* = *polarize/bem/icc*. It is a relaxation parameter defined in :ref:`(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 `. @@ -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: diff --git a/doc/src/pair_dielectric.rst b/doc/src/pair_dielectric.rst index dbdd163a76..ef1f186021 100644 --- a/doc/src/pair_dielectric.rst +++ b/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 ` commands. +at each atom, which are intended to be used by the :doc:`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 `, -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 ` 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, diff --git a/examples/PACKAGES/dielectric/README b/examples/PACKAGES/dielectric/README index 1d5f30263c..2b43c5c299 100644 --- a/examples/PACKAGES/dielectric/README +++ b/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. diff --git a/examples/PACKAGES/dielectric/in.confined b/examples/PACKAGES/dielectric/in.confined index beb5b9a2b0..f487cbc1f0 100644 --- a/examples/PACKAGES/dielectric/in.confined +++ b/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