diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt
index 5e3247af994626ac3392a424e487c56647ce75f4..63a63b373693198c9950adcf1c1bbc05cc329177 100644
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@@ -940,6 +940,8 @@ KOKKOS, o = USER-OMP, t = OPT.
"lj/charmm/coul/charmm/implicit (ko)"_pair_charmm.html,
"lj/charmm/coul/long (giko)"_pair_charmm.html,
"lj/charmm/coul/msm"_pair_charmm.html,
+"lj/charmmfsw/coul/charmmfsh"_pair_charmm.html,
+"lj/charmmfsw/coul/long"_pair_charmm.html,
"lj/class2 (gko)"_pair_class2.html,
"lj/class2/coul/cut (ko)"_pair_class2.html,
"lj/class2/coul/long (gko)"_pair_class2.html,
@@ -1148,6 +1150,7 @@ USER-OMP, t = OPT.
"zero"_dihedral_zero.html,
"hybrid"_dihedral_hybrid.html,
"charmm (ko)"_dihedral_charmm.html,
+"charmmfsh"_dihedral_charmm.html,
"class2 (ko)"_dihedral_class2.html,
"harmonic (io)"_dihedral_harmonic.html,
"helix (o)"_dihedral_helix.html,
diff --git a/doc/src/Section_example.txt b/doc/src/Section_example.txt
index 0fadaf11fd430f5111e5096c64605c5dc00fa69e..26dc3b96985afbf2c3fe7956d7f395b66be07229 100644
--- a/doc/src/Section_example.txt
+++ b/doc/src/Section_example.txt
@@ -51,9 +51,11 @@ Lowercase directories :h4
accelerate: run with various acceleration options (OpenMP, GPU, Phi)
balance: dynamic load balancing, 2d system
body: body particles, 2d system
+cmap: CMAP 5-body contributions to CHARMM force field
colloid: big colloid particles in a small particle solvent, 2d system
comb: models using the COMB potential
coreshell: core/shell model using CORESHELL package
+controller: use of fix controller as a thermostat
crack: crack propagation in a 2d solid
deposit: deposit atoms and molecules on a surface
dipole: point dipolar particles, 2d system
@@ -62,6 +64,8 @@ eim: NaCl using the EIM potential
ellipse: ellipsoidal particles in spherical solvent, 2d system
flow: Couette and Poiseuille flow in a 2d channel
friction: frictional contact of spherical asperities between 2d surfaces
+gcmc: Grand Canonical Monte Carlo (GCMC) via the fix gcmc command
+granregion: use of fix wall/region/gran as boundary on granular particles
hugoniostat: Hugoniostat shock dynamics
indent: spherical indenter into a 2d solid
kim: use of potentials in Knowledge Base for Interatomic Models (KIM)
@@ -69,6 +73,7 @@ meam: MEAM test for SiC and shear (same as shear examples)
melt: rapid melt of 3d LJ system
micelle: self-assembly of small lipid-like molecules into 2d bilayers
min: energy minimization of 2d LJ melt
+mscg: parameterize a multi-scale coarse-graining (MSCG) model
msst: MSST shock dynamics
nb3b: use of nonbonded 3-body harmonic pair style
neb: nudged elastic band (NEB) calculation for barrier finding
@@ -87,7 +92,8 @@ snap: NVE dynamics for BCC tantalum crystal using SNAP potential
srd: stochastic rotation dynamics (SRD) particles as solvent
streitz: use of Streitz/Mintmire potential with charge equilibration
tad: temperature-accelerated dynamics of vacancy diffusion in bulk Si
-vashishta: use of the Vashishta potential :tb(s=:)
+vashishta: use of the Vashishta potential
+voronoi: Voronoi tesselation via compute voronoi/atom command :tb(s=:)
Here is how you can run and visualize one of the sample problems:
diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt
index e83f5399c0a116c086fd3df5bc09414c7398d4b0..a378d8d088c9f708c4d2d944d5e8ef2cbf16a995 100644
--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@@ -204,7 +204,10 @@ documentation for the formula it computes.
"bond_style"_bond_harmonic.html harmonic
"angle_style"_angle_charmm.html charmm
+"dihedral_style"_dihedral_charmm.html charmmfsh
"dihedral_style"_dihedral_charmm.html charmm
+"pair_style"_pair_charmm.html lj/charmmfsw/coul/charmmfsh
+"pair_style"_pair_charmm.html lj/charmmfsw/coul/long
"pair_style"_pair_charmm.html lj/charmm/coul/charmm
"pair_style"_pair_charmm.html lj/charmm/coul/charmm/implicit
"pair_style"_pair_charmm.html lj/charmm/coul/long :ul
@@ -212,6 +215,12 @@ documentation for the formula it computes.
"special_bonds"_special_bonds.html charmm
"special_bonds"_special_bonds.html amber :ul
+NOTE: For CHARMM, the newer {charmmfsw} or {charmmfsh} styles were
+released in March 2017. We recommend they be used instead of the
+older {charmm} styles. See discussion of the differences on the "pair
+charmm"_pair_charmm.html and "dihedral charmm"_dihedral_charmm.html
+doc pages.
+
DREIDING is a generic force field developed by the "Goddard
group"_http://www.wag.caltech.edu at Caltech and is useful for
predicting structures and dynamics of organic, biological and
diff --git a/doc/src/compute_modify.txt b/doc/src/compute_modify.txt
index 637f9b5e417da481cb173670cf3eb21c2a7fb929..c26b5c2a5fe20273230b2c019982012efd048e48 100644
--- a/doc/src/compute_modify.txt
+++ b/doc/src/compute_modify.txt
@@ -19,7 +19,7 @@ keyword = {extra/dof} or {extra} or {dynamic/dof} or {dynamic} :l
N = # of extra degrees of freedom to subtract
{extra} syntax is identical to {extra/dof}, will be disabled at some point
{dynamic/dof} value = {yes} or {no}
- yes/no = do or do not recompute the number of atoms contributing to the temperature
+ yes/no = do or do not recompute the number of degrees of freedom (DOF) contributing to the temperature
{dynamic} syntax is identical to {dynamic/dof}, will be disabled at some point :pre
:ule
@@ -46,13 +46,16 @@ degrees-of-freedom. See the "compute
temp/asphere"_compute_temp_asphere.html command for an example.
The {dynamic/dof} or {dynamic} keyword determines whether the number
-of atoms N in the compute group is re-computed each time a temperature
-is computed. Only compute styles that calculate a temperature use
-this option. By default, N is assumed to be constant. If you are
-adding atoms to the system (see the "fix pour"_fix_pour.html, "fix
-deposit"_fix_deposit.html and "fix gcmc"_fix_gcmc.html commands) or
-expect atoms to be lost (e.g. due to evaporation), then this option
-should be used to insure the temperature is correctly normalized.
+of atoms N in the compute group and their associated degrees of
+freedom are re-computed each time a temperature is computed. Only
+compute styles that calculate a temperature use this option. By
+default, N and their DOF are assumed to be constant. If you are
+adding atoms or molecules to the system (see the "fix
+pour"_fix_pour.html, "fix deposit"_fix_deposit.html, and "fix
+gcmc"_fix_gcmc.html commands) or expect atoms or molecules to be lost
+(e.g. due to exiting the simulation box or via "fix
+evaporation"_fix_evaporation.html), then this option should be used to
+insure the temperature is correctly normalized.
NOTE: The {extra} and {dynamic} keywords should not be used as they
are deprecated (March 2017) and will eventually be disabled. Instead,
diff --git a/doc/src/dihedral_charmm.txt b/doc/src/dihedral_charmm.txt
index 87322cb0af29ad9b55b668821f4f2510c994806c..802cf2c4adc28bac4a9e53e2cd9c46be8d3044b4 100644
--- a/doc/src/dihedral_charmm.txt
+++ b/doc/src/dihedral_charmm.txt
@@ -10,21 +10,25 @@ dihedral_style charmm command :h3
dihedral_style charmm/intel command :h3
dihedral_style charmm/kk command :h3
dihedral_style charmm/omp command :h3
+dihedral_style charmmfsh command :h3
[Syntax:]
-dihedral_style charmm :pre
+dihedral_style style :pre
+
+style = {charmm} or {charmmfsh} :ul
[Examples:]
dihedral_style charmm
+dihedral_style charmmfsh
dihedral_coeff 1 0.2 1 180 1.0
dihedral_coeff 2 1.8 1 0 1.0
dihedral_coeff 1 3.1 2 180 0.5 :pre
[Description:]
-The {charmm} dihedral style uses the potential
+The {charmm} and {charmmfsh} dihedral styles use the potential
:c,image(Eqs/dihedral_charmm.jpg)
@@ -34,6 +38,11 @@ field (see comment on weighting factors below). See
"(Cornell)"_#dihedral-Cornell for a description of the AMBER force
field.
+NOTE: The newer {charmmfsh} style was released in March 2017. We
+recommend it be used instead of the older {charmm} style when running
+a simulation with the CHARMM force field. See the discussion below
+and more details on the "pair_style charmm"_pair_charmm.html doc page.
+
The following coefficients must be defined for each dihedral type via the
"dihedral_coeff"_dihedral_coeff.html command as in the example above, or in
the data file or restart files read by the "read_data"_read_data.html
@@ -73,13 +82,23 @@ special_bonds 1-4 scaling factor to 0.0 (which is the
default). Otherwise 1-4 non-bonded interactions in dihedrals will be
computed twice.
-Also note that for AMBER force fields, which use pair styles with
-"lj/cut", the special_bonds 1-4 scaling factor should be set to the
-AMBER defaults (1/2 and 5/6) and all the dihedral weighting factors
-(4th coeff above) must be set to 0.0. In this case, you can use any
-pair style you wish, since the dihedral does not need any
-Lennard-Jones parameter information and will not compute any 1-4
-non-bonded interactions.
+For simulations using the CHARMM force field, the difference between
+the {charmm} and {charmmfsh} styles is in the computation of the 1-4
+non-bond interactions, if the distance between the two atoms is within
+the switching distance of the pairwise potential defined by the
+corresponding CHARMM pair style, i.e. between the inner and outer
+cutoffs specified for the pair style. See the discussion on the
+"CHARMM pair_style"_pair_charmm.html doc page for details.
+
+Note that for AMBER force fields, which use pair styles with "lj/cut",
+the special_bonds 1-4 scaling factor should be set to the AMBER
+defaults (1/2 and 5/6) and all the dihedral weighting factors (4th
+coeff above) must be set to 0.0. In this case, you can use any pair
+style you wish, since the dihedral does not need any Lennard-Jones
+parameter information and will not compute any 1-4 non-bonded
+interactions. Likewise the {charmm} or {charmmfsh} styles are
+identical in this case since no 1-4 non-bonded interactions are
+computed.
:line
diff --git a/doc/src/kspace_style.txt b/doc/src/kspace_style.txt
index 5c4909f83ab31ed5343a0717e4fc850c1177570a..371540bd68b6e5c95c2473c66dc39c8aae9826da 100644
--- a/doc/src/kspace_style.txt
+++ b/doc/src/kspace_style.txt
@@ -182,13 +182,13 @@ currently support the -DFFT_SINGLE compiler switch.
:line
The {msm} style invokes a multi-level summation method MSM solver,
-"(Hardy)"_#Hardy2006 or "(Hardy2)"_#Hardy2009, which maps atom charge to a 3d
-mesh, and uses a multi-level hierarchy of coarser and coarser meshes
-on which direct coulomb solves are done. This method does not use
-FFTs and scales as N. It may therefore be faster than the other
+"(Hardy)"_#Hardy2006 or "(Hardy2)"_#Hardy2009, which maps atom charge
+to a 3d mesh, and uses a multi-level hierarchy of coarser and coarser
+meshes on which direct coulomb solves are done. This method does not
+use FFTs and scales as N. It may therefore be faster than the other
K-space solvers for relatively large problems when running on large
-core counts. MSM can also be used for non-periodic boundary conditions and
-for mixed periodic and non-periodic boundaries.
+core counts. MSM can also be used for non-periodic boundary conditions
+and for mixed periodic and non-periodic boundaries.
MSM is most competitive versus Ewald and PPPM when only relatively
low accuracy forces, about 1e-4 relative error or less accurate,
diff --git a/doc/src/pair_charmm.txt b/doc/src/pair_charmm.txt
index 3a1daf6c0340e66bb268faeeb219e109eab99766..dbaff7d13dc5c3dcc2ef83302446d2c4362cd08b 100644
--- a/doc/src/pair_charmm.txt
+++ b/doc/src/pair_charmm.txt
@@ -17,12 +17,14 @@ pair_style lj/charmm/coul/long/opt command :h3
pair_style lj/charmm/coul/long/omp command :h3
pair_style lj/charmm/coul/msm command :h3
pair_style lj/charmm/coul/msm/omp command :h3
+pair_style lj/charmmfsw/coul/charmmfsh command :h3
+pair_style lj/charmmfsw/coul/long command :h3
[Syntax:]
pair_style style args :pre
-style = {lj/charmm/coul/charmm} or {lj/charmm/coul/charmm/implicit} or {lj/charmm/coul/long} or {lj/charmm/coul/msm}
+style = {lj/charmm/coul/charmm} or {lj/charmm/coul/charmm/implicit} or {lj/charmm/coul/long} or {lj/charmm/coul/msm} or {lj/charmmfsw/coul/charmmfsh} or {lj/charmmfsw/coul/long}
args = list of arguments for a particular style :ul
{lj/charmm/coul/charmm} args = inner outer (inner2) (outer2)
inner, outer = global switching cutoffs for Lennard Jones (and Coulombic if only 2 args)
@@ -35,12 +37,20 @@ args = list of arguments for a particular style :ul
cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args)
{lj/charmm/coul/msm} args = inner outer (cutoff)
inner, outer = global switching cutoffs for LJ (and Coulombic if only 2 args)
+ cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args)
+ {lj/charmmfsw/coul/charmmfsh} args = inner outer (cutoff)
+ inner, outer = global cutoffs for LJ (and Coulombic if only 2 args)
+ cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args)
+ {lj/charmmfsw/coul/long} args = inner outer (cutoff)
+ inner, outer = global cutoffs for LJ (and Coulombic if only 2 args)
cutoff = global cutoff for Coulombic (optional, outer is Coulombic cutoff if only 2 args) :pre
[Examples:]
pair_style lj/charmm/coul/charmm 8.0 10.0
pair_style lj/charmm/coul/charmm 8.0 10.0 7.0 9.0
+pair_style lj/charmmfsw/coul/charmmfsh 8.0 10.0
+pair_style lj/charmmfsw/coul/charmmfsh 8.0 10.0 7.0 9.0
pair_coeff * * 100.0 2.0
pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
@@ -51,6 +61,8 @@ pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
pair_style lj/charmm/coul/long 8.0 10.0
pair_style lj/charmm/coul/long 8.0 10.0 9.0
+pair_style lj/charmmfsw/coul/long 8.0 10.0
+pair_style lj/charmmfsw/coul/long 8.0 10.0 9.0
pair_coeff * * 100.0 2.0
pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
@@ -61,21 +73,55 @@ pair_coeff 1 1 100.0 2.0 150.0 3.5 :pre
[Description:]
-The {lj/charmm} styles compute LJ and Coulombic interactions with an
-additional switching function S(r) that ramps the energy and force
-smoothly to zero between an inner and outer cutoff. It is a widely
-used potential in the "CHARMM"_http://www.scripps.edu/brooks MD code.
-See "(MacKerell)"_#pair-MacKerell for a description of the CHARMM force
-field.
+These pair styles compute Lennard Jones (LJ) and Coulombic
+interactions with additional switching or shifting functions that ramp
+the energy and/or force smoothly to zero between an inner and outer
+cutoff. They are implementations of the widely used CHARMM force
+field used in the "CHARMM"_http://www.scripps.edu/brooks MD code (and
+others). See "(MacKerell)"_#pair-MacKerell for a description of the
+CHARMM force field.
+
+The styles with {charmm} (not {charmmfsw} or {charmmfsh}) in their
+name are the older, original LAMMPS implemetations. They compute the
+LJ and Coulombic interactions with an energy switching function (esw,
+a cubic polynomial), which ramps the energy smoothly to zero between
+the inner and outer cutoff. This can cause irregularities in
+pair-wise forces (due to the discontinuous 2nd derivative of energy at
+the boundaries of the switching region), which in some cases can
+result in substantial artifacts in an MD simulation.
+
+The newer styles with {charmmfsw} or {charmmfsh} in their name replace
+the energy switching with force switching (fsw) and force shifting
+(fsh) functions, for LJ and Coulombic interactions respectively.
+These follow the formulas and description given in
+"(Steinbach)"_#Steinbach and "(Brooks)"_#Brooks to minimize these
+artifacts.
+
+NOTE: The newer {charmmfsw} or {charmmfsh} styles were released in
+March 2017. We recommend they be used instead of the older {charmm}
+styles. Eventually code from the new styles will propagate into the
+related pair styles (e.g. implicit, accelerator, free energy
+variants).
+
+The general CHARMM formulas are as follows
:c,image(Eqs/pair_charmm.jpg)
-Both the LJ and Coulombic terms require an inner and outer cutoff.
-They can be the same for both formulas or different depending on
-whether 2 or 4 arguments are used in the pair_style command. In each
-case, the inner cutoff distance must be less than the outer cutoff.
-It it typical to make the difference between the 2 cutoffs about 1.0
-Angstrom.
+where the specific form of the switching or shifting function S(r) is
+given in the "(Steinbach)"_#Steinbach paper.
+
+When using the {lj/charmm/coul/charmm styles}, both the LJ and
+Coulombic terms require an inner and outer cutoff. They can be the
+same for both formulas or different depending on whether 2 or 4
+arguments are used in the pair_style command. For the
+{lj/charmmfsw/coul/charmmfsh} style, the LJ term requires both an
+inner and outer cutoff, while the Coulombic term requires only one
+cutoff. If the Coulomb cutoff is not specified (2 instead of 3
+arguments), the LJ outer cutoff is used for the Coulombic cutoff. In
+all cases where an inner and outer cutoff are specified, the inner
+cutoff distance must be less than the outer cutoff. It is typical to
+make the difference between the inner and outer cutoffs about 2.0
+Angstroms.
Style {lj/charmm/coul/charmm/implicit} computes the same formulas as
style {lj/charmm/coul/charmm} except that an additional 1/r term is
@@ -86,12 +132,18 @@ screening. It is designed for use in a simulation of an unsolvated
biomolecule (no explicit water molecules).
Styles {lj/charmm/coul/long} and {lj/charmm/coul/msm} compute the same
-formulas as style {lj/charmm/coul/charmm} except that an additional
-damping factor is applied to the Coulombic term, as described for the
-"lj/cut"_pair_lj.html pair styles. Only one Coulombic cutoff is
-specified for {lj/charmm/coul/long} and {lj/charmm/coul/msm}; if only
-2 arguments are used in the pair_style command, then the outer LJ
-cutoff is used as the single Coulombic cutoff.
+formulas as style {lj/charmm/coul/charmm} and style
+{lj/charmmfsw/coul/long} computes the same formulas as style
+{lj/charmmfsw/coul/charmmfsh}, except that an additional damping
+factor is applied to the Coulombic term, so it can be used in
+conjunction with the "kspace_style"_kspace_style.html command and its
+{ewald} or {pppm} or {msm} option. Only one Coulombic cutoff is
+specified for these styles; if only 2 arguments are used in the
+pair_style command, then the outer LJ cutoff is used as the single
+Coulombic cutoff. The Coulombic cutoff specified for these styles
+means that pairwise interactions within this distance are computed
+directly; interactions outside that distance are computed in
+reciprocal space.
The following coefficients must be defined for each pair of atoms
types via the "pair_coeff"_pair_coeff.html command as in the examples
@@ -111,7 +163,7 @@ sigma.
The latter 2 coefficients are optional. If they are specified, they
are used in the LJ formula between 2 atoms of these types which are
also first and fourth atoms in any dihedral. No cutoffs are specified
-because this CHARMM force field does not allow varying cutoffs for
+because the CHARMM force field does not allow varying cutoffs for
individual atom pairs; all pairs use the global cutoff(s) specified in
the pair_style command.
@@ -148,38 +200,41 @@ mixed. The default mix value is {arithmetic} to coincide with the
usual settings for the CHARMM force field. See the "pair_modify"
command for details.
-None of the lj/charmm pair styles support the
+None of the {lj/charmm} or {lj/charmmfsw} pair styles support the
"pair_modify"_pair_modify.html shift option, since the Lennard-Jones
portion of the pair interaction is smoothed to 0.0 at the cutoff.
-The {lj/charmm/coul/long} style supports the
-"pair_modify"_pair_modify.html table option since it can tabulate the
-short-range portion of the long-range Coulombic interaction.
+The {lj/charmm/coul/long} and {lj/charmmfsw/coul/long} styles support
+the "pair_modify"_pair_modify.html table option since they can
+tabulate the short-range portion of the long-range Coulombic
+interaction.
-None of the lj/charmm pair styles support the
+None of the {lj/charmm} or {lj/charmmfsw} pair styles support the
"pair_modify"_pair_modify.html tail option for adding long-range tail
corrections to energy and pressure, since the Lennard-Jones portion of
the pair interaction is smoothed to 0.0 at the cutoff.
-All of the lj/charmm pair styles write their information to "binary
-restart files"_restart.html, so pair_style and pair_coeff commands do
-not need to be specified in an input script that reads a restart file.
+All of the {lj/charmm} and {lj/charmmfsw} pair styles write their
+information to "binary restart files"_restart.html, so pair_style and
+pair_coeff commands do not need to be specified in an input script
+that reads a restart file.
-The lj/charmm/coul/long pair style supports the use of the {inner},
-{middle}, and {outer} keywords of the "run_style respa"_run_style.html
-command, meaning the pairwise forces can be partitioned by distance at
-different levels of the rRESPA hierarchy. The other styles only
-support the {pair} keyword of run_style respa. See the
-"run_style"_run_style.html command for details.
+The {lj/charmm/coul/long} and {lj/charmmfsw/coul/long} pair styles
+support the use of the {inner}, {middle}, and {outer} keywords of the
+"run_style respa"_run_style.html command, meaning the pairwise forces
+can be partitioned by distance at different levels of the rRESPA
+hierarchy. The other styles only support the {pair} keyword of
+run_style respa. See the "run_style"_run_style.html command for
+details.
:line
[Restrictions:]
-The {lj/charmm/coul/charmm} and {lj/charmm/coul/charmm/implicit}
-styles are part of the MOLECULE package. The {lj/charmm/coul/long}
-style is part of the KSPACE package. They are only enabled if LAMMPS
-was built with those packages. See the "Making
+All the styles with {coul/charmm} or {coul/charmmfsh} styles are part
+of the MOLECULE package. All the styles with {coul/long} style are
+part of the KSPACE package. They are only enabled if LAMMPS was built
+with those packages. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info. Note that
the MOLECULE and KSPACE packages are installed by default.
@@ -191,6 +246,13 @@ the MOLECULE and KSPACE packages are installed by default.
:line
+:link(Brooks)
+[(Brooks)] Brooks, et al, J Comput Chem, 30, 1545 (2009).
+
:link(pair-MacKerell)
[(MacKerell)] MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
+
+:link(Stenibach)
+[(Steinbach)] Steinbach, Brooks, J Comput Chem, 15, 667 (1994).
+
diff --git a/doc/src/pair_momb.txt b/doc/src/pair_momb.txt
index ff822242a2135eb728cb12b2dbbd7fedc0496228..77c4f184d9380f11848bb7e0deedc719802d130d 100644
--- a/doc/src/pair_momb.txt
+++ b/doc/src/pair_momb.txt
@@ -25,18 +25,19 @@ pair_coeff 1 2 momb 0.0 1.0 1.0 10.2847 2.361 :pre
[Description:]
-Style {momb} computes pairwise van der Waals (vdW) and short-range interactions using the
-Morse potential and "(Grimme)"_#Grimme method implemented in the Many-Body
-Metal-Organic (MOMB) force field described comprehensively in "(Fichthorn)"_#Fichthorn and
-"(Zhou)"_#Zhou4. Grimme's method is widely used to correct for dispersion
-in density functional theory calculations.
+Style {momb} computes pairwise van der Waals (vdW) and short-range
+interactions using the Morse potential and "(Grimme)"_#Grimme method
+implemented in the Many-Body Metal-Organic (MOMB) force field
+described comprehensively in "(Fichthorn)"_#Fichthorn and
+"(Zhou)"_#Zhou4. Grimme's method is widely used to correct for
+dispersion in density functional theory calculations.
:c,image(Eqs/pair_momb.jpg)
-For the {momb} pair style, the following coefficients must be defined for each
-pair of atoms types via the "pair_coeff"_pair_coeff.html command as in the
-examples above, or in the data file or restart files read by the
-"read_data"_read_data.html as described below:
+For the {momb} pair style, the following coefficients must be defined
+for each pair of atoms types via the "pair_coeff"_pair_coeff.html
+command as in the examples above, or in the data file or restart files
+read by the "read_data"_read_data.html as described below:
D0 (energy units)
alpha (1/distance units)
diff --git a/examples/README b/examples/README
index 5880bb78d76762a56f261505f21683533afa07b5..090ed733ac792397cf0003e34e1b7faf81aca52e 100644
--- a/examples/README
+++ b/examples/README
@@ -63,8 +63,8 @@ balance: dynamic load balancing, 2d system
body: body particles, 2d system
cmap: CMAP 5-body contributions to CHARMM force field
colloid: big colloid particles in a small particle solvent, 2d system
-coreshell: adiabatic core/shell model
comb: models using the COMB potential
+coreshell: adiabatic core/shell model
controller: use of fix controller as a thermostat
crack: crack propagation in a 2d solid
deposit: deposition of atoms and molecules onto a 3d substrate
@@ -74,6 +74,7 @@ eim: NaCl using the EIM potential
ellipse: ellipsoidal particles in spherical solvent, 2d system
flow: Couette and Poiseuille flow in a 2d channel
friction: frictional contact of spherical asperities between 2d surfaces
+gcmc: Grand Canonical Monte Carlo (GCMC) via the fix gcmc command
granregion: use of fix wall/region/gran as boundary on granular particles
hugoniostat: Hugoniostat shock dynamics
indent: spherical indenter into a 2d solid
@@ -103,7 +104,7 @@ snap: NVE dynamics for BCC tantalum crystal using SNAP potential
streitz: Streitz-Mintmire potential for Al2O3
tad: temperature-accelerated dynamics of vacancy diffusion in bulk Si
vashishta: models using the Vashishta potential
-voronoi: test of Voronoi tesselation in compute voronoi/atom
+voronoi: Voronoi tesselation via compute voronoi/atom command
Here is a src/Make.py command which will perform a parallel build of a
LAMMPS executable "lmp_mpi" with all the packages needed by all the
diff --git a/examples/gcmc/CO2.txt b/examples/gcmc/CO2.txt
new file mode 100644
index 0000000000000000000000000000000000000000..72e593e3f2b1d9da040ba3e7bc9424e759f24fb7
--- /dev/null
+++ b/examples/gcmc/CO2.txt
@@ -0,0 +1,44 @@
+# CO2 molecule file. TraPPE model.
+
+3 atoms
+2 bonds
+1 angles
+
+Coords
+
+1 0.0 0.0 0.0
+2 -1.16 0.0 0.0
+3 1.16 0.0 0.0
+
+Types
+
+1 1
+2 2
+3 2
+
+Charges
+
+1 0.7
+2 -0.35
+3 -0.35
+
+Bonds
+
+1 1 1 2
+2 1 1 3
+
+Angles
+
+1 1 2 1 3
+
+Special Bond Counts
+
+1 2 0 0
+2 1 1 0
+3 1 1 0
+
+Special Bonds
+
+1 2 3
+2 1 3
+3 1 2
diff --git a/examples/gcmc/in.gcmc.co2 b/examples/gcmc/in.gcmc.co2
new file mode 100644
index 0000000000000000000000000000000000000000..0961e2b5568496b6c50e700eb635d83907ebb7d2
--- /dev/null
+++ b/examples/gcmc/in.gcmc.co2
@@ -0,0 +1,83 @@
+# GCMC for CO2 molecular fluid, rigid/small/nvt dynamics
+# Rigid CO2 TraPPE model
+# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
+# mixtures containing alkanes, carbon dioxide and
+# nitrogen AIChE J., 47,1676-1682 (2001)].
+
+# variables available on command line
+
+variable mu index -8.1
+variable disp index 0.5
+variable temp index 338.0
+variable lbox index 10.0
+variable spacing index 5.0
+
+# global model settings
+
+units real
+atom_style full
+boundary p p p
+pair_style lj/cut/coul/long 14
+pair_modify mix arithmetic tail yes
+kspace_style ewald 0.0001
+bond_style harmonic
+angle_style harmonic
+
+# box, start molecules on simple cubic lattice
+
+lattice sc ${spacing}
+region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
+create_box 2 box &
+ bond/types 1 &
+ angle/types 1 &
+ extra/bond/per/atom 2 &
+ extra/angle/per/atom 1 &
+ extra/special/per/atom 2
+molecule co2mol CO2.txt
+create_atoms 0 box mol co2mol 464563 units box
+
+# rigid CO2 TraPPE model
+
+pair_coeff 1 1 0.053649 2.8
+pair_coeff 2 2 0.156973 3.05
+bond_coeff 1 0 1.16
+angle_coeff 1 0 180
+
+# masses
+
+mass 1 12.0107
+mass 2 15.9994
+
+# MD settings
+
+group co2 type 1 2
+neighbor 2.0 bin
+neigh_modify every 1 delay 10 check yes
+velocity all create ${temp} 54654
+timestep 1.0
+
+# rigid constraints with thermostat
+
+fix myrigidnvt all rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
+fix_modify myrigidnvt dynamic/dof no
+
+# gcmc
+
+variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
+fix mygcmc all gcmc 100 100 100 0 54341 ${temp} ${mu} ${disp} mol &
+ co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+
+# output
+
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
+compute_modify thermo_temp dynamic/dof yes
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc
+thermo 1000
+
+# run
+
+run 20000
+
diff --git a/examples/gcmc/in.gcmc.lj b/examples/gcmc/in.gcmc.lj
new file mode 100644
index 0000000000000000000000000000000000000000..9b1f465f045ce2dc52e9880e326e4cde26bbf2a3
--- /dev/null
+++ b/examples/gcmc/in.gcmc.lj
@@ -0,0 +1,42 @@
+# GCMC for LJ simple fluid, no dynamics
+
+# variables available on command line
+
+variable mu index -21.0
+variable disp index 1.0
+variable temp index 2.0
+variable lbox index 10.0
+
+# global model settings
+
+units lj
+atom_style atomic
+pair_style lj/cut 3.0
+pair_modify tail yes
+
+# box
+
+region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
+create_box 1 box
+
+# lj parameters
+
+pair_coeff * * 1.0 1.0
+mass * 1.0
+
+# gcmc
+
+fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
+
+# output
+
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+compute_modify thermo_temp dynamic yes
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
+thermo 100
+
+# run
+
+run 1000
diff --git a/examples/gcmc/log.24Mar17.gcmc.co2.g++.1 b/examples/gcmc/log.24Mar17.gcmc.co2.g++.1
new file mode 100644
index 0000000000000000000000000000000000000000..7562476bf33d4639532dc538e00925adc9844975
--- /dev/null
+++ b/examples/gcmc/log.24Mar17.gcmc.co2.g++.1
@@ -0,0 +1,177 @@
+LAMMPS (17 Mar 2017)
+# GCMC for CO2 molecular fluid, rigid/small/nvt dynamics
+# Rigid CO2 TraPPE model
+# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
+# mixtures containing alkanes, carbon dioxide and
+# nitrogen AIChE J., 47,1676-1682 (2001)].
+
+# variables available on command line
+
+variable mu index -8.1
+variable disp index 0.5
+variable temp index 338.0
+variable lbox index 10.0
+variable spacing index 5.0
+
+# global model settings
+
+units real
+atom_style full
+boundary p p p
+pair_style lj/cut/coul/long 14
+pair_modify mix arithmetic tail yes
+kspace_style ewald 0.0001
+bond_style harmonic
+angle_style harmonic
+
+# box, start molecules on simple cubic lattice
+
+lattice sc ${spacing}
+lattice sc 5.0
+Lattice spacing in x,y,z = 5 5 5
+region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
+region box block 0 10.0 0 ${lbox} 0 ${lbox} units box
+region box block 0 10.0 0 10.0 0 ${lbox} units box
+region box block 0 10.0 0 10.0 0 10.0 units box
+create_box 2 box bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
+Created orthogonal box = (0 0 0) to (10 10 10)
+ 1 by 1 by 1 MPI processor grid
+molecule co2mol CO2.txt
+Read molecule co2mol:
+ 3 atoms with 2 types
+ 2 bonds with 1 types
+ 1 angles with 1 types
+ 0 dihedrals with 0 types
+ 0 impropers with 0 types
+create_atoms 0 box mol co2mol 464563 units box
+Created 24 atoms
+
+# rigid CO2 TraPPE model
+
+pair_coeff 1 1 0.053649 2.8
+pair_coeff 2 2 0.156973 3.05
+bond_coeff 1 0 1.16
+angle_coeff 1 0 180
+
+# masses
+
+mass 1 12.0107
+mass 2 15.9994
+
+# MD settings
+
+group co2 type 1 2
+24 atoms in group co2
+neighbor 2.0 bin
+neigh_modify every 1 delay 10 check yes
+velocity all create ${temp} 54654
+velocity all create 338.0 54654
+timestep 1.0
+
+# rigid constraints with thermostat
+
+fix myrigidnvt all rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
+fix myrigidnvt all rigid/nvt/small molecule temp 338.0 ${temp} 100 mol co2mol
+fix myrigidnvt all rigid/nvt/small molecule temp 338.0 338.0 100 mol co2mol
+8 rigid bodies with 24 atoms
+ 1.16 = max distance from body owner to body atom
+fix_modify myrigidnvt dynamic/dof no
+
+# gcmc
+
+variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
+fix mygcmc all gcmc 100 100 100 0 54341 ${temp} ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert 1.66666666666667 group co2 rigid myrigidnvt
+
+# output
+
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
+compute_modify thermo_temp dynamic/dof yes
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc
+thermo 1000
+
+# run
+
+run 20000
+Ewald initialization ...
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+ G vector (1/distance) = 0.164636
+ estimated absolute RMS force accuracy = 0.0332064
+ estimated relative force accuracy = 0.0001
+ KSpace vectors: actual max1d max3d = 16 2 62
+ kxmax kymax kzmax = 2 2 2
+WARNING: Fix gcmc using full_energy option (../fix_gcmc.cpp:439)
+0 atoms in group FixGCMC:gcmc_exclusion_group:mygcmc
+0 atoms in group FixGCMC:rotation_gas_atoms:mygcmc
+WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (../neighbor.cpp:472)
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 16
+ ghost atom cutoff = 16
+ binsize = 8, bins = 2 2 2
+ 1 neighbor lists, perpetual/occasional/extra = 1 0 0
+ (1) pair lj/cut/coul/long, perpetual
+ attributes: half, newton on
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+Per MPI rank memory allocation (min/avg/max) = 15.61 | 15.61 | 15.61 Mbytes
+Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_racc
+ 0 364.27579 4238.8631 -9.6809388 13.391989 0.5846359 24 0 0 0 0
+ 1000 311.39835 -327.93481 -8.6795381 9.9010062 0.51155641 21 0.13302848 0.12331626 0.6894397 0.90997852
+WARNING: Using kspace solver on system with no charge (../kspace.cpp:289)
+ 2000 905.66812 319.43347 -0.50350961 6.2991241 0.14615898 6 0.20952183 0.20430213 0.71797992 0.92626683
+ 3000 275.57393 -719.89718 -26.534978 14.238181 0.80387436 33 0.21291069 0.20460696 0.72899202 0.9133259
+ 4000 254.70771 -245.01902 -20.981537 13.160079 0.80387436 33 0.17245726 0.16974613 0.70145764 0.90542759
+ 5000 96.073601 -517.98124 -34.019065 5.441166 0.87695385 36 0.14174575 0.13607057 0.6776754 0.90155771
+ 6000 397.57265 148.92645 -7.2012893 10.665797 0.43847693 18 0.12299956 0.1202471 0.66165464 0.90274793
+ 7000 455.4271 -347.44181 -5.9244703 12.217875 0.43847693 18 0.15182038 0.14791307 0.67904236 0.90560829
+ 8000 301.03124 -627.45324 -13.251012 11.066909 0.5846359 24 0.16687346 0.16315516 0.6936719 0.91129375
+ 9000 256.5747 -565.67983 -17.814128 11.981874 0.73079488 30 0.15458482 0.15131825 0.68966283 0.90993975
+ 10000 443.60076 89.586912 -6.077863 11.900606 0.43847693 18 0.16092552 0.16020353 0.69882461 0.91422145
+ 11000 436.43777 64.412921 -6.7128469 11.708443 0.43847693 18 0.17453966 0.17480683 0.70679243 0.91369445
+ 12000 594.42207 849.07743 -3.3708621 10.040536 0.29231795 12 0.17461606 0.17568622 0.71175869 0.91333367
+ 13000 426.85849 -1093.1334 -17.524618 17.813377 0.65771539 27 0.17742896 0.17792831 0.71363306 0.91450124
+ 14000 317.75995 336.31107 -10.46774 11.681912 0.5846359 24 0.18331181 0.18427921 0.71715557 0.91652256
+ 15000 272.65129 317.50536 -26.428336 14.087176 0.80387436 33 0.17449167 0.175957 0.71122398 0.91528038
+ 16000 344.28567 -577.91079 -18.177927 16.077919 0.73079488 30 0.1661682 0.16781514 0.70485136 0.91508882
+ 17000 134.55928 -193.5668 -30.297136 7.6208177 0.87695385 36 0.15965609 0.1605036 0.69658104 0.9140445
+ 18000 231.87302 -446.07671 -14.875027 9.6763722 0.65771539 27 0.15270985 0.15351831 0.69002918 0.91372795
+ 19000 328.6835 -280.22365 -20.001303 16.982214 0.80387436 33 0.15201017 0.15272181 0.69023195 0.91272534
+ 20000 0 -20.39554 -0.14872889 -0 0 0 0.15600204 0.15750795 0.69503275 0.9138765
+Loop time of 30.9008 on 1 procs for 20000 steps with 0 atoms
+
+Performance: 55.921 ns/day, 0.429 hours/ns, 647.233 timesteps/s
+99.8% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 2.1985 | 2.1985 | 2.1985 | 0.0 | 7.11
+Bond | 0.029596 | 0.029596 | 0.029596 | 0.0 | 0.10
+Kspace | 0.23123 | 0.23123 | 0.23123 | 0.0 | 0.75
+Neigh | 0.16141 | 0.16141 | 0.16141 | 0.0 | 0.52
+Comm | 0.20628 | 0.20628 | 0.20628 | 0.0 | 0.67
+Output | 0.00068831 | 0.00068831 | 0.00068831 | 0.0 | 0.00
+Modify | 28.022 | 28.022 | 28.022 | 0.0 | 90.69
+Other | | 0.05058 | | | 0.16
+
+Nlocal: 0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost: 0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs: 0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 0
+Neighbor list builds = 40367
+Dangerous builds = 118
+
+Total wall time: 0:00:30
diff --git a/examples/gcmc/log.24Mar17.gcmc.co2.g++.4 b/examples/gcmc/log.24Mar17.gcmc.co2.g++.4
new file mode 100644
index 0000000000000000000000000000000000000000..65504b8d465cbbdb05ce8f355d81496ca8d2dd79
--- /dev/null
+++ b/examples/gcmc/log.24Mar17.gcmc.co2.g++.4
@@ -0,0 +1,179 @@
+LAMMPS (17 Mar 2017)
+# GCMC for CO2 molecular fluid, rigid/small/nvt dynamics
+# Rigid CO2 TraPPE model
+# [Potoff and J.I. Siepmann, Vapor-liquid equilibria of
+# mixtures containing alkanes, carbon dioxide and
+# nitrogen AIChE J., 47,1676-1682 (2001)].
+
+# variables available on command line
+
+variable mu index -8.1
+variable disp index 0.5
+variable temp index 338.0
+variable lbox index 10.0
+variable spacing index 5.0
+
+# global model settings
+
+units real
+atom_style full
+boundary p p p
+pair_style lj/cut/coul/long 14
+pair_modify mix arithmetic tail yes
+kspace_style ewald 0.0001
+bond_style harmonic
+angle_style harmonic
+
+# box, start molecules on simple cubic lattice
+
+lattice sc ${spacing}
+lattice sc 5.0
+Lattice spacing in x,y,z = 5 5 5
+region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} units box
+region box block 0 10.0 0 ${lbox} 0 ${lbox} units box
+region box block 0 10.0 0 10.0 0 ${lbox} units box
+region box block 0 10.0 0 10.0 0 10.0 units box
+create_box 2 box bond/types 1 angle/types 1 extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
+Created orthogonal box = (0 0 0) to (10 10 10)
+ 1 by 2 by 2 MPI processor grid
+molecule co2mol CO2.txt
+Read molecule co2mol:
+ 3 atoms with 2 types
+ 2 bonds with 1 types
+ 1 angles with 1 types
+ 0 dihedrals with 0 types
+ 0 impropers with 0 types
+create_atoms 0 box mol co2mol 464563 units box
+Created 24 atoms
+
+# rigid CO2 TraPPE model
+
+pair_coeff 1 1 0.053649 2.8
+pair_coeff 2 2 0.156973 3.05
+bond_coeff 1 0 1.16
+angle_coeff 1 0 180
+
+# masses
+
+mass 1 12.0107
+mass 2 15.9994
+
+# MD settings
+
+group co2 type 1 2
+24 atoms in group co2
+neighbor 2.0 bin
+neigh_modify every 1 delay 10 check yes
+velocity all create ${temp} 54654
+velocity all create 338.0 54654
+timestep 1.0
+
+# rigid constraints with thermostat
+
+fix myrigidnvt all rigid/nvt/small molecule temp ${temp} ${temp} 100 mol co2mol
+fix myrigidnvt all rigid/nvt/small molecule temp 338.0 ${temp} 100 mol co2mol
+fix myrigidnvt all rigid/nvt/small molecule temp 338.0 338.0 100 mol co2mol
+8 rigid bodies with 24 atoms
+ 1.16 = max distance from body owner to body atom
+fix_modify myrigidnvt dynamic/dof no
+
+# gcmc
+
+variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)
+fix mygcmc all gcmc 100 100 100 0 54341 ${temp} ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 ${mu} ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 ${disp} mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert ${tfac} group co2 rigid myrigidnvt
+fix mygcmc all gcmc 100 100 100 0 54341 338.0 -8.1 0.5 mol co2mol tfac_insert 1.66666666666667 group co2 rigid myrigidnvt
+
+# output
+
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
+compute_modify thermo_temp dynamic/dof yes
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc
+thermo 1000
+
+# run
+
+run 20000
+Ewald initialization ...
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+ G vector (1/distance) = 0.164636
+ estimated absolute RMS force accuracy = 0.0332064
+ estimated relative force accuracy = 0.0001
+ KSpace vectors: actual max1d max3d = 16 2 62
+ kxmax kymax kzmax = 2 2 2
+WARNING: Fix gcmc using full_energy option (../fix_gcmc.cpp:439)
+0 atoms in group FixGCMC:gcmc_exclusion_group:mygcmc
+0 atoms in group FixGCMC:rotation_gas_atoms:mygcmc
+WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (../neighbor.cpp:472)
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 16
+ ghost atom cutoff = 16
+ binsize = 8, bins = 2 2 2
+ 1 neighbor lists, perpetual/occasional/extra = 1 0 0
+ (1) pair lj/cut/coul/long, perpetual
+ attributes: half, newton on
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+Per MPI rank memory allocation (min/avg/max) = 15.4 | 15.4 | 15.4 Mbytes
+Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_racc
+ 0 386.52184 23582.465 -3.2433417 14.209828 0.5846359 24 0 0 0 0
+WARNING: Using kspace solver on system with no charge (../kspace.cpp:289)
+ 1000 760.80877 -39.270882 -3.5239626 12.851016 0.29231795 12 0.24161633 0.22984103 0.71087092 0.85283311
+ 2000 308.0739 -255.061 -20.411926 14.386853 0.73079488 30 0.26075352 0.24898725 0.73128383 0.88590474
+ 3000 432.34358 -1361.3278 -12.644057 15.894387 0.5846359 24 0.21121583 0.21051229 0.70036003 0.86735027
+ 4000 631.524 -63.488785 -3.6517158 13.804656 0.36539744 15 0.22486443 0.22886173 0.72358173 0.87172606
+ 5000 730.61244 -1029.284 -6.2144028 19.600352 0.43847693 18 0.23017182 0.22740779 0.72281887 0.87820845
+ 6000 752.43412 503.4547 -3.7053679 16.447663 0.36539744 15 0.22943971 0.226183 0.71450085 0.87447436
+ 7000 660.68448 828.51735 -10.592278 21.006666 0.51155641 21 0.24702096 0.24218506 0.71815602 0.8740222
+ 8000 331.58822 -621.22187 -5.3705759 7.2482776 0.36539744 15 0.23211903 0.22906813 0.70281376 0.86269411
+ 9000 413.91538 869.51669 -11.28701 15.216905 0.5846359 24 0.23246466 0.22923961 0.70832684 0.86244176
+ 10000 242.20861 -808.23311 -5.4533937 5.2945044 0.36539744 15 0.22024676 0.22031775 0.70785097 0.85712561
+ 11000 348.20046 -372.16895 -3.4663358 7.6114092 0.36539744 15 0.2252033 0.22688969 0.71513402 0.86123263
+ 12000 251.99682 303.30092 -18.58289 11.768089 0.73079488 30 0.20916844 0.21068047 0.694787 0.84635875
+ 13000 306.83592 -1582.0137 -20.810287 14.329041 0.73079488 30 0.19494837 0.196527 0.67554784 0.83056119
+ 14000 476.57411 268.94927 -14.185859 19.888076 0.65771539 27 0.19779631 0.20016859 0.67957528 0.83375167
+ 15000 267.03534 730.71183 -9.3348616 9.8171066 0.5846359 24 0.19468305 0.19814971 0.68032974 0.83810439
+ 16000 639.83235 2190.3244 -9.6666503 26.701062 0.65771539 27 0.19520687 0.19848997 0.68514387 0.84100361
+ 17000 2237.1203 -222.59057 -0.18248834 4.4456205 0.073079488 3 0.20412446 0.20757814 0.69175318 0.8434939
+ 18000 754.44841 205.54694 -10.501943 27.736031 0.5846359 24 0.2129422 0.21508015 0.69665031 0.84758331
+ 19000 1610.1148 1293.6003 -0.20849836 3.1996309 0.073079488 3 0.22061668 0.22356929 0.69949369 0.84851405
+ 20000 231.61458 -39.696514 -4.6452226 5.0629266 0.36539744 15 0.21984893 0.22246517 0.69914635 0.85058457
+Loop time of 21.1019 on 4 procs for 20000 steps with 15 atoms
+
+Performance: 81.888 ns/day, 0.293 hours/ns, 947.781 timesteps/s
+98.9% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 0.31897 | 0.41973 | 0.49748 | 10.1 | 1.99
+Bond | 0.014808 | 0.015063 | 0.015289 | 0.2 | 0.07
+Kspace | 0.3813 | 0.46228 | 0.56585 | 9.8 | 2.19
+Neigh | 0.049173 | 0.050043 | 0.050868 | 0.3 | 0.24
+Comm | 0.9755 | 0.99686 | 1.0205 | 1.9 | 4.72
+Output | 0.0014546 | 0.0015051 | 0.0016098 | 0.2 | 0.01
+Modify | 19.043 | 19.062 | 19.085 | 0.4 | 90.33
+Other | | 0.09438 | | | 0.45
+
+Nlocal: 3.75 ave 6 max 3 min
+Histogram: 3 0 0 0 0 0 0 0 0 1
+Nghost: 876.5 ave 937 max 818 min
+Histogram: 1 1 0 0 0 0 0 0 1 1
+Neighs: 490.5 ave 647 max 363 min
+Histogram: 1 0 1 0 0 1 0 0 0 1
+
+Total # of neighbors = 1962
+Ave neighs/atom = 130.8
+Ave special neighs/atom = 2
+Neighbor list builds = 40070
+Dangerous builds = 115
+
+Total wall time: 0:00:21
diff --git a/examples/gcmc/log.24Mar17.gcmc.lj.g++.1 b/examples/gcmc/log.24Mar17.gcmc.lj.g++.1
new file mode 100644
index 0000000000000000000000000000000000000000..5f967192e47dac52a4db4ead6e120bccc79fa62e
--- /dev/null
+++ b/examples/gcmc/log.24Mar17.gcmc.lj.g++.1
@@ -0,0 +1,103 @@
+LAMMPS (17 Mar 2017)
+# GCMC for LJ simple fluid, no dynamics
+
+# variables available on command line
+
+variable mu index -21.0
+variable disp index 1.0
+variable temp index 2.0
+variable lbox index 10.0
+
+# global model settings
+
+units lj
+atom_style atomic
+pair_style lj/cut 3.0
+pair_modify tail yes
+
+# box
+
+region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
+region box block 0 10.0 0 ${lbox} 0 ${lbox}
+region box block 0 10.0 0 10.0 0 ${lbox}
+region box block 0 10.0 0 10.0 0 10.0
+create_box 1 box
+Created orthogonal box = (0 0 0) to (10 10 10)
+ 1 by 1 by 1 MPI processor grid
+
+# lj parameters
+
+pair_coeff * * 1.0 1.0
+mass * 1.0
+
+# gcmc
+
+fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
+
+# output
+
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+compute_modify thermo_temp dynamic yes
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
+thermo 100
+
+# run
+
+run 1000
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 3.3
+ ghost atom cutoff = 3.3
+ binsize = 1.65, bins = 7 7 7
+ 1 neighbor lists, perpetual/occasional/extra = 1 0 0
+ (1) pair lj/cut, perpetual
+ attributes: half, newton on
+ pair build: half/bin/atomonly/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+Per MPI rank memory allocation (min/avg/max) = 0.4369 | 0.4369 | 0.4369 Mbytes
+Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
+ 0 0 0 0 -0 0 0 0 0 0
+ 100 1.9042848 0.39026453 -1.7692765 2.8466449 0.292 292 0.3619855 0.30247792 0.40278761
+ 200 1.8651924 0.47815517 -1.8494955 2.7886155 0.305 305 0.34021109 0.30357196 0.37759189
+ 300 2.0626994 0.52068504 -1.8197295 3.0834166 0.291 291 0.32055605 0.3003043 0.36103862
+ 400 2.0394818 0.53751435 -1.7636699 3.0482184 0.278 278 0.31698808 0.29995864 0.35441275
+ 500 1.9628066 0.54594742 -1.7145336 2.9339513 0.287 287 0.31211861 0.29724228 0.35161407
+ 600 1.9845913 0.40846162 -1.8199325 2.9669308 0.299 299 0.30976643 0.29612711 0.34933559
+ 700 1.8582606 0.53445462 -1.7869306 2.777974 0.296 296 0.30642103 0.29446478 0.34633665
+ 800 2.0340641 0.66057698 -1.7075279 3.0403148 0.283 283 0.30730979 0.29746793 0.34768045
+ 900 2.0830765 0.63731971 -1.894775 3.114911 0.322 322 0.30636338 0.29737705 0.34737644
+ 1000 1.9688933 0.50024802 -1.7013944 2.9428299 0.281 281 0.3053174 0.29772245 0.34788254
+Loop time of 3.98286 on 1 procs for 1000 steps with 281 atoms
+
+Performance: 108464.750 tau/day, 251.076 timesteps/s
+99.9% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 0.10563 | 0.10563 | 0.10563 | 0.0 | 2.65
+Neigh | 0.33428 | 0.33428 | 0.33428 | 0.0 | 8.39
+Comm | 0.027969 | 0.027969 | 0.027969 | 0.0 | 0.70
+Output | 0.00017285 | 0.00017285 | 0.00017285 | 0.0 | 0.00
+Modify | 3.5096 | 3.5096 | 3.5096 | 0.0 | 88.12
+Other | | 0.005197 | | | 0.13
+
+Nlocal: 281 ave 281 max 281 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost: 977 ave 977 max 977 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs: 5902 ave 5902 max 5902 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 5902
+Ave neighs/atom = 21.0036
+Neighbor list builds = 1000
+Dangerous builds = 0
+Total wall time: 0:00:03
diff --git a/examples/gcmc/log.24Mar17.gcmc.lj.g++.4 b/examples/gcmc/log.24Mar17.gcmc.lj.g++.4
new file mode 100644
index 0000000000000000000000000000000000000000..db4cc7fd9ac95c5dc2188f036e45341c7d1b574f
--- /dev/null
+++ b/examples/gcmc/log.24Mar17.gcmc.lj.g++.4
@@ -0,0 +1,103 @@
+LAMMPS (17 Mar 2017)
+# GCMC for LJ simple fluid, no dynamics
+
+# variables available on command line
+
+variable mu index -21.0
+variable disp index 1.0
+variable temp index 2.0
+variable lbox index 10.0
+
+# global model settings
+
+units lj
+atom_style atomic
+pair_style lj/cut 3.0
+pair_modify tail yes
+
+# box
+
+region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
+region box block 0 10.0 0 ${lbox} 0 ${lbox}
+region box block 0 10.0 0 10.0 0 ${lbox}
+region box block 0 10.0 0 10.0 0 10.0
+create_box 1 box
+Created orthogonal box = (0 0 0) to (10 10 10)
+ 1 by 2 by 2 MPI processor grid
+
+# lj parameters
+
+pair_coeff * * 1.0 1.0
+mass * 1.0
+
+# gcmc
+
+fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
+
+# output
+
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+compute_modify thermo_temp dynamic yes
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
+thermo 100
+
+# run
+
+run 1000
+Neighbor list info ...
+ update every 1 steps, delay 10 steps, check yes
+ max neighbors/atom: 2000, page size: 100000
+ master list distance cutoff = 3.3
+ ghost atom cutoff = 3.3
+ binsize = 1.65, bins = 7 7 7
+ 1 neighbor lists, perpetual/occasional/extra = 1 0 0
+ (1) pair lj/cut, perpetual
+ attributes: half, newton on
+ pair build: half/bin/atomonly/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+Per MPI rank memory allocation (min/avg/max) = 0.434 | 0.434 | 0.434 Mbytes
+Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
+ 0 0 0 0 -0 0 0 0 0 0
+ 100 2.0328045 0.58661762 -1.6812724 3.0385824 0.287 287 0.35917318 0.30067507 0.38663622
+ 200 1.9594279 0.50682399 -1.7308396 2.9287927 0.284 284 0.33788365 0.30337335 0.37300293
+ 300 2.0602937 0.7028247 -1.9278541 3.0806296 0.315 315 0.31882007 0.29697498 0.36167185
+ 400 1.995183 0.4328246 -1.8715454 2.983026 0.307 307 0.31527654 0.29681901 0.35673374
+ 500 2.1390101 0.48232215 -1.554138 3.1960306 0.257 257 0.31372975 0.30003067 0.35558858
+ 600 2.0584244 0.4929049 -1.6995569 3.0767263 0.283 283 0.31114213 0.29801665 0.35160109
+ 700 1.9155066 0.49654243 -1.5770611 2.8624174 0.265 265 0.31056419 0.29944173 0.35157337
+ 800 2.0883562 0.52731947 -1.8261112 3.1220925 0.3 300 0.30730979 0.29704354 0.34898892
+ 900 2.0470677 0.5605993 -2.0130053 3.0610656 0.322 322 0.30484441 0.29586719 0.34678883
+ 1000 2.004135 0.50642204 -1.6956257 2.9955798 0.283 283 0.30396929 0.29634309 0.34770304
+Loop time of 3.688 on 4 procs for 1000 steps with 283 atoms
+
+Performance: 117136.751 tau/day, 271.150 timesteps/s
+99.2% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section | min time | avg time | max time |%varavg| %total
+---------------------------------------------------------------
+Pair | 0.024644 | 0.026027 | 0.027483 | 0.6 | 0.71
+Neigh | 0.085449 | 0.088998 | 0.092893 | 0.9 | 2.41
+Comm | 0.045756 | 0.051296 | 0.056578 | 1.7 | 1.39
+Output | 0.00028491 | 0.00030857 | 0.00035262 | 0.0 | 0.01
+Modify | 3.5189 | 3.5191 | 3.5194 | 0.0 | 95.42
+Other | | 0.002221 | | | 0.06
+
+Nlocal: 70.75 ave 77 max 68 min
+Histogram: 1 2 0 0 0 0 0 0 0 1
+Nghost: 514.25 ave 520 max 507 min
+Histogram: 1 0 0 0 1 0 0 1 0 1
+Neighs: 1483.5 ave 1715 max 1359 min
+Histogram: 2 0 0 1 0 0 0 0 0 1
+
+Total # of neighbors = 5934
+Ave neighs/atom = 20.9682
+Neighbor list builds = 1000
+Dangerous builds = 0
+Total wall time: 0:00:03
diff --git a/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp b/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..11c7a147e76f41b005986a1685f8140d168ea426
--- /dev/null
+++ b/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
@@ -0,0 +1,1078 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, 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 author: Paul Crozier (SNL)
+ The lj-fsw sections (force-switched) were provided by
+ Robert Meissner and Lucio Colombi Ciacchi of
+ Bremen University, Bremen, Germany, with
+ additional assistance from Robert A. Latour, Clemson University
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_lj_charmmfsw_coul_long.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "kspace.h"
+#include "update.h"
+#include "integrate.h"
+#include "respa.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define EWALD_F 1.12837917
+#define EWALD_P 0.3275911
+#define A1 0.254829592
+#define A2 -0.284496736
+#define A3 1.421413741
+#define A4 -1.453152027
+#define A5 1.061405429
+
+/* ---------------------------------------------------------------------- */
+
+PairLJCharmmfswCoulLong::PairLJCharmmfswCoulLong(LAMMPS *lmp) : Pair(lmp)
+{
+ respa_enable = 1;
+ ewaldflag = pppmflag = 1;
+ ftable = NULL;
+ implicit = 0;
+ mix_flag = ARITHMETIC;
+ writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairLJCharmmfswCoulLong::~PairLJCharmmfswCoulLong()
+{
+ if (!copymode) {
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(epsilon);
+ memory->destroy(sigma);
+ memory->destroy(eps14);
+ memory->destroy(sigma14);
+ memory->destroy(lj1);
+ memory->destroy(lj2);
+ memory->destroy(lj3);
+ memory->destroy(lj4);
+ memory->destroy(lj14_1);
+ memory->destroy(lj14_2);
+ memory->destroy(lj14_3);
+ memory->destroy(lj14_4);
+ }
+ if (ftable) free_tables();
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype,itable;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwl12,evdwl6,ecoul,fpair;
+ double fraction,table;
+ double r,rinv,r2inv,r3inv,r6inv,rsq,forcecoul,forcelj,factor_coul,factor_lj;
+ double grij,expm2,prefactor,t,erfc;
+ double switch1;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = ecoul = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cut_bothsq) {
+ r2inv = 1.0/rsq;
+
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = qqrd2e * qtmp*q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ union_int_float_t rsq_lookup;
+ rsq_lookup.f = rsq;
+ itable = rsq_lookup.i & ncoulmask;
+ itable >>= ncoulshiftbits;
+ fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+ table = ftable[itable] + fraction*dftable[itable];
+ forcecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq) {
+ r6inv = r2inv*r2inv*r2inv;
+ jtype = type[j];
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+ } else forcelj = 0.0;
+
+ fpair = (forcecoul + factor_lj*forcelj) * r2inv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) {
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq)
+ ecoul = prefactor*erfc;
+ else {
+ table = etable[itable] + fraction*detable[itable];
+ ecoul = qtmp*q[j] * table;
+ }
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else ecoul = 0.0;
+
+ if (rsq < cut_ljsq) {
+ if (rsq > cut_lj_innersq) {
+ r = sqrt(rsq);
+ rinv = 1.0/r;
+ r3inv = rinv*rinv*rinv;
+ evdwl12 = lj3[itype][jtype]*cut_lj6*denom_lj12 *
+ (r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
+ evdwl6 = -lj4[itype][jtype]*cut_lj3*denom_lj6 *
+ (r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);;
+ evdwl = evdwl12 + evdwl6;
+ } else {
+ evdwl12 = r6inv*lj3[itype][jtype]*r6inv -
+ lj3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
+ evdwl6 = -lj4[itype][jtype]*r6inv +
+ lj4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
+ evdwl = evdwl12 + evdwl6;
+ }
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,ecoul,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::compute_inner()
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
+ double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double rsw;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = listinner->inum;
+ ilist = listinner->ilist;
+ numneigh = listinner->numneigh;
+ firstneigh = listinner->firstneigh;
+
+ double cut_out_on = cut_respa[0];
+ double cut_out_off = cut_respa[1];
+
+ double cut_out_diff = cut_out_off - cut_out_on;
+ double cut_out_on_sq = cut_out_on*cut_out_on;
+ double cut_out_off_sq = cut_out_off*cut_out_off;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cut_out_off_sq) {
+ r2inv = 1.0/rsq;
+ forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
+
+ r6inv = r2inv*r2inv*r2inv;
+ jtype = type[j];
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+
+ fpair = (forcecoul + factor_lj*forcelj) * r2inv;
+
+ if (rsq > cut_out_on_sq) {
+ rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+ fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+ }
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::compute_middle()
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
+ double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double switch1;
+ double rsw;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = listmiddle->inum;
+ ilist = listmiddle->ilist;
+ numneigh = listmiddle->numneigh;
+ firstneigh = listmiddle->firstneigh;
+
+ double cut_in_off = cut_respa[0];
+ double cut_in_on = cut_respa[1];
+ double cut_out_on = cut_respa[2];
+ double cut_out_off = cut_respa[3];
+
+ double cut_in_diff = cut_in_on - cut_in_off;
+ double cut_out_diff = cut_out_off - cut_out_on;
+ double cut_in_off_sq = cut_in_off*cut_in_off;
+ double cut_in_on_sq = cut_in_on*cut_in_on;
+ double cut_out_on_sq = cut_out_on*cut_out_on;
+ double cut_out_off_sq = cut_out_off*cut_out_off;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
+ r2inv = 1.0/rsq;
+ forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
+
+ r6inv = r2inv*r2inv*r2inv;
+ jtype = type[j];
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+
+ fpair = (forcecoul + factor_lj*forcelj) * r2inv;
+ if (rsq < cut_in_on_sq) {
+ rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+ fpair *= rsw*rsw*(3.0 - 2.0*rsw);
+ }
+ if (rsq > cut_out_on_sq) {
+ rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+ fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
+ }
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::compute_outer(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype,itable;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwl6,evdwl12,ecoul,fpair;
+ double fraction,table;
+ double r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double grij,expm2,prefactor,t,erfc;
+ double switch1;
+ double rsw;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double rsq;
+
+ evdwl = ecoul = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = listouter->inum;
+ ilist = listouter->ilist;
+ numneigh = listouter->numneigh;
+ firstneigh = listouter->firstneigh;
+
+ double cut_in_off = cut_respa[2];
+ double cut_in_on = cut_respa[3];
+
+ double cut_in_diff = cut_in_on - cut_in_off;
+ double cut_in_off_sq = cut_in_off*cut_in_off;
+ double cut_in_on_sq = cut_in_on*cut_in_on;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ if (rsq < cut_bothsq) {
+ r2inv = 1.0/rsq;
+
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = qqrd2e * qtmp*q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
+ if (rsq > cut_in_off_sq) {
+ if (rsq < cut_in_on_sq) {
+ rsw = (r - cut_in_off)/cut_in_diff;
+ forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw);
+ if (factor_coul < 1.0)
+ forcecoul -=
+ (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw);
+ } else {
+ forcecoul += prefactor;
+ if (factor_coul < 1.0)
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else {
+ union_int_float_t rsq_lookup;
+ rsq_lookup.f = rsq;
+ itable = rsq_lookup.i & ncoulmask;
+ itable >>= ncoulshiftbits;
+ fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+ table = ftable[itable] + fraction*dftable[itable];
+ forcecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq && rsq > cut_in_off_sq) {
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+ if (rsq < cut_in_on_sq) {
+ rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+ forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
+ }
+ } else forcelj = 0.0;
+
+ fpair = (forcecoul + forcelj) * r2inv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) {
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ ecoul = prefactor*erfc;
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ table = etable[itable] + fraction*detable[itable];
+ ecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ptable[itable] + fraction*dptable[itable];
+ prefactor = qtmp*q[j] * table;
+ ecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else ecoul = 0.0;
+
+ if (rsq < cut_ljsq) {
+ r6inv = r2inv*r2inv*r2inv;
+ evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ rinv = 1.0/r;
+ r3inv = rinv*rinv*rinv;
+ evdwl12 = lj3[itype][jtype]*cut_lj6*denom_lj12 *
+ (r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
+ evdwl6 = -lj4[itype][jtype]*cut_lj3*denom_lj6 *
+ (r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);;
+ evdwl = evdwl12 + evdwl6;
+ } else {
+ evdwl12 = r6inv*lj3[itype][jtype]*r6inv -
+ lj3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
+ evdwl6 = -lj4[itype][jtype]*r6inv +
+ lj4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
+ evdwl = evdwl12 + evdwl6;
+ }
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (vflag) {
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ table = vtable[itable] + fraction*dvtable[itable];
+ forcecoul = qtmp*q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ptable[itable] + fraction*dptable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+
+ if (rsq <= cut_in_off_sq) {
+ r6inv = r2inv*r2inv*r2inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+ } else if (rsq <= cut_in_on_sq) {
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+ }
+
+ fpair = (forcecoul + factor_lj*forcelj) * r2inv;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,ecoul,fpair,delx,dely,delz);
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(epsilon,n+1,n+1,"pair:epsilon");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+ memory->create(eps14,n+1,n+1,"pair:eps14");
+ memory->create(sigma14,n+1,n+1,"pair:sigma14");
+ memory->create(lj1,n+1,n+1,"pair:lj1");
+ memory->create(lj2,n+1,n+1,"pair:lj2");
+ memory->create(lj3,n+1,n+1,"pair:lj3");
+ memory->create(lj4,n+1,n+1,"pair:lj4");
+ memory->create(lj14_1,n+1,n+1,"pair:lj14_1");
+ memory->create(lj14_2,n+1,n+1,"pair:lj14_2");
+ memory->create(lj14_3,n+1,n+1,"pair:lj14_3");
+ memory->create(lj14_4,n+1,n+1,"pair:lj14_4");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+ unlike other pair styles,
+ there are no individual pair settings that these override
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::settings(int narg, char **arg)
+{
+ if (narg != 2 && narg != 3) error->all(FLERR,"Illegal pair_style command");
+
+ cut_lj_inner = force->numeric(FLERR,arg[0]);
+ cut_lj = force->numeric(FLERR,arg[1]);
+ if (narg == 2) cut_coul = cut_lj;
+ else cut_coul = force->numeric(FLERR,arg[2]);
+
+ // indicates pair_style being used for dihedral_charmm
+
+ dihedflag = 1;
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::coeff(int narg, char **arg)
+{
+ if (narg != 4 && narg != 6) error->all(FLERR,"Illegal pair_coeff command");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+ double epsilon_one = force->numeric(FLERR,arg[2]);
+ double sigma_one = force->numeric(FLERR,arg[3]);
+ double eps14_one = epsilon_one;
+ double sigma14_one = sigma_one;
+ if (narg == 6) {
+ eps14_one = force->numeric(FLERR,arg[4]);
+ sigma14_one = force->numeric(FLERR,arg[5]);
+ }
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ epsilon[i][j] = epsilon_one;
+ sigma[i][j] = sigma_one;
+ eps14[i][j] = eps14_one;
+ sigma14[i][j] = sigma14_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::init_style()
+{
+ if (!atom->q_flag)
+ error->all(FLERR,
+ "Pair style lj/charmmfsw/coul/long requires atom attribute q");
+
+ // request regular or rRESPA neighbor lists
+
+ int irequest;
+
+ if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
+ int respa = 0;
+ if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
+ if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
+
+ if (respa == 0) irequest = neighbor->request(this,instance_me);
+ else if (respa == 1) {
+ irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->id = 1;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->respainner = 1;
+ irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->id = 3;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->respaouter = 1;
+ } else {
+ irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->id = 1;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->respainner = 1;
+ irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->id = 2;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->respamiddle = 1;
+ irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->id = 3;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->respaouter = 1;
+ }
+
+ } else irequest = neighbor->request(this,instance_me);
+
+ // require cut_lj_inner < cut_lj
+
+ if (cut_lj_inner >= cut_lj)
+ error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
+
+ cut_lj_innersq = cut_lj_inner * cut_lj_inner;
+ cut_ljsq = cut_lj * cut_lj;
+ cut_ljinv = 1.0/cut_lj;
+ cut_lj_innerinv = 1.0/cut_lj_inner;
+ cut_lj3 = cut_lj * cut_lj * cut_lj;
+ cut_lj3inv = cut_ljinv * cut_ljinv * cut_ljinv;
+ cut_lj_inner3inv = cut_lj_innerinv * cut_lj_innerinv * cut_lj_innerinv;
+ cut_lj_inner3 = cut_lj_inner * cut_lj_inner * cut_lj_inner;
+ cut_lj6 = cut_ljsq * cut_ljsq * cut_ljsq;
+ cut_lj6inv = cut_lj3inv * cut_lj3inv;
+ cut_lj_inner6inv = cut_lj_inner3inv * cut_lj_inner3inv;
+ cut_lj_inner6 = cut_lj_innersq * cut_lj_innersq * cut_lj_innersq;
+ cut_coulsq = cut_coul * cut_coul;
+ cut_bothsq = MAX(cut_ljsq,cut_coulsq);
+
+ denom_lj = (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
+ (cut_ljsq-cut_lj_innersq);
+ denom_lj12 = 1.0/(cut_lj6 - cut_lj_inner6);
+ denom_lj6 = 1.0/(cut_lj3 - cut_lj_inner3);
+
+ // set & error check interior rRESPA cutoffs
+
+ if (strstr(update->integrate_style,"respa") &&
+ ((Respa *) update->integrate)->level_inner >= 0) {
+ cut_respa = ((Respa *) update->integrate)->cutoff;
+ if (MIN(cut_lj,cut_coul) < cut_respa[3])
+ error->all(FLERR,"Pair cutoff < Respa interior cutoff");
+ if (cut_lj_inner < cut_respa[1])
+ error->all(FLERR,"Pair inner cutoff < Respa interior cutoff");
+ } else cut_respa = NULL;
+
+ // insure use of KSpace long-range solver, set g_ewald
+
+ if (force->kspace == NULL)
+ error->all(FLERR,"Pair style requires a KSpace style");
+ g_ewald = force->kspace->g_ewald;
+
+ // setup force tables
+
+ if (ncoultablebits) init_tables(cut_coul,cut_respa);
+}
+
+/* ----------------------------------------------------------------------
+ neighbor callback to inform pair style of neighbor list to use
+ regular or rRESPA
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::init_list(int id, NeighList *ptr)
+{
+ if (id == 0) list = ptr;
+ else if (id == 1) listinner = ptr;
+ else if (id == 2) listmiddle = ptr;
+ else if (id == 3) listouter = ptr;
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairLJCharmmfswCoulLong::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) {
+ epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+ sigma[i][i],sigma[j][j]);
+ sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
+ eps14[i][j] = mix_energy(eps14[i][i],eps14[j][j],
+ sigma14[i][i],sigma14[j][j]);
+ sigma14[i][j] = mix_distance(sigma14[i][i],sigma14[j][j]);
+ }
+
+ double cut = MAX(cut_lj,cut_coul);
+
+ lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ lj14_1[i][j] = 48.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
+ lj14_2[i][j] = 24.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
+ lj14_3[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
+ lj14_4[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
+
+ lj1[j][i] = lj1[i][j];
+ lj2[j][i] = lj2[i][j];
+ lj3[j][i] = lj3[i][j];
+ lj4[j][i] = lj4[i][j];
+ lj14_1[j][i] = lj14_1[i][j];
+ lj14_2[j][i] = lj14_2[i][j];
+ lj14_3[j][i] = lj14_3[i][j];
+ lj14_4[j][i] = lj14_4[i][j];
+
+ return cut;
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) {
+ fwrite(&epsilon[i][j],sizeof(double),1,fp);
+ fwrite(&sigma[i][j],sizeof(double),1,fp);
+ fwrite(&eps14[i][j],sizeof(double),1,fp);
+ fwrite(&sigma14[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) {
+ fread(&epsilon[i][j],sizeof(double),1,fp);
+ fread(&sigma[i][j],sizeof(double),1,fp);
+ fread(&eps14[i][j],sizeof(double),1,fp);
+ fread(&sigma14[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&eps14[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma14[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::write_restart_settings(FILE *fp)
+{
+ fwrite(&cut_lj_inner,sizeof(double),1,fp);
+ fwrite(&cut_lj,sizeof(double),1,fp);
+ fwrite(&cut_coul,sizeof(double),1,fp);
+ fwrite(&offset_flag,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+ fwrite(&ncoultablebits,sizeof(int),1,fp);
+ fwrite(&tabinner,sizeof(double),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::read_restart_settings(FILE *fp)
+{
+ if (comm->me == 0) {
+ fread(&cut_lj_inner,sizeof(double),1,fp);
+ fread(&cut_lj,sizeof(double),1,fp);
+ fread(&cut_coul,sizeof(double),1,fp);
+ fread(&offset_flag,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ fread(&ncoultablebits,sizeof(int),1,fp);
+ fread(&tabinner,sizeof(double),1,fp);
+ }
+ MPI_Bcast(&cut_lj_inner,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_lj,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
+ MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
+}
+
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::write_data(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",
+ i,epsilon[i][i],sigma[i][i],eps14[i][i],sigma14[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulLong::write_data_all(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ for (int j = i; j <= atom->ntypes; j++)
+ fprintf(fp,"%d %d %g %g %g %g\n",i,j,
+ epsilon[i][j],sigma[i][j],eps14[i][j],sigma14[i][j]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJCharmmfswCoulLong::single(int i, int j, int itype, int jtype,
+ double rsq,
+ double factor_coul, double factor_lj,
+ double &fforce)
+{
+ double r,rinv,r2inv,r3inv,r6inv,grij,expm2,t,erfc,prefactor;
+ double switch1,fraction,table,forcecoul,forcelj,phicoul,philj,philj12,philj6;
+ int itable;
+
+ r = sqrt(rsq);
+ rinv = 1.0/r;
+ r2inv = 1.0/rsq;
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else {
+ union_int_float_t rsq_lookup;
+ rsq_lookup.f = rsq;
+ itable = rsq_lookup.i & ncoulmask;
+ itable >>= ncoulshiftbits;
+ fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+ table = ftable[itable] + fraction*dftable[itable];
+ forcecoul = atom->q[i]*atom->q[j] * table;
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = atom->q[i]*atom->q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 0.0;
+ if (rsq < cut_ljsq) {
+ r3inv = rinv*rinv*rinv;
+ r6inv = r3inv*r3inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+ } else forcelj = 0.0;
+ fforce = (forcecoul + factor_lj*forcelj) * r2inv;
+
+ double eng = 0.0;
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq)
+ phicoul = prefactor*erfc;
+ else {
+ table = etable[itable] + fraction*detable[itable];
+ phicoul = atom->q[i]*atom->q[j] * table;
+ }
+ if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
+ eng += phicoul;
+ }
+
+ if (rsq < cut_ljsq) {
+ if (rsq > cut_lj_innersq) {
+ philj12 = lj3[itype][jtype]*cut_lj6*denom_lj12 *
+ (r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
+ philj6 = -lj4[itype][jtype]*cut_lj3*denom_lj6 *
+ (r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);;
+ philj = philj12 + philj6;
+ } else {
+ philj12 = r6inv*lj3[itype][jtype]*r6inv -
+ lj3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
+ philj6 = -lj4[itype][jtype]*r6inv +
+ lj4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
+ philj = philj12 + philj6;
+ }
+ eng += factor_lj*philj;
+ }
+
+ return eng;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairLJCharmmfswCoulLong::extract(const char *str, int &dim)
+{
+ dim = 2;
+ if (strcmp(str,"lj14_1") == 0) return (void *) lj14_1;
+ if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2;
+ if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3;
+ if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4;
+
+ dim = 0;
+ if (strcmp(str,"implicit") == 0) return (void *) &implicit;
+
+ // info extracted by dihedral_charmmfsw
+
+ if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
+ if (strcmp(str,"cut_lj_inner") == 0) return (void *) &cut_lj_inner;
+ if (strcmp(str,"cut_lj") == 0) return (void *) &cut_lj;
+ if (strcmp(str,"dihedflag") == 0) return (void *) &dihedflag;
+
+ return NULL;
+}
diff --git a/src/KSPACE/pair_lj_charmmfsw_coul_long.h b/src/KSPACE/pair_lj_charmmfsw_coul_long.h
new file mode 100644
index 0000000000000000000000000000000000000000..650a908e4851dfe7ddb2102fa6b69961cdc727dc
--- /dev/null
+++ b/src/KSPACE/pair_lj_charmmfsw_coul_long.h
@@ -0,0 +1,110 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, 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 PAIR_CLASS
+
+PairStyle(lj/charmmfsw/coul/long,PairLJCharmmfswCoulLong)
+
+#else
+
+#ifndef LMP_PAIR_LJ_CHARMMFSW_COUL_LONG_H
+#define LMP_PAIR_LJ_CHARMMFSW_COUL_LONG_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairLJCharmmfswCoulLong : public Pair {
+ public:
+ PairLJCharmmfswCoulLong(class LAMMPS *);
+ virtual ~PairLJCharmmfswCoulLong();
+
+ virtual void compute(int, int);
+ virtual void settings(int, char **);
+ void coeff(int, char **);
+ virtual void init_style();
+ void init_list(int, class NeighList *);
+ virtual double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ void write_restart_settings(FILE *);
+ void read_restart_settings(FILE *);
+ void write_data(FILE *);
+ void write_data_all(FILE *);
+ virtual double single(int, int, int, int, double, double, double, double &);
+
+ void compute_inner();
+ void compute_middle();
+ virtual void compute_outer(int, int);
+ virtual void *extract(const char *, int &);
+
+ protected:
+ int implicit;
+ int dihedflag;
+
+ double cut_lj_inner,cut_lj,cut_ljinv,cut_lj_innerinv;
+ double cut_lj_innersq,cut_ljsq;
+ double cut_lj3inv,cut_lj_inner3inv,cut_lj3,cut_lj_inner3;
+ double cut_lj6inv,cut_lj_inner6inv,cut_lj6,cut_lj_inner6;
+ double cut_coul,cut_coulsq;
+ double cut_bothsq;
+ double denom_lj,denom_lj12,denom_lj6;
+ double **epsilon,**sigma,**eps14,**sigma14;
+ double **lj1,**lj2,**lj3,**lj4,**offset;
+ double **lj14_1,**lj14_2,**lj14_3,**lj14_4;
+ double *cut_respa;
+ double g_ewald;
+
+ virtual void allocate();
+};
+
+}
+
+#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.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair style lj/charmmfsw/coul/long requires atom attribute q
+
+The atom style defined does not have these attributes.
+
+E: Pair inner cutoff >= Pair outer cutoff
+
+The specified cutoffs for the pair style are inconsistent.
+
+E: Pair cutoff < Respa interior cutoff
+
+One or more pairwise cutoffs are too short to use with the specified
+rRESPA cutoffs.
+
+E: Pair inner cutoff < Respa interior cutoff
+
+One or more pairwise cutoffs are too short to use with the specified
+rRESPA cutoffs.
+
+E: Pair style requires a KSpace style
+
+No kspace style is defined.
+
+*/
diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp
index bf48c5d9aeba46eac255f6397923d21ef926d6bb..abf75c85c534ae3cdc3df88c3e28f513c0da3e52 100644
--- a/src/MANYBODY/pair_airebo.cpp
+++ b/src/MANYBODY/pair_airebo.cpp
@@ -3183,7 +3183,9 @@ double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj,
for (i=0; i<64; i++) coeffs[i]=0.0;
if (typei==0 && typej==0) {
- //if the inputs are out of bounds set them back to a point in bounds
+
+ // if the inputs are out of bounds set them back to a point in bounds
+
if (Nij<piCCdom[0][0]) Nij=piCCdom[0][0];
if (Nij>piCCdom[0][1]) Nij=piCCdom[0][1];
if (Nji<piCCdom[1][0]) Nji=piCCdom[1][0];
@@ -3213,10 +3215,10 @@ double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj,
}
}
-
// CH interaction
if ((typei==0 && typej==1) || (typei==1 && typej==0)) {
+
// if the inputs are out of bounds set them back to a point in bounds
if (Nij<piCHdom[0][0] || Nij>piCHdom[0][1] ||
diff --git a/src/MOLECULE/dihedral_charmmfsh.cpp b/src/MOLECULE/dihedral_charmmfsh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..93c1853fe53124bd6d9583be6317645631c1a147
--- /dev/null
+++ b/src/MOLECULE/dihedral_charmmfsh.cpp
@@ -0,0 +1,482 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, 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 author: Paul Crozier (SNL)
+ The force-shifted sections were provided by Robert Meissner
+ and Lucio Colombi Ciacchi of Bremen University, Bremen, Germany,
+ with additional assistance from Robert A. Latour, Clemson University
+------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include <math.h>
+#include <stdlib.h>
+#include "dihedral_charmmfsh.h"
+#include "atom.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "domain.h"
+#include "force.h"
+#include "pair.h"
+#include "update.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define TOLERANCE 0.05
+
+/* ---------------------------------------------------------------------- */
+
+DihedralCharmmfsh::DihedralCharmmfsh(LAMMPS *lmp) : Dihedral(lmp)
+{
+ weightflag = 0;
+ writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+DihedralCharmmfsh::~DihedralCharmmfsh()
+{
+ if (allocated && !copymode) {
+ memory->destroy(setflag);
+ memory->destroy(k);
+ memory->destroy(multiplicity);
+ memory->destroy(shift);
+ memory->destroy(cos_shift);
+ memory->destroy(sin_shift);
+ memory->destroy(weight);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DihedralCharmmfsh::compute(int eflag, int vflag)
+{
+ int i1,i2,i3,i4,i,m,n,type;
+ double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
+ double edihedral,f1[3],f2[3],f3[3],f4[3];
+ double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
+ double df,df1,ddf1,fg,hg,fga,hgb,gaa,gbb;
+ double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz;
+ double c,s,p,sx2,sy2,sz2;
+ int itype,jtype;
+ double delx,dely,delz,rsq,r2inv,r6inv,r;
+ double forcecoul,forcelj,fpair,ecoul,evdwl;
+
+ edihedral = evdwl = ecoul = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = 0;
+
+ // insure pair->ev_tally() will use 1-4 virial contribution
+
+ if (weightflag && vflag_global == 2)
+ force->pair->vflag_either = force->pair->vflag_global = 1;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *atomtype = atom->type;
+ int **dihedrallist = neighbor->dihedrallist;
+ int ndihedrallist = neighbor->ndihedrallist;
+ int nlocal = atom->nlocal;
+ int newton_bond = force->newton_bond;
+ double qqrd2e = force->qqrd2e;
+
+ for (n = 0; n < ndihedrallist; n++) {
+ i1 = dihedrallist[n][0];
+ i2 = dihedrallist[n][1];
+ i3 = dihedrallist[n][2];
+ i4 = dihedrallist[n][3];
+ type = dihedrallist[n][4];
+
+ // 1st bond
+
+ vb1x = x[i1][0] - x[i2][0];
+ vb1y = x[i1][1] - x[i2][1];
+ vb1z = x[i1][2] - x[i2][2];
+
+ // 2nd bond
+
+ vb2x = x[i3][0] - x[i2][0];
+ vb2y = x[i3][1] - x[i2][1];
+ vb2z = x[i3][2] - x[i2][2];
+
+ vb2xm = -vb2x;
+ vb2ym = -vb2y;
+ vb2zm = -vb2z;
+
+ // 3rd bond
+
+ vb3x = x[i4][0] - x[i3][0];
+ vb3y = x[i4][1] - x[i3][1];
+ vb3z = x[i4][2] - x[i3][2];
+
+ ax = vb1y*vb2zm - vb1z*vb2ym;
+ ay = vb1z*vb2xm - vb1x*vb2zm;
+ az = vb1x*vb2ym - vb1y*vb2xm;
+ bx = vb3y*vb2zm - vb3z*vb2ym;
+ by = vb3z*vb2xm - vb3x*vb2zm;
+ bz = vb3x*vb2ym - vb3y*vb2xm;
+
+ rasq = ax*ax + ay*ay + az*az;
+ rbsq = bx*bx + by*by + bz*bz;
+ rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
+ rg = sqrt(rgsq);
+
+ rginv = ra2inv = rb2inv = 0.0;
+ if (rg > 0) rginv = 1.0/rg;
+ if (rasq > 0) ra2inv = 1.0/rasq;
+ if (rbsq > 0) rb2inv = 1.0/rbsq;
+ rabinv = sqrt(ra2inv*rb2inv);
+
+ c = (ax*bx + ay*by + az*bz)*rabinv;
+ s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
+
+ // error check
+
+ if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
+ int me;
+ MPI_Comm_rank(world,&me);
+ if (screen) {
+ char str[128];
+ sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
+ TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT " " TAGINT_FORMAT,
+ me,update->ntimestep,
+ atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
+ error->warning(FLERR,str,0);
+ fprintf(screen," 1st atom: %d %g %g %g\n",
+ me,x[i1][0],x[i1][1],x[i1][2]);
+ fprintf(screen," 2nd atom: %d %g %g %g\n",
+ me,x[i2][0],x[i2][1],x[i2][2]);
+ fprintf(screen," 3rd atom: %d %g %g %g\n",
+ me,x[i3][0],x[i3][1],x[i3][2]);
+ fprintf(screen," 4th atom: %d %g %g %g\n",
+ me,x[i4][0],x[i4][1],x[i4][2]);
+ }
+ }
+
+ if (c > 1.0) c = 1.0;
+ if (c < -1.0) c = -1.0;
+
+ m = multiplicity[type];
+ p = 1.0;
+ ddf1 = df1 = 0.0;
+
+ for (i = 0; i < m; i++) {
+ ddf1 = p*c - df1*s;
+ df1 = p*s + df1*c;
+ p = ddf1;
+ }
+
+ p = p*cos_shift[type] + df1*sin_shift[type];
+ df1 = df1*cos_shift[type] - ddf1*sin_shift[type];
+ df1 *= -m;
+ p += 1.0;
+
+ if (m == 0) {
+ p = 1.0 + cos_shift[type];
+ df1 = 0.0;
+ }
+
+ if (eflag) edihedral = k[type] * p;
+
+ fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm;
+ hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm;
+ fga = fg*ra2inv*rginv;
+ hgb = hg*rb2inv*rginv;
+ gaa = -ra2inv*rg;
+ gbb = rb2inv*rg;
+
+ dtfx = gaa*ax;
+ dtfy = gaa*ay;
+ dtfz = gaa*az;
+ dtgx = fga*ax - hgb*bx;
+ dtgy = fga*ay - hgb*by;
+ dtgz = fga*az - hgb*bz;
+ dthx = gbb*bx;
+ dthy = gbb*by;
+ dthz = gbb*bz;
+
+ df = -k[type] * df1;
+
+ sx2 = df*dtgx;
+ sy2 = df*dtgy;
+ sz2 = df*dtgz;
+
+ f1[0] = df*dtfx;
+ f1[1] = df*dtfy;
+ f1[2] = df*dtfz;
+
+ f2[0] = sx2 - f1[0];
+ f2[1] = sy2 - f1[1];
+ f2[2] = sz2 - f1[2];
+
+ f4[0] = df*dthx;
+ f4[1] = df*dthy;
+ f4[2] = df*dthz;
+
+ f3[0] = -sx2 - f4[0];
+ f3[1] = -sy2 - f4[1];
+ f3[2] = -sz2 - f4[2];
+
+ // apply force to each of 4 atoms
+
+ if (newton_bond || i1 < nlocal) {
+ f[i1][0] += f1[0];
+ f[i1][1] += f1[1];
+ f[i1][2] += f1[2];
+ }
+
+ if (newton_bond || i2 < nlocal) {
+ f[i2][0] += f2[0];
+ f[i2][1] += f2[1];
+ f[i2][2] += f2[2];
+ }
+
+ if (newton_bond || i3 < nlocal) {
+ f[i3][0] += f3[0];
+ f[i3][1] += f3[1];
+ f[i3][2] += f3[2];
+ }
+
+ if (newton_bond || i4 < nlocal) {
+ f[i4][0] += f4[0];
+ f[i4][1] += f4[1];
+ f[i4][2] += f4[2];
+ }
+
+ if (evflag)
+ ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4,
+ vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
+
+ // 1-4 LJ and Coulomb interactions
+ // tally energy/virial in pair, using newton_bond as newton flag
+
+ if (weight[type] > 0.0) {
+ itype = atomtype[i1];
+ jtype = atomtype[i4];
+
+ delx = x[i1][0] - x[i4][0];
+ dely = x[i1][1] - x[i4][1];
+ delz = x[i1][2] - x[i4][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+
+ // modifying coul and LJ force and energies to apply
+ // force_shift and force_switch as in CHARMM pairwise
+ // LJ interactions between 1-4 atoms should usually be
+ // for r < cut_inner, so switching not applied
+
+ r = sqrt(rsq);
+ if (implicit) forcecoul = qqrd2e * q[i1]*q[i4]*r2inv;
+ else if (dihedflag) forcecoul = qqrd2e * q[i1]*q[i4]*sqrt(r2inv);
+ else forcecoul = qqrd2e * q[i1]*q[i4]*(sqrt(r2inv) -
+ r*cut_coulinv14*cut_coulinv14);
+ forcelj = r6inv * (lj14_1[itype][jtype]*r6inv - lj14_2[itype][jtype]);
+ fpair = weight[type] * (forcelj+forcecoul)*r2inv;
+
+ if (eflag) {
+ if (dihedflag) ecoul = weight[type] * forcecoul;
+ else ecoul = weight[type] * qqrd2e * q[i1]*q[i4] *
+ (sqrt(r2inv) + r*cut_coulinv14*cut_coulinv14 -
+ 2.0*cut_coulinv14);
+ evdwl14_12 = r6inv*lj14_3[itype][jtype]*r6inv -
+ lj14_3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
+ evdwl14_6 = -lj14_4[itype][jtype]*r6inv +
+ lj14_4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
+ evdwl = evdwl14_12 + evdwl14_6;
+ evdwl *= weight[type];
+ }
+
+ if (newton_bond || i1 < nlocal) {
+ f[i1][0] += delx*fpair;
+ f[i1][1] += dely*fpair;
+ f[i1][2] += delz*fpair;
+ }
+ if (newton_bond || i4 < nlocal) {
+ f[i4][0] -= delx*fpair;
+ f[i4][1] -= dely*fpair;
+ f[i4][2] -= delz*fpair;
+ }
+
+ if (evflag) force->pair->ev_tally(i1,i4,nlocal,newton_bond,
+ evdwl,ecoul,fpair,delx,dely,delz);
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DihedralCharmmfsh::allocate()
+{
+ allocated = 1;
+ int n = atom->ndihedraltypes;
+
+ memory->create(k,n+1,"dihedral:k");
+ memory->create(multiplicity,n+1,"dihedral:k");
+ memory->create(shift,n+1,"dihedral:shift");
+ memory->create(cos_shift,n+1,"dihedral:cos_shift");
+ memory->create(sin_shift,n+1,"dihedral:sin_shift");
+ memory->create(weight,n+1,"dihedral:weight");
+
+ memory->create(setflag,n+1,"dihedral:setflag");
+ for (int i = 1; i <= n; i++) setflag[i] = 0;
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one type
+------------------------------------------------------------------------- */
+
+void DihedralCharmmfsh::coeff(int narg, char **arg)
+{
+ if (narg != 5) error->all(FLERR,"Incorrect args for dihedral coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi;
+ force->bounds(FLERR,arg[0],atom->ndihedraltypes,ilo,ihi);
+
+ // require integer values of shift for backwards compatibility
+ // arbitrary phase angle shift could be allowed, but would break
+ // backwards compatibility and is probably not needed
+
+ double k_one = force->numeric(FLERR,arg[1]);
+ int multiplicity_one = force->inumeric(FLERR,arg[2]);
+ int shift_one = force->inumeric(FLERR,arg[3]);
+ double weight_one = force->numeric(FLERR,arg[4]);
+
+ if (multiplicity_one < 0)
+ error->all(FLERR,"Incorrect multiplicity arg for dihedral coefficients");
+ if (weight_one < 0.0 || weight_one > 1.0)
+ error->all(FLERR,"Incorrect weight arg for dihedral coefficients");
+ if (weight_one > 0.0) weightflag=1;
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ k[i] = k_one;
+ shift[i] = shift_one;
+ cos_shift[i] = cos(MY_PI*shift_one/180.0);
+ sin_shift[i] = sin(MY_PI*shift_one/180.0);
+ multiplicity[i] = multiplicity_one;
+ weight[i] = weight_one;
+ setflag[i] = 1;
+ count++;
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ error check and initialize all values needed for force computation
+------------------------------------------------------------------------- */
+
+void DihedralCharmmfsh::init_style()
+{
+ // insure use of CHARMM pair_style if any weight factors are non-zero
+ // set local ptrs to LJ 14 arrays setup by Pair
+
+ if (weightflag) {
+ int itmp;
+ if (force->pair == NULL)
+ error->all(FLERR,"Dihedral charmmfsh is incompatible with Pair style");
+ lj14_1 = (double **) force->pair->extract("lj14_1",itmp);
+ lj14_2 = (double **) force->pair->extract("lj14_2",itmp);
+ lj14_3 = (double **) force->pair->extract("lj14_3",itmp);
+ lj14_4 = (double **) force->pair->extract("lj14_4",itmp);
+ int *ptr = (int *) force->pair->extract("implicit",itmp);
+ if (!lj14_1 || !lj14_2 || !lj14_3 || !lj14_4 || !ptr)
+ error->all(FLERR,"Dihedral charmmfsh is incompatible with Pair style");
+ implicit = *ptr;
+ }
+
+ // constants for applying force switch (LJ) and force_shift (coul)
+ // to 1/4 dihedral atoms to match CHARMM pairwise interactions
+
+ int itmp;
+ int *p_dihedflag = (int *) force->pair->extract("dihedflag",itmp);
+ double *p_cutljinner = (double *) force->pair->extract("cut_lj_inner",itmp);
+ double *p_cutlj = (double *) force->pair->extract("cut_lj",itmp);
+ double *p_cutcoul = (double *) force->pair->extract("cut_coul",itmp);
+
+ if (p_cutcoul == NULL || p_cutljinner == NULL ||
+ p_cutlj == NULL || p_dihedflag == NULL)
+ error->all(FLERR,"Dihedral charmmfsh is incompatible with Pair style");
+
+ dihedflag = *p_dihedflag;
+ cut_coul14 = *p_cutcoul;
+ cut_lj_inner14 = *p_cutljinner;
+ cut_lj14 = *p_cutlj;
+
+ cut_coulinv14 = 1/cut_coul14;
+ cut_lj_inner3inv = (1/cut_lj_inner14) * (1/cut_lj_inner14) *
+ (1/cut_lj_inner14);
+ cut_lj_inner6inv = cut_lj_inner3inv * cut_lj_inner3inv;
+ cut_lj3inv = (1/cut_lj14) * (1/cut_lj14) * (1/cut_lj14);
+ cut_lj6inv = cut_lj3inv * cut_lj3inv;
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes out coeffs to restart file
+------------------------------------------------------------------------- */
+
+void DihedralCharmmfsh::write_restart(FILE *fp)
+{
+ fwrite(&k[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp);
+ fwrite(&shift[1],sizeof(int),atom->ndihedraltypes,fp);
+ fwrite(&weight[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&weightflag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads coeffs from restart file, bcasts them
+------------------------------------------------------------------------- */
+
+void DihedralCharmmfsh::read_restart(FILE *fp)
+{
+ allocate();
+
+ if (comm->me == 0) {
+ fread(&k[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&multiplicity[1],sizeof(int),atom->ndihedraltypes,fp);
+ fread(&shift[1],sizeof(int),atom->ndihedraltypes,fp);
+ fread(&weight[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&weightflag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&multiplicity[1],atom->ndihedraltypes,MPI_INT,0,world);
+ MPI_Bcast(&shift[1],atom->ndihedraltypes,MPI_INT,0,world);
+ MPI_Bcast(&weight[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&weightflag,1,MPI_INT,0,world);
+
+ for (int i = 1; i <= atom->ndihedraltypes; i++) {
+ setflag[i] = 1;
+ cos_shift[i] = cos(MY_PI*shift[i]/180.0);
+ sin_shift[i] = sin(MY_PI*shift[i]/180.0);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void DihedralCharmmfsh::write_data(FILE *fp)
+{
+ for (int i = 1; i <= atom->ndihedraltypes; i++)
+ fprintf(fp,"%d %g %d %d %g\n",i,k[i],multiplicity[i],shift[i],weight[i]);
+}
+
diff --git a/src/MOLECULE/dihedral_charmmfsh.h b/src/MOLECULE/dihedral_charmmfsh.h
new file mode 100644
index 0000000000000000000000000000000000000000..44ea9b265860de23625ad4daf0d2aabd074ffc78
--- /dev/null
+++ b/src/MOLECULE/dihedral_charmmfsh.h
@@ -0,0 +1,81 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, 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 DIHEDRAL_CLASS
+
+DihedralStyle(charmmfsh,DihedralCharmmfsh)
+
+#else
+
+#ifndef LMP_DIHEDRAL_CHARMMFSH_H
+#define LMP_DIHEDRAL_CHARMMFSH_H
+
+#include <stdio.h>
+#include "dihedral.h"
+
+namespace LAMMPS_NS {
+
+class DihedralCharmmfsh : public Dihedral {
+ public:
+ DihedralCharmmfsh(class LAMMPS *);
+ virtual ~DihedralCharmmfsh();
+ virtual void compute(int, int);
+ virtual void coeff(int, char **);
+ virtual void init_style();
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ void write_data(FILE *);
+
+ protected:
+ int implicit,weightflag,dihedflag;
+ double cut_lj_inner14,cut_lj14,cut_coul14;
+ double evdwl14_12,evdwl14_6,cut_coulinv14;
+ double cut_lj_inner3inv,cut_lj_inner6inv,cut_lj3inv,cut_lj6inv;
+
+ double *k,*weight,*cos_shift,*sin_shift;
+ int *multiplicity,*shift;
+ double **lj14_1,**lj14_2,**lj14_3,**lj14_4;
+
+ virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+W: Dihedral problem: %d %ld %d %d %d %d
+
+Conformation of the 4 listed dihedral atoms is extreme; you may want
+to check your simulation geometry.
+
+E: Incorrect args for dihedral coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Incorrect multiplicity arg for dihedral coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Incorrect weight arg for dihedral coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Dihedral charmmfsh is incompatible with Pair style
+
+Dihedral style charmmfsh must be used with a pair style charmm
+in order for the 1-4 epsilon/sigma parameters to be defined.
+
+*/
diff --git a/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp b/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c75da63caedcaf42d5fcc8b079cd2a0a457b57b6
--- /dev/null
+++ b/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp
@@ -0,0 +1,546 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, 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 author: Paul Crozier (SNL)
+ The lj-fsw/coul-fsh (force-switched and force-shifted) sections
+ were provided by Robert Meissner
+ and Lucio Colombi Ciacchi of Bremen University, Bremen, Germany,
+ with additional assistance from Robert A. Latour, Clemson University
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_lj_charmmfsw_coul_charmmfsh.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairLJCharmmfswCoulCharmmfsh::PairLJCharmmfswCoulCharmmfsh(LAMMPS *lmp) :
+ Pair(lmp)
+{
+ implicit = 0;
+ mix_flag = ARITHMETIC;
+ writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairLJCharmmfswCoulCharmmfsh::~PairLJCharmmfswCoulCharmmfsh()
+{
+ if (!copymode) {
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(epsilon);
+ memory->destroy(sigma);
+ memory->destroy(eps14);
+ memory->destroy(sigma14);
+ memory->destroy(lj1);
+ memory->destroy(lj2);
+ memory->destroy(lj3);
+ memory->destroy(lj4);
+ memory->destroy(lj14_1);
+ memory->destroy(lj14_2);
+ memory->destroy(lj14_3);
+ memory->destroy(lj14_4);
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwl12,evdwl6,ecoul,fpair;
+ double r,rinv,r3inv,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double switch1;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = ecoul = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cut_bothsq) {
+ r2inv = 1.0/rsq;
+ r = sqrt(rsq);
+
+ if (rsq < cut_coulsq) {
+ forcecoul = qqrd2e * qtmp*q[j]*
+ (sqrt(r2inv) - r*cut_coulinv*cut_coulinv);
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq) {
+ r6inv = r2inv*r2inv*r2inv;
+ jtype = type[j];
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+ } else forcelj = 0.0;
+
+ fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) {
+ if (rsq < cut_coulsq) {
+ ecoul = qqrd2e * qtmp*q[j]*
+ (sqrt(r2inv) + cut_coulinv*cut_coulinv*r - 2.0*cut_coulinv);
+ ecoul *= factor_coul;
+ } else ecoul = 0.0;
+ if (rsq < cut_ljsq) {
+ if (rsq > cut_lj_innersq) {
+ rinv = 1.0/r;
+ r3inv = rinv*rinv*rinv;
+ evdwl12 = lj3[itype][jtype]*cut_lj6*denom_lj12 *
+ (r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
+ evdwl6 = -lj4[itype][jtype]*cut_lj3*denom_lj6 *
+ (r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);;
+ evdwl = evdwl12 + evdwl6;
+ } else {
+ evdwl12 = r6inv*lj3[itype][jtype]*r6inv -
+ lj3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
+ evdwl6 = -lj4[itype][jtype]*r6inv +
+ lj4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
+ evdwl = evdwl12 + evdwl6;
+ }
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,ecoul,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(epsilon,n+1,n+1,"pair:epsilon");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+ memory->create(eps14,n+1,n+1,"pair:eps14");
+ memory->create(sigma14,n+1,n+1,"pair:sigma14");
+ memory->create(lj1,n+1,n+1,"pair:lj1");
+ memory->create(lj2,n+1,n+1,"pair:lj2");
+ memory->create(lj3,n+1,n+1,"pair:lj3");
+ memory->create(lj4,n+1,n+1,"pair:lj4");
+ memory->create(lj14_1,n+1,n+1,"pair:lj14_1");
+ memory->create(lj14_2,n+1,n+1,"pair:lj14_2");
+ memory->create(lj14_3,n+1,n+1,"pair:lj14_3");
+ memory->create(lj14_4,n+1,n+1,"pair:lj14_4");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+ unlike other pair styles,
+ there are no individual pair settings that these override
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::settings(int narg, char **arg)
+{
+ if (narg != 2 && narg != 3)
+ error->all(FLERR,"Illegal pair_style command");
+
+ cut_lj_inner = force->numeric(FLERR,arg[0]);
+ cut_lj = force->numeric(FLERR,arg[1]);
+ if (narg == 2) {
+ cut_coul = cut_lj;
+ } else {
+ cut_coul = force->numeric(FLERR,arg[2]);
+ }
+
+ // indicates pair_style being used for dihedral_charmm
+
+ dihedflag = 0;
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::coeff(int narg, char **arg)
+{
+ if (narg != 4 && narg != 6)
+ error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+ double epsilon_one = force->numeric(FLERR,arg[2]);
+ double sigma_one = force->numeric(FLERR,arg[3]);
+ double eps14_one = epsilon_one;
+ double sigma14_one = sigma_one;
+ if (narg == 6) {
+ eps14_one = force->numeric(FLERR,arg[4]);
+ sigma14_one = force->numeric(FLERR,arg[5]);
+ }
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ epsilon[i][j] = epsilon_one;
+ sigma[i][j] = sigma_one;
+ eps14[i][j] = eps14_one;
+ sigma14[i][j] = sigma14_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::init_style()
+{
+ if (!atom->q_flag)
+ error->all(FLERR,"Pair style lj/charmmfsw/coul/charmmfsh "
+ "requires atom attribute q");
+
+ neighbor->request(this,instance_me);
+
+ // require cut_lj_inner < cut_lj
+
+ if (cut_lj_inner >= cut_lj)
+ error->all(FLERR,"Pair inner lj cutoff >= Pair outer lj cutoff");
+
+ cut_lj_innersq = cut_lj_inner * cut_lj_inner;
+ cut_ljsq = cut_lj * cut_lj;
+ cut_ljinv = 1.0/cut_lj;
+ cut_lj_innerinv = 1.0/cut_lj_inner;
+ cut_lj3 = cut_lj * cut_lj * cut_lj;
+ cut_lj3inv = cut_ljinv * cut_ljinv * cut_ljinv;
+ cut_lj_inner3inv = cut_lj_innerinv * cut_lj_innerinv * cut_lj_innerinv;
+ cut_lj_inner3 = cut_lj_inner * cut_lj_inner * cut_lj_inner;
+ cut_lj6 = cut_ljsq * cut_ljsq * cut_ljsq;
+ cut_lj6inv = cut_lj3inv * cut_lj3inv;
+ cut_lj_inner6inv = cut_lj_inner3inv * cut_lj_inner3inv;
+ cut_lj_inner6 = cut_lj_innersq * cut_lj_innersq * cut_lj_innersq;
+ cut_coulsq = cut_coul * cut_coul;
+ cut_coulinv = 1.0/cut_coul;
+ cut_bothsq = MAX(cut_ljsq,cut_coulsq);
+
+ denom_lj = (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
+ (cut_ljsq-cut_lj_innersq);
+ denom_lj12 = 1.0/(cut_lj6 - cut_lj_inner6);
+ denom_lj6 = 1.0/(cut_lj3 - cut_lj_inner3);
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairLJCharmmfswCoulCharmmfsh::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) {
+ epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+ sigma[i][i],sigma[j][j]);
+ sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
+ eps14[i][j] = mix_energy(eps14[i][i],eps14[j][j],
+ sigma14[i][i],sigma14[j][j]);
+ sigma14[i][j] = mix_distance(sigma14[i][i],sigma14[j][j]);
+ }
+
+ double cut = MAX(cut_lj,cut_coul);
+
+ lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ lj14_1[i][j] = 48.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
+ lj14_2[i][j] = 24.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
+ lj14_3[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],12.0);
+ lj14_4[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],6.0);
+
+ lj1[j][i] = lj1[i][j];
+ lj2[j][i] = lj2[i][j];
+ lj3[j][i] = lj3[i][j];
+ lj4[j][i] = lj4[i][j];
+ lj14_1[j][i] = lj14_1[i][j];
+ lj14_2[j][i] = lj14_2[i][j];
+ lj14_3[j][i] = lj14_3[i][j];
+ lj14_4[j][i] = lj14_4[i][j];
+
+ return cut;
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::write_data(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ fprintf(fp,"%d %g %g %g %g\n",
+ i,epsilon[i][i],sigma[i][i],eps14[i][i],sigma14[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::write_data_all(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ for (int j = i; j <= atom->ntypes; j++)
+ fprintf(fp,"%d %d %g %g %g %g\n",i,j,
+ epsilon[i][j],sigma[i][j],eps14[i][j],sigma14[i][j]);
+}
+
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) {
+ fwrite(&epsilon[i][j],sizeof(double),1,fp);
+ fwrite(&sigma[i][j],sizeof(double),1,fp);
+ fwrite(&eps14[i][j],sizeof(double),1,fp);
+ fwrite(&sigma14[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) {
+ fread(&epsilon[i][j],sizeof(double),1,fp);
+ fread(&sigma[i][j],sizeof(double),1,fp);
+ fread(&eps14[i][j],sizeof(double),1,fp);
+ fread(&sigma14[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&eps14[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma14[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::write_restart_settings(FILE *fp)
+{
+ fwrite(&cut_lj_inner,sizeof(double),1,fp);
+ fwrite(&cut_lj,sizeof(double),1,fp);
+ fwrite(&cut_coul,sizeof(double),1,fp);
+ fwrite(&offset_flag,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJCharmmfswCoulCharmmfsh::read_restart_settings(FILE *fp)
+{
+ if (comm->me == 0) {
+ fread(&cut_lj_inner,sizeof(double),1,fp);
+ fread(&cut_lj,sizeof(double),1,fp);
+ fread(&cut_coul,sizeof(double),1,fp);
+ fread(&offset_flag,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&cut_lj_inner,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_lj,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJCharmmfswCoulCharmmfsh::
+single(int i, int j, int itype, int jtype,
+ double rsq, double factor_coul, double factor_lj, double &fforce)
+{
+ double r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;
+ double phicoul,philj,philj12,philj6;
+ double switch1;
+
+ r2inv = 1.0/rsq;
+ r = sqrt(rsq);
+ rinv = 1.0/r;
+ if (rsq < cut_coulsq) {
+ forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] *
+ (sqrt(r2inv) - r*cut_coulinv*cut_coulinv);
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq) {
+ r6inv = r2inv*r2inv*r2inv;
+ r3inv = rinv*rinv*rinv;
+ forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ if (rsq > cut_lj_innersq) {
+ switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+ (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+ forcelj = forcelj*switch1;
+ }
+ } else forcelj = 0.0;
+
+ fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
+
+ double eng = 0.0;
+ if (rsq < cut_coulsq) {
+ phicoul = force->qqrd2e * atom->q[i]*atom->q[j] *
+ (sqrt(r2inv) + cut_coulinv*cut_coulinv*r - 2.0*cut_coulinv);
+ eng += factor_coul*phicoul;
+ }
+ if (rsq < cut_ljsq) {
+ if (rsq > cut_lj_innersq) {
+ philj12 = lj3[itype][jtype]*cut_lj6*denom_lj12 *
+ (r6inv - cut_lj6inv)*(r6inv - cut_lj6inv);
+ philj6 = -lj4[itype][jtype]*cut_lj3*denom_lj6 *
+ (r3inv - cut_lj3inv)*(r3inv - cut_lj3inv);;
+ philj = philj12 + philj6;
+ } else {
+ philj12 = r6inv*lj3[itype][jtype]*r6inv -
+ lj3[itype][jtype]*cut_lj_inner6inv*cut_lj6inv;
+ philj6 = -lj4[itype][jtype]*r6inv +
+ lj4[itype][jtype]*cut_lj_inner3inv*cut_lj3inv;
+ philj = philj12 + philj6;
+ }
+ eng += factor_lj*philj;
+ }
+
+ return eng;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairLJCharmmfswCoulCharmmfsh::extract(const char *str, int &dim)
+{
+ dim = 2;
+ if (strcmp(str,"lj14_1") == 0) return (void *) lj14_1;
+ if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2;
+ if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3;
+ if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4;
+
+ dim = 0;
+ if (strcmp(str,"implicit") == 0) return (void *) &implicit;
+
+ // info extracted by dihedral_charmmf
+
+ if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
+ if (strcmp(str,"cut_lj_inner") == 0) return (void *) &cut_lj_inner;
+ if (strcmp(str,"cut_lj") == 0) return (void *) &cut_lj;
+ if (strcmp(str,"dihedflag") == 0) return (void *) &dihedflag;
+
+ return NULL;
+}
diff --git a/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.h b/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.h
new file mode 100644
index 0000000000000000000000000000000000000000..c18036caea43bd5fcd059c9c35e556bdc6b360a5
--- /dev/null
+++ b/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.h
@@ -0,0 +1,88 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, 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 PAIR_CLASS
+
+PairStyle(lj/charmmfsw/coul/charmmfsh,PairLJCharmmfswCoulCharmmfsh)
+
+#else
+
+#ifndef LMP_PAIR_LJ_CHARMMFSW_COUL_CHARMMFSH_H
+#define LMP_PAIR_LJ_CHARMMFSW_COUL_CHARMMFSH_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairLJCharmmfswCoulCharmmfsh : public Pair {
+ public:
+ PairLJCharmmfswCoulCharmmfsh(class LAMMPS *);
+ virtual ~PairLJCharmmfswCoulCharmmfsh();
+ virtual void compute(int, int);
+ virtual void settings(int, char **);
+ void coeff(int, char **);
+ virtual void init_style();
+ virtual double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ void write_restart_settings(FILE *);
+ void read_restart_settings(FILE *);
+ void write_data(FILE *);
+ void write_data_all(FILE *);
+ virtual double single(int, int, int, int, double, double, double, double &);
+ virtual void *extract(const char *, int &);
+
+ protected:
+ int implicit;
+ int dihedflag;
+
+ double cut_lj_inner,cut_lj,cut_coul,cut_coulinv,cut_ljinv,cut_lj_innerinv;
+ double cut_lj_innersq,cut_ljsq,cut_coulsq,cut_bothsq;
+ double cut_lj3inv,cut_lj_inner3inv,cut_lj3,cut_lj_inner3;
+ double cut_lj6inv,cut_lj_inner6inv,cut_lj6,cut_lj_inner6;
+ double denom_lj,denom_lj12,denom_lj6;
+
+ double **epsilon,**sigma,**eps14,**sigma14;
+ double **lj1,**lj2,**lj3,**lj4;
+ double **lj14_1,**lj14_2,**lj14_3,**lj14_4;
+
+ virtual void allocate();
+};
+
+}
+
+#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.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair style lj/charmmfsw/coul/charmmfsh requires atom attribute q
+
+The atom style defined does not have these attributes.
+
+E: Pair inner cutoff >= Pair outer cutoff
+
+The specified cutoffs for the pair style are inconsistent.
+
+*/
+
diff --git a/src/USER-MISC/pair_momb.cpp b/src/USER-MISC/pair_momb.cpp
index 8e7111bbec179dc533b51c30c6aa3eac5a04fc2b..b7337c17a8025feb08eefd258d7ad6517d4c2623 100644
--- a/src/USER-MISC/pair_momb.cpp
+++ b/src/USER-MISC/pair_momb.cpp
@@ -9,8 +9,11 @@
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
- -------------------------------------------------------------------------
- Contributed by Kristen Fichthorn, Tonnam Balankura, Ya Zhou @ Penn State University
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing authors: Kristen Fichthorn, Tonnam Balankura,
+ Ya Zhou (Penn State University)
------------------------------------------------------------------------- */
#include <math.h>