diff --git a/doc/src/atom_style.rst b/doc/src/atom_style.rst index b753fdbe5f..c3b16c6ad1 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/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/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..6a997aa7ea 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 @@ -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 ` 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 ` and -:doc:`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 `). +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 ` 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 +` and :doc:`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 `). -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 +`) should be set to 40; there is no need to manually rescale +charge. This will produce the proper force for any :doc:`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 ` -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 ` 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) ` to solve :math:`\sigma_b`. +Fix *polarize/bem/gmres* employs the Generalized Minimum Residual +(GMRES) as described in :ref:`(Barros) ` to solve +:math:`\sigma_b`. Fix *polarize/bem/icc* employs the successive over-relaxation algorithm as described in :ref:`(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) ` to solve :math:`\sigma_b`. +Fix *polarize/functional* employs the energy functional variation +approach as described in :ref:`(Jadhao) ` to solve +:math:`\sigma_b`. - -More details on the implementation of these fixes and their recommended use -are described in :ref:`(NguyenTD) `. +More details on the implementation of these fixes and their recommended +use are described in :ref:`(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 `. -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) -*polarize/bem/gmres* or *polarize/bem/icc* compute a global 2-element vector -which can be accessed by various :doc:`output commands `. -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) ` 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 ` command. + +---------- + +*polarize/bem/gmres* or *polarize/bem/icc* compute a global 2-element +vector which can be accessed by various :doc:`output commands +`. 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 ` page for more info. -Note that the *polarize/bem/gmres* and *polarize/bem/icc* fixes only support -:doc:`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 ` *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: diff --git a/doc/src/pair_dielectric.rst b/doc/src/pair_dielectric.rst index dbdd163a76..0b60d202ea 100644 --- a/doc/src/pair_dielectric.rst +++ b/doc/src/pair_dielectric.rst @@ -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 ` 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 +` 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. +These pair styles should be used with :doc:`atom_style dielectric +`. 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 `" +command for details. The :doc:`pair_modify ` table option is not relevant for this pair style. -This pair style writes its information to :doc:`binary restart files `, 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 +`, 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 ` 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 ` 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 +` page for more info. Related commands """""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 1e6add76fc..df93ebf761 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -600,6 +600,7 @@ Cummins Cundall cundall Curk +curv Cusentino customIDs cutbond diff --git a/examples/PACKAGES/dielectric/README b/examples/PACKAGES/dielectric/README index 1d5f30263c..ef775e479f 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 @@ -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. diff --git a/examples/PACKAGES/dielectric/in.confined b/examples/PACKAGES/dielectric/in.confined index beb5b9a2b0..0f9dab7bba 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 variable data index data.confined diff --git a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp b/src/DIELECTRIC/fix_polarize_bem_gmres.cpp index 7d29cc889c..07303addea 100644 --- a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp +++ b/src/DIELECTRIC/fix_polarize_bem_gmres.cpp @@ -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(force->kspace))->efield; - else if (strcmp(force->kspace_style, "msm/dielectric") == 0) - efield_kspace = (dynamic_cast(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(force->kspace))->efield; + else if (strcmp(force->kspace_style, "msm/dielectric") == 0) + efield_kspace = (dynamic_cast(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) diff --git a/src/DIELECTRIC/fix_polarize_bem_icc.cpp b/src/DIELECTRIC/fix_polarize_bem_icc.cpp index 3d451256ff..9b778122af 100644 --- a/src/DIELECTRIC/fix_polarize_bem_icc.cpp +++ b/src/DIELECTRIC/fix_polarize_bem_icc.cpp @@ -175,10 +175,9 @@ void FixPolarizeBEMICC::setup(int /*vflag*/) efield_kspace = (dynamic_cast(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; } }