11 Mar 2013 version
+14 Mar 2013 version
Version info:
diff --git a/doc/Manual.txt b/doc/Manual.txt
index 696ac69109..476bd69ac2 100644
--- a/doc/Manual.txt
+++ b/doc/Manual.txt
@@ -18,7 +18,7 @@
LAMMPS-ICMS Documentation :c,h3
-11 Mar 2013 version :c,h4
+14 Mar 2013 version :c,h4
Version info: :h4
diff --git a/doc/Section_howto.html b/doc/Section_howto.html
index b03e117e76..7d8677f8e2 100644
--- a/doc/Section_howto.html
+++ b/doc/Section_howto.html
@@ -445,41 +445,34 @@ set to 0.0, it corresponds to the original 1983 TIP3P model
(Jorgensen).
O mass = 15.9994
-H mass = 1.008
-
-O charge = -0.834
-H charge = 0.417
-
-LJ epsilon of OO = 0.1521
+H mass = 1.008
+O charge = -0.834
+H charge = 0.417
+LJ epsilon of OO = 0.1521
LJ sigma of OO = 3.1507
LJ epsilon of HH = 0.0460
LJ sigma of HH = 0.4000
LJ epsilon of OH = 0.0836
-LJ sigma of OH = 1.7753
-
-K of OH bond = 450
-r0 of OH bond = 0.9572
-
-K of HOH angle = 55
+LJ sigma of OH = 1.7753
+K of OH bond = 450
+r0 of OH bond = 0.9572
+K of HOH angle = 55
theta of HOH angle = 104.52
These are the parameters to use for TIP3P with a long-range Coulombic
-solver (Ewald or PPPM in LAMMPS), see (Price) for details:
+solver (e.g. Ewald or PPPM in LAMMPS), see (Price) for
+details:
O mass = 15.9994
-H mass = 1.008
-
-O charge = -0.830
-H charge = 0.415
-
-LJ epsilon of OO = 0.102
+H mass = 1.008
+O charge = -0.830
+H charge = 0.415
+LJ epsilon of OO = 0.102
LJ sigma of OO = 3.188
-LJ epsilon, sigma of OH, HH = 0.0
-
-K of OH bond = 450
-r0 of OH bond = 0.9572
-
-K of HOH angle = 55
+LJ epsilon, sigma of OH, HH = 0.0
+K of OH bond = 450
+r0 of OH bond = 0.9572
+K of HOH angle = 55
theta of HOH angle = 104.52
Wikipedia also has a nice article on water
@@ -496,15 +489,18 @@ This site M is located at a fixed distance away from the oxygen along
the bisector of the HOH bond angle. A bond style of harmonic and an
angle style of harmonic or charmm should also be used.
-A TIP4P model is run with LAMMPS using two commands:
+
A TIP4P model is run with LAMMPS using either this command
+for a cutoff model:
+
+pair_style lj/cut/tip4p/cut
+
+or these two commands for a long-range model:
-Note that only a TIP4P model with long-range Coulombics is currently
-implemented. A cutoff version may be added in the future. for both
-models, the bond lengths and bond angles should be held fixed using
-the fix shake command.
+
For both models, the bond lengths and bond angles should be held fixed
+using the fix shake command.
These are the additional parameters (in real units) to set for O and H
atoms and the water molecule to run a rigid TIP4P model with a cutoff
@@ -513,35 +509,58 @@ the pair_style command, not as part of the pair
coefficients.
O mass = 15.9994
-H mass = 1.008
-
-O charge = -1.040
-H charge = 0.520
-
-r0 of OH bond = 0.9572
+H mass = 1.008
+O charge = -1.040
+H charge = 0.520
+r0 of OH bond = 0.9572
theta of HOH angle = 104.52
-
-OM distance = 0.15
-
-LJ epsilon of O-O = 0.1550
+OM distance = 0.15
+LJ epsilon of O-O = 0.1550
LJ sigma of O-O = 3.1536
-LJ epsilon, sigma of OH, HH = 0.0
+LJ epsilon, sigma of OH, HH = 0.0
+Coulombic cutoff = 8.5
-These are the parameters to use for TIP4P with a long-range Coulombic
-solver (Ewald or PPPM in LAMMPS):
+
For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
+http://dx.doi.org/10.1063/1.1931662) these values can be used:
O mass = 15.9994
-H mass = 1.008
+H mass = 1.008
+O charge = -1.1794
+H charge = 0.5897
+r0 of OH bond = 0.9572
+theta of HOH angle = 104.52
+OM distance = 0.1577
+LJ epsilon of O-O = 0.21084
+LJ sigma of O-O = 3.1668
+LJ epsilon, sigma of OH, HH = 0.0
+Coulombic cutoff = 8.5
-O charge = -1.0484
-H charge = 0.5242
+
For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
+http://dx.doi.org/10.1063/1.2121687), these values can be used:
-r0 of OH bond = 0.9572
-theta of HOH angle = 104.52
+
O mass = 15.9994
+H mass = 1.008
+O charge = -1.1128
+H charge = 0.5564
+r0 of OH bond = 0.9572
+theta of HOH angle = 104.52
+OM distance = 0.1546
+LJ epsilon of O-O = 0.1852
+LJ sigma of O-O = 3.1589
+LJ epsilon, sigma of OH, HH = 0.0
+Coulombic cutoff = 8.5
-OM distance = 0.1250
+
These are the parameters to use for TIP4P with a long-range Coulombic
+solver (e.g. Ewald or PPPM in LAMMPS):
-LJ epsilon of O-O = 0.16275
+
O mass = 15.9994
+H mass = 1.008
+O charge = -1.0484
+H charge = 0.5242
+r0 of OH bond = 0.9572
+theta of HOH angle = 104.52
+OM distance = 0.1250
+LJ epsilon of O-O = 0.16275
LJ sigma of O-O = 3.16435
LJ epsilon, sigma of OH, HH = 0.0
@@ -574,16 +593,13 @@ used.
atoms and the water molecule to run a rigid SPC model.
O mass = 15.9994
-H mass = 1.008
-
-O charge = -0.820
-H charge = 0.410
-
-LJ epsilon of OO = 0.1553
+H mass = 1.008
+O charge = -0.820
+H charge = 0.410
+LJ epsilon of OO = 0.1553
LJ sigma of OO = 3.166
-LJ epsilon, sigma of OH, HH = 0.0
-
-r0 of OH bond = 1.0
+LJ epsilon, sigma of OH, HH = 0.0
+r0 of OH bond = 1.0
theta of HOH angle = 109.47
Note that as originally proposed, the SPC model was run with a 9
diff --git a/doc/Section_howto.txt b/doc/Section_howto.txt
index 9a1d99e819..f8254f90b9 100644
--- a/doc/Section_howto.txt
+++ b/doc/Section_howto.txt
@@ -440,40 +440,33 @@ set to 0.0, it corresponds to the original 1983 TIP3P model
"(Jorgensen)"_#Jorgensen.
O mass = 15.9994
-H mass = 1.008 :all(b),p
-
+H mass = 1.008
O charge = -0.834
-H charge = 0.417 :all(b),p
-
+H charge = 0.417
LJ epsilon of OO = 0.1521
LJ sigma of OO = 3.1507
LJ epsilon of HH = 0.0460
LJ sigma of HH = 0.4000
LJ epsilon of OH = 0.0836
-LJ sigma of OH = 1.7753 :all(b),p
-
+LJ sigma of OH = 1.7753
K of OH bond = 450
-r0 of OH bond = 0.9572 :all(b),p
-
+r0 of OH bond = 0.9572
K of HOH angle = 55
theta of HOH angle = 104.52 :all(b),p
These are the parameters to use for TIP3P with a long-range Coulombic
-solver (Ewald or PPPM in LAMMPS), see "(Price)"_#Price for details:
+solver (e.g. Ewald or PPPM in LAMMPS), see "(Price)"_#Price for
+details:
O mass = 15.9994
-H mass = 1.008 :all(b),p
-
+H mass = 1.008
O charge = -0.830
-H charge = 0.415 :all(b),p
-
+H charge = 0.415
LJ epsilon of OO = 0.102
LJ sigma of OO = 3.188
-LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
-
+LJ epsilon, sigma of OH, HH = 0.0
K of OH bond = 450
-r0 of OH bond = 0.9572 :all(b),p
-
+r0 of OH bond = 0.9572
K of HOH angle = 55
theta of HOH angle = 104.52 :all(b),p
@@ -491,15 +484,18 @@ This site M is located at a fixed distance away from the oxygen along
the bisector of the HOH bond angle. A bond style of {harmonic} and an
angle style of {harmonic} or {charmm} should also be used.
-A TIP4P model is run with LAMMPS using two commands:
+A TIP4P model is run with LAMMPS using either this command
+for a cutoff model:
+
+"pair_style lj/cut/tip4p/cut"_pair_lj.html
+
+or these two commands for a long-range model:
"pair_style lj/cut/tip4p/long"_pair_lj.html
"kspace_style pppm/tip4p"_kspace_style.html :ul
-Note that only a TIP4P model with long-range Coulombics is currently
-implemented. A cutoff version may be added in the future. for both
-models, the bond lengths and bond angles should be held fixed using
-the "fix shake"_fix_shake.html command.
+For both models, the bond lengths and bond angles should be held fixed
+using the "fix shake"_fix_shake.html command.
These are the additional parameters (in real units) to set for O and H
atoms and the water molecule to run a rigid TIP4P model with a cutoff
@@ -508,34 +504,57 @@ the "pair_style"_pair_style.html command, not as part of the pair
coefficients.
O mass = 15.9994
-H mass = 1.008 :all(b),p
-
+H mass = 1.008
O charge = -1.040
-H charge = 0.520 :all(b),p
-
+H charge = 0.520
r0 of OH bond = 0.9572
-theta of HOH angle = 104.52 :all(b),p
-
-OM distance = 0.15 :all(b),p
-
+theta of HOH angle = 104.52
+OM distance = 0.15
LJ epsilon of O-O = 0.1550
LJ sigma of O-O = 3.1536
-LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
+LJ epsilon, sigma of OH, HH = 0.0
+Coulombic cutoff = 8.5 :all(b),p
-These are the parameters to use for TIP4P with a long-range Coulombic
-solver (Ewald or PPPM in LAMMPS):
+For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
+http://dx.doi.org/10.1063/1.1931662) these values can be used:
O mass = 15.9994
-H mass = 1.008 :all(b),p
-
-O charge = -1.0484
-H charge = 0.5242 :all(b),p
-
+H mass = 1.008
+O charge = -1.1794
+H charge = 0.5897
r0 of OH bond = 0.9572
-theta of HOH angle = 104.52 :all(b),p
+theta of HOH angle = 104.52
+OM distance = 0.1577
+LJ epsilon of O-O = 0.21084
+LJ sigma of O-O = 3.1668
+LJ epsilon, sigma of OH, HH = 0.0
+Coulombic cutoff = 8.5 :all(b),p
-OM distance = 0.1250 :all(b),p
+For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
+http://dx.doi.org/10.1063/1.2121687), these values can be used:
+O mass = 15.9994
+H mass = 1.008
+O charge = -1.1128
+H charge = 0.5564
+r0 of OH bond = 0.9572
+theta of HOH angle = 104.52
+OM distance = 0.1546
+LJ epsilon of O-O = 0.1852
+LJ sigma of O-O = 3.1589
+LJ epsilon, sigma of OH, HH = 0.0
+Coulombic cutoff = 8.5 :all(b),p
+
+These are the parameters to use for TIP4P with a long-range Coulombic
+solver (e.g. Ewald or PPPM in LAMMPS):
+
+O mass = 15.9994
+H mass = 1.008
+O charge = -1.0484
+H charge = 0.5242
+r0 of OH bond = 0.9572
+theta of HOH angle = 104.52
+OM distance = 0.1250
LJ epsilon of O-O = 0.16275
LJ sigma of O-O = 3.16435
LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
@@ -569,15 +588,12 @@ These are the additional parameters (in real units) to set for O and H
atoms and the water molecule to run a rigid SPC model.
O mass = 15.9994
-H mass = 1.008 :all(b),p
-
+H mass = 1.008
O charge = -0.820
-H charge = 0.410 :all(b),p
-
+H charge = 0.410
LJ epsilon of OO = 0.1553
LJ sigma of OO = 3.166
-LJ epsilon, sigma of OH, HH = 0.0 :all(b),p
-
+LJ epsilon, sigma of OH, HH = 0.0
r0 of OH bond = 1.0
theta of HOH angle = 109.47 :all(b),p
diff --git a/doc/compute_temp_asphere.html b/doc/compute_temp_asphere.html
index 59b32f6951..7c70a62924 100644
--- a/doc/compute_temp_asphere.html
+++ b/doc/compute_temp_asphere.html
@@ -66,6 +66,11 @@ should have 3 dof instead of 6 in 3d (or 2 instead of 3 in 2d). A
uniaxial aspherical particle has two of its three shape parameters the
same. If it does not rotate around the axis perpendicular to its
circular cross section, then it should have 5 dof instead of 6 in 3d.
+The latter is the case for uniaxial ellipsoids in a GayBerne
+model since there is no induced torque around the
+optical axis. It will also be the case for biaxial ellipsoids when
+exactly two of the semiaxes have the same length and the corresponding
+relative well depths are equal.
The translational kinetic energy is computed the same as is described
by the compute temp command. The rotational
diff --git a/doc/compute_temp_asphere.txt b/doc/compute_temp_asphere.txt
index 019c1f72c7..ff4a90c610 100644
--- a/doc/compute_temp_asphere.txt
+++ b/doc/compute_temp_asphere.txt
@@ -58,6 +58,11 @@ should have 3 dof instead of 6 in 3d (or 2 instead of 3 in 2d). A
uniaxial aspherical particle has two of its three shape parameters the
same. If it does not rotate around the axis perpendicular to its
circular cross section, then it should have 5 dof instead of 6 in 3d.
+The latter is the case for uniaxial ellipsoids in a "GayBerne
+model"_pair_gayberne.html since there is no induced torque around the
+optical axis. It will also be the case for biaxial ellipsoids when
+exactly two of the semiaxes have the same length and the corresponding
+relative well depths are equal.
The translational kinetic energy is computed the same as is described
by the "compute temp"_compute_temp.html command. The rotational
diff --git a/doc/pair_buck_long.html b/doc/pair_buck_long.html
index fa8ee5769a..2316d7d534 100644
--- a/doc/pair_buck_long.html
+++ b/doc/pair_buck_long.html
@@ -140,9 +140,9 @@ interaction, assuming flag_buck is cut.
shift option for the energy of the Buckingham portion of the pair
interaction.
-This pair style does not support the pair_modify
-table option since a tabulation capability has not yet been added to
-this potential.
+
This pair style supports the pair_modify table and
+table/disp options since they can tabulate the short-range portion of
+the long-range Coulombic and dispersion interactions.
This pair style write its information to binary restart
files, so pair_style and pair_coeff commands do not need
diff --git a/doc/pair_buck_long.txt b/doc/pair_buck_long.txt
index 41ee33d8e4..b8e96214c2 100644
--- a/doc/pair_buck_long.txt
+++ b/doc/pair_buck_long.txt
@@ -131,9 +131,9 @@ This pair style does not support the "pair_modify"_pair_modify.html
shift option for the energy of the Buckingham portion of the pair
interaction.
-This pair style does not support the "pair_modify"_pair_modify.html
-table option since a tabulation capability has not yet been added to
-this potential.
+This pair style supports the "pair_modify"_pair_modify.html table and
+table/disp options since they can tabulate the short-range portion of
+the long-range Coulombic and dispersion interactions.
This pair style write its information to "binary restart
files"_restart.html, so pair_style and pair_coeff commands do not need
diff --git a/doc/pair_lj.html b/doc/pair_lj.html
index 0076b12c29..88da42e1db 100644
--- a/doc/pair_lj.html
+++ b/doc/pair_lj.html
@@ -136,6 +136,11 @@ pair_style lj/cut/coul/msm 10.0 8.0
pair_coeff * * 100.0 3.0
pair_coeff 1 1 100.0 3.5 9.0
+
pair_style lj/cut/tip4p/cut 1 2 7 8 0.15 12.0
+pair_style lj/cut/tip4p/cut 1 2 7 8 0.15 12.0 10.0
+pair_coeff * * 100.0 3.0
+pair_coeff 1 1 100.0 3.5 9.0
+
pair_style lj/cut/tip4p/long 1 2 7 8 0.15 12.0
pair_style lj/cut/tip4p/long 1 2 7 8 0.15 12.0 10.0
pair_coeff * * 100.0 3.0
@@ -195,14 +200,15 @@ specified for this style means that pairwise interactions within this
distance are computed directly; interactions outside that distance are
computed in reciprocal space.
-Style lj/cut/tip4p/long implements the TIP4P water model of
-(Jorgensen), which introduces a massless site located a
-short distance away from the oxygen atom along the bisector of the HOH
-angle. The atomic types of the oxygen and hydrogen atoms, the bond
-and angle types for OH and HOH interactions, and the distance to the
-massless charge site are specified as pair_style arguments. Style
-lj/cut/tip4p/cut is similar to lj/cut/tip4p/long, only it uses a
-cutoff for Coulomb interactions.
+
Styles lj/cut/tip4p/cut and lj/cut/tip4p/long implement the TIP4P
+water model of (Jorgensen), which introduces a massless
+site located a short distance away from the oxygen atom along the
+bisector of the HOH angle. The atomic types of the oxygen and
+hydrogen atoms, the bond and angle types for OH and HOH interactions,
+and the distance to the massless charge site are specified as
+pair_style arguments. Style lj/cut/tip4p/cut uses a cutoff for
+Coulomb interactions; style lj/cut/tip4p/long is for use with a
+long-range Coulombic solver (Ewald or PPPM).
IMPORTANT NOTE: For each TIP4P water molecule in your system, the atom
IDs for the O and 2 H atoms must be consecutive, with the O atom
@@ -211,14 +217,15 @@ with each O atom. For example, if the atom ID of an O atom in a TIP4P
water molecule is 500, then its 2 H atoms must have IDs 501 and 502.
See the howto section for more
-information on how to use the TIP4P pair style. Note that the
-neighobr list cutoff for Coulomb interactions is effectively extended
-by a distance 2*qdist when using the TIP4P pair style, to account for
-the offset distance of the fictitious charges on O atoms in water
-molecules. Thus it is typically best in an efficiency sense to use a
-LJ cutoff >= Coulomb cutoff + 2*qdist, to shrink the size of the
-neighbor list. This leads to slightly larger cost for the long-range
-calculation, so you can test the trade-off for your model.
+information on how to use the TIP4P pair styles and lists of
+parameters to set. Note that the neighobr list cutoff for Coulomb
+interactions is effectively extended by a distance 2*qdist when using
+the TIP4P pair style, to account for the offset distance of the
+fictitious charges on O atoms in water molecules. Thus it is
+typically best in an efficiency sense to use a LJ cutoff >= Coulomb
+cutoff + 2*qdist, to shrink the size of the neighbor list. This leads
+to slightly larger cost for the long-range calculation, so you can
+test the trade-off for your model.
For all of the lj/cut pair styles, the following coefficients must
be defined for each pair of atoms types via the
@@ -244,10 +251,11 @@ are specified, they are used as the LJ and Coulombic cutoffs for this
type pair. You cannot specify 2 cutoffs for style lj/cut, since it
has no Coulombic terms.
-For lj/cut/coul/long and lj/cut/coul/msm and lj/cut/tip4p/long
-only the LJ cutoff can be specified since a Coulombic cutoff cannot be
-specified for an individual I,J type pair. All type pairs use the
-same global Coulombic cutoff specified in the pair_style command.
+
For lj/cut/coul/long and lj/cut/coul/msm and lj/cut/tip4p/cut
+and lj/cut/tip4p/long only the LJ cutoff can be specified since a
+Coulombic cutoff cannot be specified for an individual I,J type pair.
+All type pairs use the same global Coulombic cutoff specified in the
+pair_style command.
diff --git a/doc/pair_lj.txt b/doc/pair_lj.txt
index 3b76d70129..418ed9415c 100644
--- a/doc/pair_lj.txt
+++ b/doc/pair_lj.txt
@@ -105,6 +105,11 @@ pair_style lj/cut/coul/msm 10.0 8.0
pair_coeff * * 100.0 3.0
pair_coeff 1 1 100.0 3.5 9.0 :pre
+pair_style lj/cut/tip4p/cut 1 2 7 8 0.15 12.0
+pair_style lj/cut/tip4p/cut 1 2 7 8 0.15 12.0 10.0
+pair_coeff * * 100.0 3.0
+pair_coeff 1 1 100.0 3.5 9.0 :pre
+
pair_style lj/cut/tip4p/long 1 2 7 8 0.15 12.0
pair_style lj/cut/tip4p/long 1 2 7 8 0.15 12.0 10.0
pair_coeff * * 100.0 3.0
@@ -164,14 +169,15 @@ specified for this style means that pairwise interactions within this
distance are computed directly; interactions outside that distance are
computed in reciprocal space.
-Style {lj/cut/tip4p/long} implements the TIP4P water model of
-"(Jorgensen)"_#Jorgensen, which introduces a massless site located a
-short distance away from the oxygen atom along the bisector of the HOH
-angle. The atomic types of the oxygen and hydrogen atoms, the bond
-and angle types for OH and HOH interactions, and the distance to the
-massless charge site are specified as pair_style arguments. Style
-{lj/cut/tip4p/cut} is similar to {lj/cut/tip4p/long}, only it uses a
-cutoff for Coulomb interactions.
+Styles {lj/cut/tip4p/cut} and {lj/cut/tip4p/long} implement the TIP4P
+water model of "(Jorgensen)"_#Jorgensen, which introduces a massless
+site located a short distance away from the oxygen atom along the
+bisector of the HOH angle. The atomic types of the oxygen and
+hydrogen atoms, the bond and angle types for OH and HOH interactions,
+and the distance to the massless charge site are specified as
+pair_style arguments. Style {lj/cut/tip4p/cut} uses a cutoff for
+Coulomb interactions; style {lj/cut/tip4p/long} is for use with a
+long-range Coulombic solver (Ewald or PPPM).
IMPORTANT NOTE: For each TIP4P water molecule in your system, the atom
IDs for the O and 2 H atoms must be consecutive, with the O atom
@@ -180,14 +186,15 @@ with each O atom. For example, if the atom ID of an O atom in a TIP4P
water molecule is 500, then its 2 H atoms must have IDs 501 and 502.
See the "howto section"_Section_howto.html#howto_8 for more
-information on how to use the TIP4P pair style. Note that the
-neighobr list cutoff for Coulomb interactions is effectively extended
-by a distance 2*qdist when using the TIP4P pair style, to account for
-the offset distance of the fictitious charges on O atoms in water
-molecules. Thus it is typically best in an efficiency sense to use a
-LJ cutoff >= Coulomb cutoff + 2*qdist, to shrink the size of the
-neighbor list. This leads to slightly larger cost for the long-range
-calculation, so you can test the trade-off for your model.
+information on how to use the TIP4P pair styles and lists of
+parameters to set. Note that the neighobr list cutoff for Coulomb
+interactions is effectively extended by a distance 2*qdist when using
+the TIP4P pair style, to account for the offset distance of the
+fictitious charges on O atoms in water molecules. Thus it is
+typically best in an efficiency sense to use a LJ cutoff >= Coulomb
+cutoff + 2*qdist, to shrink the size of the neighbor list. This leads
+to slightly larger cost for the long-range calculation, so you can
+test the trade-off for your model.
For all of the {lj/cut} pair styles, the following coefficients must
be defined for each pair of atoms types via the
@@ -213,10 +220,11 @@ are specified, they are used as the LJ and Coulombic cutoffs for this
type pair. You cannot specify 2 cutoffs for style {lj/cut}, since it
has no Coulombic terms.
-For {lj/cut/coul/long} and {lj/cut/coul/msm} and {lj/cut/tip4p/long}
-only the LJ cutoff can be specified since a Coulombic cutoff cannot be
-specified for an individual I,J type pair. All type pairs use the
-same global Coulombic cutoff specified in the pair_style command.
+For {lj/cut/coul/long} and {lj/cut/coul/msm} and {lj/cut/tip4p/cut}
+and {lj/cut/tip4p/long} only the LJ cutoff can be specified since a
+Coulombic cutoff cannot be specified for an individual I,J type pair.
+All type pairs use the same global Coulombic cutoff specified in the
+pair_style command.
:line
diff --git a/src/KSPACE/pair_buck_long_coul_long.cpp b/src/KSPACE/pair_buck_long_coul_long.cpp
index b7435dcc63..859dac8db0 100644
--- a/src/KSPACE/pair_buck_long_coul_long.cpp
+++ b/src/KSPACE/pair_buck_long_coul_long.cpp
@@ -51,6 +51,7 @@ PairBuckLongCoulLong::PairBuckLongCoulLong(LAMMPS *lmp) : Pair(lmp)
dispersionflag = ewaldflag = pppmflag = 1;
respa_enable = 1;
ftable = NULL;
+ fdisptable = NULL;
}
/* ----------------------------------------------------------------------
@@ -94,7 +95,7 @@ void PairBuckLongCoulLong::settings(int narg, char **arg)
if (!((ewald_order^ewald_off) & (1<<1)))
error->all(FLERR,"Coulomb cut not supported in pair_style buck/long/coul/coul");
cut_buck_global = force->numeric(*(arg++));
- if (*arg && ((ewald_order & 0x42) == 0x42))
+ if (narg == 4 && ((ewald_order & 0x42) == 0x42))
error->all(FLERR,"Only one cutoff allowed when requesting all long");
if (narg == 4) cut_coul = force->numeric(*arg);
else cut_coul = cut_buck_global;
@@ -282,10 +283,10 @@ void PairBuckLongCoulLong::init_style()
error->all(FLERR,"Pair style requires a KSpace style");
if (ewald_order&(1<<1)) g_ewald = force->kspace->g_ewald;
if (ewald_order&(1<<6)) g_ewald_6 = force->kspace->g_ewald_6;
-
// setup force tables
if (ncoultablebits) init_tables(cut_coul,cut_respa);
+ if (ndisptablebits) init_tables_disp(cut_buck_global);
}
/* ----------------------------------------------------------------------
@@ -309,7 +310,8 @@ double PairBuckLongCoulLong::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
- cut_buck[i][j] = cut_buck_read[i][j];
+ if (ewald_order&(1<<6)) cut_buck[i][j] = cut_buck_global;
+ else cut_buck[i][j] = cut_buck_read[i][j];
buck_a[i][j] = buck_a_read[i][j];
buck_c[i][j] = buck_c_read[i][j];
buck_rho[i][j] = buck_rho_read[i][j];
@@ -441,6 +443,7 @@ void PairBuckLongCoulLong::read_restart_settings(FILE *fp)
void PairBuckLongCoulLong::compute(int eflag, int vflag)
{
+
double evdwl,ecoul,fpair;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@@ -528,19 +531,36 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
register double rn = r2inv*r2inv*r2inv,
expr = exp(-r*rhoinvi[typej]);
if (order6) { // long-range
- register double x2 = g2*rsq, a2 = 1.0/x2;
- x2 = a2*exp(-x2)*buckci[typej];
- if (ni == 0) {
- force_buck =
- r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
- if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
+ if (!ndisptablebits || rsq <= tabinnerdispsq) {
+ register double x2 = g2*rsq, a2 = 1.0/x2;
+ x2 = a2*exp(-x2)*buckci[typej];
+ if (ni == 0) {
+ force_buck =
+ r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
+ if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
+ }
+ else { // special case
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej]-
+ g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
+ if (eflag) evdwl = f*expr*buckai[typej] -
+ g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
+ }
}
- else { // special case
- register double f = special_lj[ni], t = rn*(1.0-f);
- force_buck = f*r*expr*buck1i[typej]-
- g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
- if (eflag) evdwl = f*expr*buckai[typej] -
- g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
+ else { //table real space
+ register union_int_float_t disp_t;
+ disp_t.f = rsq;
+ register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+ register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+ if (ni == 0) {
+ force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej];
+ if (eflag) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
+ }
+ else { //speial case
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej];
+ if (eflag) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
+ }
}
}
else { // cut
@@ -607,15 +627,14 @@ void PairBuckLongCoulLong::compute_inner()
double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi;
vector xi, d;
- ineighn = (ineigh = listinner->ilist)+listinner->inum;
-
+ ineighn = (ineigh = list->ilist) + list->inum;
for (; ineighfirstneigh[i])+listinner->numneigh[i];
+ jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
for (; jneighilist)+listmiddle->inum;
+ ineighn = (ineigh = list->ilist)+list->inum;
for (; ineighfirstneigh[i])+listmiddle->numneigh[i];
+ jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
for (; jneighilist)+listouter->inum;
+ ineighn = (ineigh = list->ilist)+list->inum;
for (; ineighfirstneigh[i])+listouter->numneigh[i];
+ jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
for (; jneighcut_in_off_sq)&&(rsq cut_in_off_sq)) {
register double rsw = (r-cut_in_off)/cut_in_diff;
- frespa = rsw*rsw*(3.0-2.0*rsw);
+ frespa = 1-rsw*rsw*(3.0-2.0*rsw);
}
if (order1 && (rsq < cut_coulsq)) { // coulombic
@@ -835,19 +858,20 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
if (ni == 0) {
s *= g_ewald*exp(-x*x);
- force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
+ force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
if (eflag) ecoul = t;
}
else { // correct for special
- r = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
- force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r;
- if (eflag) ecoul = t-r;
+ register double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
+ force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-ri-respa_coul;
+ if (eflag) ecoul = t-ri;
}
} // table real space
else {
- if (respa_flag) respa_coul = ni == 0 ? // correct for respa
- frespa*qri*q[j]/r :
- frespa*qri*q[j]/r*special_coul[ni];
+ if (respa_flag) {
+ register double s = qri*q[j];
+ respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
+ }
register union_int_float_t t;
t.f = rsq;
register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
@@ -859,7 +883,10 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
else { // correct for special
t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
- if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
+ if (eflag) {
+ t.f = (1.0-special_coul[ni])*(ptable[k]+f*dptable[k]);
+ ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
+ }
}
}
}
@@ -872,30 +899,48 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) :
frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
if (order6) { // long-range form
- register double x2 = g2*rsq, a2 = 1.0/x2;
- x2 = a2*exp(-x2)*buckci[typej];
- if (ni == 0) {
- force_buck =
- r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
- if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
+ if (!ndisptablebits || rsq <= tabinnerdispsq) {
+ register double x2 = g2*rsq, a2 = 1.0/x2;
+ x2 = a2*exp(-x2)*buckci[typej];
+ if (ni == 0) {
+ force_buck =
+ r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq-respa_buck;
+ if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
+ }
+ else { // correct for special
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej]-
+ g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]-respa_buck;
+ if (eflag) evdwl = f*expr*buckai[typej] -
+ g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
+ }
}
- else { // correct for special
- register double f = special_lj[ni], t = rn*(1.0-f);
- force_buck = f*r*expr*buck1i[typej]-
- g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
- if (eflag) evdwl = f*expr*buckai[typej] -
- g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
+ else { // table real space
+ register union_int_float_t disp_t;
+ disp_t.f = rsq;
+ register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+ register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+ register double rn = r2inv*r2inv*r2inv;
+ if (ni == 0) {
+ force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]-respa_buck;
+ if (eflag) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
+ }
+ else { //special case
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]+t*buck2i[typej]-respa_buck;
+ if (eflag) evdwl = f*expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
+ }
}
}
else { // cut form
if (ni == 0) {
- force_buck = r*expr*buck1i[typej]-rn*buck2i[typej];
+ force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]-respa_buck;
if (eflag)
evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej];
}
else { // correct for special
register double f = special_lj[ni];
- force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]);
+ force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej])-respa_buck;
if (eflag)
evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
}
@@ -904,22 +949,24 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
else force_buck = respa_buck = evdwl = 0.0;
fpair = (force_coul+force_buck)*r2inv;
- frespa = fpair-(respa_coul+respa_buck)*r2inv;
if (newton_pair || j < nlocal) {
register double *fj = f0+(j+(j<<1)), f;
- fi[0] += f = d[0]*frespa; fj[0] -= f;
- fi[1] += f = d[1]*frespa; fj[1] -= f;
- fi[2] += f = d[2]*frespa; fj[2] -= f;
+ fi[0] += f = d[0]*fpair; fj[0] -= f;
+ fi[1] += f = d[1]*fpair; fj[1] -= f;
+ fi[2] += f = d[2]*fpair; fj[2] -= f;
}
else {
- fi[0] += d[0]*frespa;
- fi[1] += d[1]*frespa;
- fi[2] += d[2]*frespa;
+ fi[0] += d[0]*fpair;
+ fi[1] += d[1]*fpair;
+ fi[2] += d[2]*fpair;
}
- if (evflag) ev_tally(i,j,nlocal,newton_pair,
- evdwl,ecoul,fpair,d[0],d[1],d[2]);
+ if (evflag) {
+ fvirial = (force_coul + force_buck + respa_coul + respa_buck)*r2inv;
+ ev_tally(i,j,nlocal,newton_pair,
+ evdwl,ecoul,fvirial,d[0],d[1],d[2]);
+ }
}
}
}
diff --git a/src/KSPACE/pair_buck_long_coul_long.h b/src/KSPACE/pair_buck_long_coul_long.h
index 28759bef5a..c1bcebf83c 100644
--- a/src/KSPACE/pair_buck_long_coul_long.h
+++ b/src/KSPACE/pair_buck_long_coul_long.h
@@ -45,9 +45,9 @@ class PairBuckLongCoulLong : public Pair {
double single(int, int, int, int, double, double, double, double &);
void *extract(const char *, int &);
- void compute_inner();
- void compute_middle();
- void compute_outer(int, int);
+ virtual void compute_inner();
+ virtual void compute_middle();
+ virtual void compute_outer(int, int);
protected:
double cut_buck_global;
diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp
index 58cecb039f..62a22bf7a7 100644
--- a/src/KSPACE/pppm_disp.cpp
+++ b/src/KSPACE/pppm_disp.cpp
@@ -3355,11 +3355,13 @@ void PPPMDisp::compute_gf_6()
sqk = pow(qx,2.0) + pow(qy,2.0) + pow(qz,2.0);
- denominator = gf_denom(snx2,sny2,snz2, gf_b_6, order_6);
- rtsqk = sqrt(sqk);
- term = (1-2*sqk*inv2ew*inv2ew)*sx*sy*sz +
- 2*sqk*rtsqk*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtsqk*inv2ew);
- greensfn_6[n++] = numerator*term*wx*wy*wz/denominator;
+ if (sqk != 0.0) {
+ denominator = gf_denom(snx2,sny2,snz2, gf_b_6, order_6);
+ rtsqk = sqrt(sqk);
+ term = (1-2*sqk*inv2ew*inv2ew)*sx*sy*sz +
+ 2*sqk*rtsqk*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtsqk*inv2ew);
+ greensfn_6[n++] = numerator*term*wx*wy*wz/denominator;
+ } else greensfn_6[n++] = 0.0;
}
}
}
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index 0522b70146..997c791f5a 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -162,7 +162,10 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
gcmc_nmax = 0;
local_gas_list = NULL;
- model_atom = NULL;
+ atom_coord = NULL;
+ model_atom = NULL;
+ model_atom_buf = NULL;
+
}
/* ----------------------------------------------------------------------
@@ -204,6 +207,7 @@ void FixGCMC::options(int narg, char **arg)
FixGCMC::~FixGCMC()
{
+ if (regionflag) delete [] idregion;
delete random_equal;
delete random_unequal;
memory->destroy(local_gas_list);
diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
index 4bd52d293e..72f809aec7 100644
--- a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
+++ b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
@@ -26,10 +26,12 @@
#include "angle.h"
#include "bond.h"
#include "comm.h"
+#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
+using namespace MathConst;
/* ---------------------------------------------------------------------- */
@@ -408,7 +410,7 @@ void PairLJCutTIP4PCut::allocate()
void PairLJCutTIP4PCut::settings(int narg, char **arg)
{
- if (narg != 7) error->all(FLERR,"Illegal pair_style command");
+ if (narg < 6 || narg > 7) error->all(FLERR,"Illegal pair_style command");
typeO = force->inumeric(arg[0]);
typeH = force->inumeric(arg[1]);
@@ -417,7 +419,9 @@ void PairLJCutTIP4PCut::settings(int narg, char **arg)
qdist = force->numeric(arg[4]);
cut_lj_global = force->numeric(arg[5]);
- cut_coul = force->numeric(arg[6]);
+ if (narg == 6) cut_coul = cut_lj_global;
+ else cut_coul = force->numeric(arg[6]);
+
cut_coulsq = cut_coul * cut_coul;
cut_coulsqplus = (cut_coul + 2.0*qdist) * (cut_coul + 2.0*qdist);
@@ -435,7 +439,8 @@ void PairLJCutTIP4PCut::settings(int narg, char **arg)
void PairLJCutTIP4PCut::coeff(int narg, char **arg)
{
- if (narg != 4) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg < 4 || narg > 5)
+ error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@@ -445,12 +450,15 @@ void PairLJCutTIP4PCut::coeff(int narg, char **arg)
double epsilon_one = force->numeric(arg[2]);
double sigma_one = force->numeric(arg[3]);
+ double cut_lj_one = cut_lj_global;
+ if (narg == 5) cut_lj_one = force->numeric(arg[4]);
+
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;
- cut_lj[i][j] = cut_lj_global;
+ cut_lj[i][j] = cut_lj_one;
setflag[i][j] = 1;
count++;
}
@@ -522,6 +530,32 @@ double PairLJCutTIP4PCut::init_one(int i, int j)
lj4[j][i] = lj4[i][j];
offset[j][i] = offset[i][j];
+ // compute I,J contribution to long-range tail correction
+ // count total # of atoms of type I and J via Allreduce
+
+ if (tail_flag) {
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ double count[2],all[2];
+ count[0] = count[1] = 0.0;
+ for (int k = 0; k < nlocal; k++) {
+ if (type[k] == i) count[0] += 1.0;
+ if (type[k] == j) count[1] += 1.0;
+ }
+ MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
+
+ double sig2 = sigma[i][j]*sigma[i][j];
+ double sig6 = sig2*sig2*sig2;
+ double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
+ double rc6 = rc3*rc3;
+ double rc9 = rc3*rc6;
+ etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
+ sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
+ ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
+ sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
+ }
+
// check that LJ epsilon = 0.0 for water H
// set LJ cutoff to 0.0 for any interaction involving water H
// so LJ term isn't calculated in compute()
@@ -598,10 +632,12 @@ void PairLJCutTIP4PCut::write_restart_settings(FILE *fp)
fwrite(&typeB,sizeof(int),1,fp);
fwrite(&typeA,sizeof(int),1,fp);
fwrite(&qdist,sizeof(double),1,fp);
+
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul,sizeof(double),1,fp);
- fwrite(&cut_coulsq,sizeof(double),1,fp);
- fwrite(&cut_coulsqplus,sizeof(double),1,fp);
+ fwrite(&offset_flag,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+ fwrite(&tail_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
@@ -616,10 +652,12 @@ void PairLJCutTIP4PCut::read_restart_settings(FILE *fp)
fread(&typeB,sizeof(int),1,fp);
fread(&typeA,sizeof(int),1,fp);
fread(&qdist,sizeof(double),1,fp);
+
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul,sizeof(double),1,fp);
- fread(&cut_coulsq,sizeof(double),1,fp);
- fread(&cut_coulsqplus,sizeof(double),1,fp);
+ fread(&offset_flag,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ fread(&tail_flag,sizeof(int),1,fp);
}
MPI_Bcast(&typeO,1,MPI_INT,0,world);
@@ -627,10 +665,15 @@ void PairLJCutTIP4PCut::read_restart_settings(FILE *fp)
MPI_Bcast(&typeB,1,MPI_INT,0,world);
MPI_Bcast(&typeA,1,MPI_INT,0,world);
MPI_Bcast(&qdist,1,MPI_DOUBLE,0,world);
+
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut_coulsq,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut_coulsqplus,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(&tail_flag,1,MPI_INT,0,world);
+
+ cut_coulsq = cut_coul * cut_coul;
+ cut_coulsqplus = (cut_coul + 2.0*qdist) * (cut_coul + 2.0*qdist);
}
/* ----------------------------------------------------------------------
diff --git a/src/USER-OMP/msm_omp.h b/src/USER-OMP/msm_omp.h
index a086b7af92..0a41c57ac3 100644
--- a/src/USER-OMP/msm_omp.h
+++ b/src/USER-OMP/msm_omp.h
@@ -1,46 +1,46 @@
-/* -*- 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 KSPACE_CLASS
-
-KSpaceStyle(msm/omp,MSMOMP)
-
-#else
-
-#ifndef LMP_MSM_OMP_H
-#define LMP_MSM_OMP_H
-
-#include "msm.h"
-#include "thr_omp.h"
-
-namespace LAMMPS_NS {
-
- class MSMOMP : public MSM, public ThrOMP {
- public:
- MSMOMP(class LAMMPS *, int, char **);
- virtual ~MSMOMP () {};
-
- protected:
- virtual void direct(int);
- virtual void compute(int,int);
-
- private:
- template void direct_eval(int);
- template void direct_peratom(int);
-
-};
-
-}
-
-#endif
-#endif
+/* -*- 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 KSPACE_CLASS
+
+KSpaceStyle(msm/omp,MSMOMP)
+
+#else
+
+#ifndef LMP_MSM_OMP_H
+#define LMP_MSM_OMP_H
+
+#include "msm.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+ class MSMOMP : public MSM, public ThrOMP {
+ public:
+ MSMOMP(class LAMMPS *, int, char **);
+ virtual ~MSMOMP () {};
+
+ protected:
+ virtual void direct(int);
+ virtual void compute(int,int);
+
+ private:
+ template void direct_eval(int);
+ template void direct_peratom(int);
+
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
index e464472e16..ba31833d7a 100644
--- a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
+++ b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
@@ -3,23 +3,23 @@
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
- This software is distributed under the GNU General Public License.
+ 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: Axel Kohlmeyer (Temple U)
-------------------------------------------------------------------------- */
#include "math.h"
+#include "math_vector.h"
#include "pair_buck_long_coul_long_omp.h"
#include "atom.h"
#include "comm.h"
-#include "math_vector.h"
-#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
+#include "force.h"
#include "suffix.h"
using namespace LAMMPS_NS;
@@ -34,11 +34,11 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
-PairBuckLongCoulLongOMP::PairBuckLongCoulLongOMP(LAMMPS *lmp) :
+PairBuckLongCoulLongOMP::PairBuckLongCoulLongOMP(LAMMPS *lmp) :
PairBuckLongCoulLong(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
- respa_enable = 0;
+ respa_enable = 1;
cut_respa = NULL;
}
@@ -49,6 +49,8 @@ void PairBuckLongCoulLongOMP::compute(int eflag, int vflag)
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 0;
+ const int order1 = ewald_order&(1<<1);
+ const int order6 = ewald_order&(1<<6);
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
@@ -64,18 +66,243 @@ void PairBuckLongCoulLongOMP::compute(int eflag, int vflag)
ThrData *thr = fix->get_thr(tid);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
- if (evflag) {
- if (eflag) {
- if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
- else eval<1,1,0>(ifrom, ito, thr);
+ if (order6) {
+ if (order1) {
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,0,1,1>(ifrom, ito, thr);
+ else eval<1,1,0,0,0,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,0,1,1>(ifrom, ito, thr);
+ else eval<1,0,0,0,0,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,0,1,1>(ifrom, ito, thr);
+ else eval<0,0,0,0,0,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,0,1,1>(ifrom, ito, thr);
+ else eval<1,1,0,1,0,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,0,1,1>(ifrom, ito, thr);
+ else eval<1,0,0,1,0,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,0,1,1>(ifrom, ito, thr);
+ else eval<0,0,0,1,0,1,1>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0,0,1,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,1,1,1>(ifrom, ito, thr);
+ else eval<1,0,0,0,1,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,1,1,1>(ifrom, ito, thr);
+ else eval<0,0,0,0,1,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0,1,1,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,1,1,1>(ifrom, ito, thr);
+ else eval<1,0,0,1,1,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,1,1,1>(ifrom, ito, thr);
+ else eval<0,0,0,1,1,1,1>(ifrom, ito, thr);
+ }
+ }
+ }
} else {
- if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
- else eval<1,0,0>(ifrom, ito, thr);
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,0,0,1>(ifrom, ito, thr);
+ else eval<1,1,0,0,0,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,0,0,1>(ifrom, ito, thr);
+ else eval<1,0,0,0,0,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0,0,0,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,0,0,1>(ifrom, ito, thr);
+ else eval<1,1,0,1,0,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,0,0,1>(ifrom, ito, thr);
+ else eval<1,0,0,1,0,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0,1,0,0,1>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,1,0,1>(ifrom, ito, thr);
+ else eval<1,1,0,0,1,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0,0,1,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,1,0,1>(ifrom, ito, thr);
+ else eval<0,0,0,0,1,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,1,0,1>(ifrom, ito, thr);
+ else eval<1,1,0,1,1,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0,1,1,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,1,0,1>(ifrom, ito, thr);
+ else eval<0,0,0,1,1,0,1>(ifrom, ito, thr);
+ }
+ }
+ }
}
} else {
- if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
- else eval<0,0,0>(ifrom, ito, thr);
- }
+ if (order1) {
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,0,1,0>(ifrom, ito, thr);
+ else eval<1,1,0,0,0,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,0,1,0>(ifrom, ito, thr);
+ else eval<1,0,0,0,0,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,0,1,0>(ifrom, ito, thr);
+ else eval<0,0,0,0,0,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,0,1,0>(ifrom, ito, thr);
+ else eval<1,1,0,1,0,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,0,1,0>(ifrom, ito, thr);
+ else eval<1,0,0,1,0,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,0,1,0>(ifrom, ito, thr);
+ else eval<0,0,0,1,0,1,0>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,1,1,0>(ifrom, ito, thr);
+ else eval<1,1,0,0,1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,1,1,0>(ifrom, ito, thr);
+ else eval<1,0,0,0,1,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,1,1,0>(ifrom, ito, thr);
+ else eval<0,0,0,0,1,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,1,1,0>(ifrom, ito, thr);
+ else eval<1,1,0,1,1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,1,1,0>(ifrom, ito, thr);
+ else eval<1,0,0,1,1,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,1,1,0>(ifrom, ito, thr);
+ else eval<0,0,0,1,1,1,0>(ifrom, ito, thr);
+ }
+ }
+ }
+ } else {
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,0,0,0>(ifrom, ito, thr);
+ else eval<1,1,0,0,0,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,0,0,0>(ifrom, ito, thr);
+ else eval<1,0,0,0,0,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,0,0,0>(ifrom, ito, thr);
+ else eval<0,0,0,0,0,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,0,0,0>(ifrom, ito, thr);
+ else eval<1,1,0,1,0,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,0,0,0>(ifrom, ito, thr);
+ else eval<1,0,0,1,0,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,0,0,0>(ifrom, ito, thr);
+ else eval<0,0,0,1,0,0,0>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,0,1,0,0>(ifrom, ito, thr);
+ else eval<1,1,0,0,1,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,0,1,0,0>(ifrom, ito, thr);
+ else eval<1,0,0,0,1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,0,1,0,0>(ifrom, ito, thr);
+ else eval<0,0,0,0,1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1,1,1,0,0>(ifrom, ito, thr);
+ else eval<1,1,0,1,1,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1,1,1,0,0>(ifrom, ito, thr);
+ else eval<1,0,0,1,1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1,1,1,0,0>(ifrom, ito, thr);
+ else eval<0,0,0,1,1,0,0>(ifrom, ito, thr);
+ }
+ }
+ }
+ }
+ }
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
@@ -83,20 +310,335 @@ void PairBuckLongCoulLongOMP::compute(int eflag, int vflag)
/* ---------------------------------------------------------------------- */
-template
+/* ---------------------------------------------------------------------- */
+
+void PairBuckLongCoulLongOMP::compute_inner()
+{
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(0, 0, nall, 0, 0, thr);
+ eval_inner(ifrom, ito, thr);
+
+ } // end of omp parallel region
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBuckLongCoulLongOMP::compute_middle()
+{
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(0, 0, nall, 0, 0, thr);
+ eval_middle(ifrom, ito, thr);
+
+ } // end of omp parallel region
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBuckLongCoulLongOMP::compute_outer(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+ const int order1 = ewald_order&(1<<1);
+ const int order6 = ewald_order&(1<<6);
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ if (order6) {
+ if (order1) {
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,0,1,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,0,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,0,1,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,0,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,0,1,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,0,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,0,1,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,0,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,0,1,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,0,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,0,1,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,0,1,1>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,1,1,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,1,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,1,1,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,1,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,1,1,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,1,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,1,1,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,1,1,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,1,1,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,1,1,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,1,1,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,1,1,1>(ifrom, ito, thr);
+ }
+ }
+ }
+ } else {
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,0,0,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,0,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,0,0,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,0,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,0,0,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,0,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,0,0,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,0,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,0,0,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,0,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,0,0,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,0,0,1>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,1,0,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,1,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,1,0,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,1,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,1,0,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,1,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,1,0,1>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,1,0,1>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,1,0,1>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,1,0,1>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,1,0,1>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,1,0,1>(ifrom, ito, thr);
+ }
+ }
+ }
+ }
+ } else {
+ if (order1) {
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,0,1,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,0,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,0,1,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,0,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,0,1,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,0,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,0,1,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,0,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,0,1,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,0,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,0,1,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,0,1,0>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,1,1,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,1,1,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,1,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,1,1,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,1,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,1,1,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,1,1,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,1,1,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,1,1,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,1,1,0>(ifrom, ito, thr);
+ }
+ }
+ }
+ } else {
+ if (!ndisptablebits) {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,0,0,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,0,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,0,0,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,0,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,0,0,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,0,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,0,0,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,0,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,0,0,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,0,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,0,0,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,0,0,0>(ifrom, ito, thr);
+ }
+ }
+ } else {
+ if (!ncoultablebits) {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,0,1,0,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,0,1,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,0,1,0,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,0,1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,0,1,0,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,0,1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval_outer<1,1,1,1,1,0,0>(ifrom, ito, thr);
+ else eval_outer<1,1,0,1,1,0,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval_outer<1,0,1,1,1,0,0>(ifrom, ito, thr);
+ else eval_outer<1,0,0,1,1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval_outer<0,0,1,1,1,0,0>(ifrom, ito, thr);
+ else eval_outer<0,0,0,1,1,0,0>(ifrom, ito, thr);
+ }
+ }
+ }
+ }
+ }
+
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+/* ---------------------------------------------------------------------- */
+
+template < const int EVFLAG, const int EFLAG,
+ const int NEWTON_PAIR, const int CTABLE, const int DISPTABLE, const int ORDER1, const int ORDER6 >
void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
{
+
double evdwl,ecoul,fpair;
evdwl = ecoul = 0.0;
const double * const * const x = atom->x;
double * const * const f = thr->get_f();
const double * const q = atom->q;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- double *special_coul = force->special_coul;
- double *special_lj = force->special_lj;
- double qqrd2e = force->qqrd2e;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+ const double * const special_coul = force->special_coul;
+ const double * const special_lj = force->special_lj;
+ const double qqrd2e = force->qqrd2e;
const double *x0 = x[0];
double *f0 = f[0], *fi = f0;
@@ -105,17 +647,17 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
// loop over neighbors of my atoms
- int i, ii, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6);
- int *jneigh, *jneighn, typei, typej, ni;
- double qi, qri, *cutsqi, *cut_bucksqi,
+ int i, ii, j;
+ int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni;
+ double qi = 0.0, qri = 0.0, *cutsqi, *cut_bucksqi,
*buck1i, *buck2i, *buckai, *buckci, *rhoinvi, *offseti;
double r, rsq, r2inv, force_coul, force_buck;
- double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2;
+ double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2;
vector xi, d;
for (ii = iifrom; ii < iito; ++ii) { // loop over my atoms
- i = ilist[ii]; fi = f0+3*i;
- if (order1) qri = (qi = q[i])*qqrd2e; // initialize constants
+ i = ilist[ii]; fi=f0+3*i;
+ if (ORDER1) qri = (qi = q[i])*qqrd2e; // initialize constants
offseti = offset[typei = type[i]];
buck1i = buck1[typei]; buck2i = buck2[typei];
buckai = buck_a[typei]; buckci = buck_c[typei], rhoinvi = rhoinv[typei];
@@ -128,7 +670,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
ni = sbmask(j);
j &= NEIGHMASK;
- { const register double *xj = x0+(j+(j<<1));
+ { register const double *xj = x0+(j+(j<<1));
d[0] = xi[0] - xj[0]; // pair vector
d[1] = xi[1] - xj[1];
d[2] = xi[2] - xj[2]; }
@@ -137,21 +679,23 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
r2inv = 1.0/rsq;
r = sqrt(rsq);
- if (order1 && (rsq < cut_coulsq)) { // coulombic
- if (!ncoultablebits || rsq <= tabinnersq) { // series real space
+ if (ORDER1 && (rsq < cut_coulsq)) { // coulombic
+ if (!CTABLE || rsq <= tabinnersq) { // series real space
register double x = g_ewald*r;
register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
if (ni == 0) {
s *= g_ewald*exp(-x*x);
force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
if (EFLAG) ecoul = t;
- } else { // special case
+ }
+ else { // special case
register double f = s*(1.0-special_coul[ni])/r;
s *= g_ewald*exp(-x*x);
force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
if (EFLAG) ecoul = t-f;
- } // table real space
- } else {
+ }
+ } // table real space
+ else {
register union_int_float_t t;
t.f = rsq;
register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
@@ -166,38 +710,60 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
}
}
- } else force_coul = ecoul = 0.0;
+ }
+ else force_coul = ecoul = 0.0;
if (rsq < cut_bucksqi[typej]) { // buckingham
register double rn = r2inv*r2inv*r2inv,
expr = exp(-r*rhoinvi[typej]);
- if (order6) { // long-range
- register double x2 = g2*rsq, a2 = 1.0/x2;
- x2 = a2*exp(-x2)*buckci[typej];
- if (ni == 0) {
- force_buck =
- r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
- if (EFLAG) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
- } else { // special case
- register double f = special_lj[ni], t = rn*(1.0-f);
- force_buck = f*r*expr*buck1i[typej]-
- g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
- if (EFLAG) evdwl = f*expr*buckai[typej] -
- g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
+ if (ORDER6) { // long-range
+ if (!DISPTABLE || rsq <= tabinnerdispsq) {
+ register double x2 = g2*rsq, a2 = 1.0/x2;
+ x2 = a2*exp(-x2)*buckci[typej];
+ if (ni == 0) {
+ force_buck =
+ r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
+ if (EFLAG) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
+ }
+ else { // special case
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej]-
+ g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
+ if (EFLAG) evdwl = f*expr*buckai[typej] -
+ g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
+ }
}
- } else { // cut
+ else { //table real space
+ register union_int_float_t disp_t;
+ disp_t.f = rsq;
+ register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+ register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+ if (ni == 0) {
+ force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej];
+ if (EFLAG) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
+ }
+ else { //speial case
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej];
+ if (EFLAG) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
+ }
+ }
+ }
+ else { // cut
if (ni == 0) {
force_buck = r*expr*buck1i[typej]-rn*buck2i[typej];
if (EFLAG) evdwl = expr*buckai[typej] -
rn*buckci[typej]-offseti[typej];
- } else { // special case
+ }
+ else { // special case
register double f = special_lj[ni];
force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]);
if (EFLAG)
evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
}
}
- } else force_buck = evdwl = 0.0;
+ }
+ else force_buck = evdwl = 0.0;
fpair = (force_coul+force_buck)*r2inv;
@@ -206,24 +772,404 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
fi[0] += f = d[0]*fpair; fj[0] -= f;
fi[1] += f = d[1]*fpair; fj[1] -= f;
fi[2] += f = d[2]*fpair; fj[2] -= f;
- } else {
+ }
+ else {
fi[0] += d[0]*fpair;
fi[1] += d[1]*fpair;
fi[2] += d[2]*fpair;
}
if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
- evdwl,ecoul,fpair,d[0],d[1],d[2],thr);
+ evdwl,ecoul,fpair,d[0],d[1],d[2],thr);
}
}
}
/* ---------------------------------------------------------------------- */
-double PairBuckLongCoulLongOMP::memory_usage()
+void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr)
{
- double bytes = memory_usage_thr();
- bytes += PairBuckLongCoulLong::memory_usage();
+ double r, rsq, r2inv, force_coul = 0.0, force_buck, fpair;
- return bytes;
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ const double * const q = atom->q;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+ const double * const special_coul = force->special_coul;
+ const double * const special_lj = force->special_lj;
+ const double qqrd2e = force->qqrd2e;
+
+ const double *x0 = x[0];
+ double *f0 = f[0], *fi = 0;
+
+ int *ilist = list->ilist;
+
+ const int newton_pair = force->newton_pair;
+
+ const double cut_out_on = cut_respa[0];
+ const double cut_out_off = cut_respa[1];
+
+ const double cut_out_diff = cut_out_off - cut_out_on;
+ const double cut_out_on_sq = cut_out_on*cut_out_on;
+ const double cut_out_off_sq = cut_out_off*cut_out_off;
+
+
+ int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni;
+ const int order1 = (ewald_order|(ewald_off^-1))&(1<<1);
+ int i, j, ii;
+ double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi;
+ vector xi, d;
+
+ for (ii = iifrom; ii < iito; ++ii) { // loop over my atoms
+ i = ilist[ii]; fi = f0+3*i;
+ if (order1) qri = qqrd2e*q[i];
+ memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
+ cut_bucksqi = cut_bucksq[typei = type[i]];
+ buck1i = buck1[typei]; buck2i = buck2[typei]; rhoinvi = rhoinv[typei];
+ jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+
+ for (; jneigh= cut_out_off_sq) continue;
+ r2inv = 1.0/rsq;
+ r = sqrt(rsq);
+
+ if (order1 && (rsq < cut_coulsq)) // coulombic
+ force_coul = ni == 0 ?
+ qri*q[j]/r : qri*q[j]/r*special_coul[ni];
+
+ if (rsq < cut_bucksqi[typej = type[j]]) { // buckingham
+ register double rn = r2inv*r2inv*r2inv,
+ expr = exp(-r*rhoinvi[typej]);
+ force_buck = ni == 0 ?
+ (r*expr*buck1i[typej]-rn*buck2i[typej]) :
+ (r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
+ }
+ else force_buck = 0.0;
+
+ fpair = (force_coul + force_buck) * r2inv;
+
+ if (rsq > cut_out_on_sq) { // switching
+ register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+ fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+ }
+
+ if (newton_pair || j < nlocal) { // force update
+ register double *fj = f0+(j+(j<<1)), f;
+ fi[0] += f = d[0]*fpair; fj[0] -= f;
+ fi[1] += f = d[1]*fpair; fj[1] -= f;
+ fi[2] += f = d[2]*fpair; fj[2] -= f;
+ }
+ else {
+ fi[0] += d[0]*fpair;
+ fi[1] += d[1]*fpair;
+ fi[2] += d[2]*fpair;
+ }
+ }
+ }
}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const thr)
+{
+ double r, rsq, r2inv, force_coul = 0.0, force_buck, fpair;
+
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ const double * const q = atom->q;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+ const double * const special_coul = force->special_coul;
+ const double * const special_lj = force->special_lj;
+ const double qqrd2e = force->qqrd2e;
+
+ const double *x0 = x[0];
+ double *f0 = f[0], *fi = 0;
+
+ int *ilist = list->ilist;
+
+ const int newton_pair = force->newton_pair;
+
+ const double cut_in_off = cut_respa[0];
+ const double cut_in_on = cut_respa[1];
+ const double cut_out_on = cut_respa[2];
+ const double cut_out_off = cut_respa[3];
+
+ const double cut_in_diff = cut_in_on - cut_in_off;
+ const double cut_out_diff = cut_out_off - cut_out_on;
+ const double cut_in_off_sq = cut_in_off*cut_in_off;
+ const double cut_in_on_sq = cut_in_on*cut_in_on;
+ const double cut_out_on_sq = cut_out_on*cut_out_on;
+ const double cut_out_off_sq = cut_out_off*cut_out_off;
+
+ int *jneigh, *jneighn, typei, typej, ni;
+ const int order1 = (ewald_order|(ewald_off^-1))&(1<<1);
+ int i, j, ii;
+ double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi;
+ vector xi, d;
+
+ for (ii = iifrom; ii < iito; ++ii) { // loop over my atoms
+ i = ilist[ii]; fi = f0+3*i;
+ if (order1) qri = qqrd2e*q[i];
+ memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
+ cut_bucksqi = cut_bucksq[typei = type[i]];
+ buck1i = buck1[typei]; buck2i = buck2[typei]; rhoinvi = rhoinv[typei];
+ jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+
+ for (; jneigh= cut_out_off_sq) continue;
+ if (rsq <= cut_in_off_sq) continue;
+ r2inv = 1.0/rsq;
+ r = sqrt(rsq);
+
+ if (order1 && (rsq < cut_coulsq)) // coulombic
+ force_coul = ni == 0 ?
+ qri*q[j]/r : qri*q[j]/r*special_coul[ni];
+
+ if (rsq < cut_bucksqi[typej = type[j]]) { // buckingham
+ register double rn = r2inv*r2inv*r2inv,
+ expr = exp(-r*rhoinvi[typej]);
+ force_buck = ni == 0 ?
+ (r*expr*buck1i[typej]-rn*buck2i[typej]) :
+ (r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
+ }
+ else force_buck = 0.0;
+
+ fpair = (force_coul + force_buck) * r2inv;
+
+ if (rsq < cut_in_on_sq) { // switching
+ register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+ fpair *= rsw*rsw*(3.0 - 2.0*rsw);
+ }
+ if (rsq > cut_out_on_sq) {
+ register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+ fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+ }
+
+ if (newton_pair || j < nlocal) { // force update
+ register double *fj = f0+(j+(j<<1)), f;
+ fi[0] += f = d[0]*fpair; fj[0] -= f;
+ fi[1] += f = d[1]*fpair; fj[1] -= f;
+ fi[2] += f = d[2]*fpair; fj[2] -= f;
+ }
+ else {
+ fi[0] += d[0]*fpair;
+ fi[1] += d[1]*fpair;
+ fi[2] += d[2]*fpair;
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+template < const int EVFLAG, const int EFLAG,
+ const int NEWTON_PAIR, const int CTABLE, const int DISPTABLE, const int ORDER1, const int ORDER6 >
+void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr)
+{
+ double evdwl,ecoul,fpair,fvirial;
+ evdwl = ecoul = 0.0;
+
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ const double * const q = atom->q;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+ const double * const special_coul = force->special_coul;
+ const double * const special_lj = force->special_lj;
+ const double qqrd2e = force->qqrd2e;
+
+ const double *x0 = x[0];
+ double *f0 = f[0], *fi = f0;
+
+ int *ilist = list->ilist;
+
+ int i, j, ii;
+ int *jneigh, *jneighn, typei, typej, ni, respa_flag;
+ double qi = 0.0, qri = 0.0;
+ double *cutsqi, *cut_bucksqi, *buck1i, *buck2i, *buckai, *buckci, *rhoinvi, *offseti;
+ double r, rsq, r2inv, force_coul, force_buck;
+ double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2;
+ double respa_buck = 0.0, respa_coul = 0.0, frespa = 0.0;
+ vector xi, d;
+
+ const double cut_in_off = cut_respa[2];
+ const double cut_in_on = cut_respa[3];
+
+ const double cut_in_diff = cut_in_on - cut_in_off;
+ const double cut_in_off_sq = cut_in_off*cut_in_off;
+ const double cut_in_on_sq = cut_in_on*cut_in_on;
+
+ for (ii = iiform; ii < iito; ++ii) { // loop over my atoms
+ i = ilist[ii]; fi = f0+3*i;
+ if (ORDER1) qri = (qi = q[i])*qqrd2e; // initialize constants
+ offseti = offset[typei = type[i]];
+ buck1i = buck1[typei]; buck2i = buck2[typei];
+ buckai = buck_a[typei]; buckci = buck_c[typei]; rhoinvi = rhoinv[typei];
+ cutsqi = cutsq[typei]; cut_bucksqi = cut_bucksq[typei];
+ memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
+ jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+
+ for (; jneigh= cutsqi[typej = type[j]]) continue;
+ r2inv = 1.0/rsq;
+ r = sqrt(rsq);
+
+ frespa = 1.0; //check whether and how to compute respa corrections
+ respa_coul = 0.0;
+ respa_buck = 0.0;
+ respa_flag = rsq < cut_in_on_sq ? 1 : 0;
+ if (respa_flag && (rsq > cut_in_off_sq)) {
+ register double rsw = (r-cut_in_off)/cut_in_diff;
+ frespa = 1-rsw*rsw*(3.0-2.0*rsw);
+ }
+
+ if (ORDER1 && (rsq < cut_coulsq)) { // coulombic
+ if (!CTABLE || rsq <= tabinnersq) { // series real space
+ register double s = qri*q[j];
+ if (respa_flag) // correct for respa
+ respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
+ register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+ if (ni == 0) {
+ s *= g_ewald*exp(-x*x);
+ force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
+ if (EFLAG) ecoul = t;
+ }
+ else { // correct for special
+ register double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
+ force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-ri-respa_coul;
+ if (EFLAG) ecoul = t-ri;
+ }
+ } // table real space
+ else {
+ if (respa_flag) {
+ register double s = qri*q[j];
+ respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
+ }
+ register union_int_float_t t;
+ t.f = rsq;
+ register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+ register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+ if (ni == 0) {
+ force_coul = qiqj*(ftable[k]+f*dftable[k]);
+ if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
+ }
+ else { // correct for special
+ t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
+ force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
+ if (EFLAG) {
+ t.f = (1.0-special_coul[ni])*(ptable[k]+f*dptable[k]);
+ ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
+ }
+ }
+ }
+ }
+ else force_coul = respa_coul = ecoul = 0.0;
+
+ if (rsq < cut_bucksqi[typej]) { // buckingham
+ register double rn = r2inv*r2inv*r2inv,
+ expr = exp(-r*rhoinvi[typej]);
+ if (respa_flag) respa_buck = ni == 0 ? // correct for respa
+ frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) :
+ frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
+ if (ORDER6) { // long-range form
+ if (!DISPTABLE || rsq <= tabinnerdispsq) {
+ register double x2 = g2*rsq, a2 = 1.0/x2;
+ x2 = a2*exp(-x2)*buckci[typej];
+ if (ni == 0) {
+ force_buck =
+ r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq-respa_buck;
+ if (EFLAG) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
+ }
+ else { // correct for special
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej]-
+ g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]-respa_buck;
+ if (EFLAG) evdwl = f*expr*buckai[typej] -
+ g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
+ }
+ }
+ else { // table real space
+ register union_int_float_t disp_t;
+ disp_t.f = rsq;
+ register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+ register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+ register double rn = r2inv*r2inv*r2inv;
+ if (ni == 0) {
+ force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]-respa_buck;
+ if (EFLAG) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
+ }
+ else { //special case
+ register double f = special_lj[ni], t = rn*(1.0-f);
+ force_buck = f*r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]+t*buck2i[typej]-respa_buck;
+ if (EFLAG) evdwl = f*expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
+ }
+ }
+ }
+ else { // cut form
+ if (ni == 0) {
+ force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]-respa_buck;
+ if (EFLAG)
+ evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej];
+ }
+ else { // correct for special
+ register double f = special_lj[ni];
+ force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej])-respa_buck;
+ if (EFLAG)
+ evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
+ }
+ }
+ }
+ else force_buck = respa_buck = evdwl = 0.0;
+
+ fpair = (force_coul+force_buck)*r2inv;
+
+ if (NEWTON_PAIR || j < nlocal) {
+ register double *fj = f0+(j+(j<<1)), f;
+ fi[0] += f = d[0]*fpair; fj[0] -= f;
+ fi[1] += f = d[1]*fpair; fj[1] -= f;
+ fi[2] += f = d[2]*fpair; fj[2] -= f;
+ }
+ else {
+ fi[0] += d[0]*fpair;
+ fi[1] += d[1]*fpair;
+ fi[2] += d[2]*fpair;
+ }
+
+ if (EVFLAG) {
+ fvirial = (force_coul + force_buck + respa_coul + respa_buck)*r2inv;
+ ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
+ evdwl,ecoul,fvirial,d[0],d[1],d[2],thr);
+ }
+ }
+ }
+}
+
diff --git a/src/USER-OMP/pair_buck_long_coul_long_omp.h b/src/USER-OMP/pair_buck_long_coul_long_omp.h
index ebb94ec4c5..2dc017d8fb 100644
--- a/src/USER-OMP/pair_buck_long_coul_long_omp.h
+++ b/src/USER-OMP/pair_buck_long_coul_long_omp.h
@@ -11,10 +11,6 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
-/* ----------------------------------------------------------------------
- Contributing author: Axel Kohlmeyer (Temple U)
-------------------------------------------------------------------------- */
-
#ifdef PAIR_CLASS
PairStyle(buck/long/coul/long/omp,PairBuckLongCoulLongOMP)
@@ -35,14 +31,83 @@ class PairBuckLongCoulLongOMP : public PairBuckLongCoulLong, public ThrOMP {
PairBuckLongCoulLongOMP(class LAMMPS *);
virtual void compute(int, int);
- virtual double memory_usage();
+ virtual void compute_inner();
+ virtual void compute_middle();
+ virtual void compute_outer(int, int);
private:
- template
- void eval(int ifrom, int ito, ThrData * const thr);
+ template
+ void eval(int, int, ThrData * const);
+
+ template
+ void eval_outer(int, int, ThrData * const);
+
+
+ void eval_inner(int, int, ThrData *const);
+ void eval_middle(int, int, ThrData *const);
};
}
#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.
+
+W: Geometric mixing assumed for 1/r^6 coefficients
+
+Self-explanatory.
+
+W: Using largest cutoff for buck/long/coul/long
+
+Self-exlanatory.
+
+E: Cutoffs missing in pair_style buck/long/coul/long
+
+Self-exlanatory.
+
+E: LJ6 off not supported in pair_style buck/long/coul/long
+
+Self-exlanatory.
+
+E: Coulomb cut not supported in pair_style buck/long/coul/coul
+
+Must use long-range Coulombic interactions.
+
+E: Only one cutoff allowed when requesting all long
+
+Self-explanatory.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair style buck/long/coul/long requires atom attribute q
+
+The atom style defined does not have this attribute.
+
+E: Pair style requires a KSpace style
+
+No kspace style is defined.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+E: Pair cutoff < Respa interior cutoff
+
+One or more pairwise cutoffs are too short to use with the specified
+rRESPA cutoffs.
+
+*/
diff --git a/src/USER-OMP/pppm_disp_omp.cpp b/src/USER-OMP/pppm_disp_omp.cpp
index 8a571115e7..7a60ad57f2 100644
--- a/src/USER-OMP/pppm_disp_omp.cpp
+++ b/src/USER-OMP/pppm_disp_omp.cpp
@@ -284,11 +284,13 @@ void PPPMDispOMP::compute_gf_6()
sqk = pow(qx,2.0) + pow(qy,2.0) + pow(qz,2.0);
- denominator = gf_denom(snx2,sny2,snz2, gf_b_6, order_6);
- rtsqk = sqrt(sqk);
- term = (1-2*sqk*inv2ew*inv2ew)*sx*sy*sz +
- 2*sqk*rtsqk*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtsqk*inv2ew);
- greensfn_6[nn] = numerator*term*wx*wy*wz/denominator;
+ if (sqk != 0.0) {
+ denominator = gf_denom(snx2,sny2,snz2, gf_b_6, order_6);
+ rtsqk = sqrt(sqk);
+ term = (1-2*sqk*inv2ew*inv2ew)*sx*sy*sz +
+ 2*sqk*rtsqk*inv2ew*inv2ew*inv2ew*rtpi*erfc(rtsqk*inv2ew);
+ greensfn_6[nn] = numerator*term*wx*wy*wz/denominator;
+ } else greensfn_6[nn] = 0.0;
}
}
}
diff --git a/src/pair_dpd_tstat.cpp b/src/pair_dpd_tstat.cpp
index 257f91d976..e09d332c8e 100644
--- a/src/pair_dpd_tstat.cpp
+++ b/src/pair_dpd_tstat.cpp
@@ -169,7 +169,8 @@ void PairDPDTstat::settings(int narg, char **arg)
void PairDPDTstat::coeff(int narg, char **arg)
{
- if (narg < 3 || narg > 4) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg < 3 || narg > 4)
+ error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
diff --git a/src/pair_lj_cut_coul_cut.cpp b/src/pair_lj_cut_coul_cut.cpp
index ce5e631bf9..cebae3d975 100644
--- a/src/pair_lj_cut_coul_cut.cpp
+++ b/src/pair_lj_cut_coul_cut.cpp
@@ -209,7 +209,8 @@ void PairLJCutCoulCut::settings(int narg, char **arg)
void PairLJCutCoulCut::coeff(int narg, char **arg)
{
- if (narg < 4 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg < 4 || narg > 6)
+ error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
diff --git a/src/version.h b/src/version.h
index 0aef617682..ce1579d1af 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "11 Mar 2013"
+#define LAMMPS_VERSION "14 Mar 2013"
diff --git a/tools/restart2data.cpp b/tools/restart2data.cpp
index 30edfe59e6..042c89170e 100644
--- a/tools/restart2data.cpp
+++ b/tools/restart2data.cpp
@@ -2153,6 +2153,7 @@ void pair(FILE *fp, Data &data, char *style, int flag)
(strcmp(style,"lj/cut/coul/debye") == 0) ||
(strcmp(style,"lj/cut/coul/long") == 0) ||
(strcmp(style,"lj/cut/coul/long/proxy") == 0) ||
+ (strcmp(style,"lj/cut/tip4p/cut") == 0) ||
(strcmp(style,"lj/cut/tip4p/long") == 0) ||
(strcmp(style,"lj/long/coullong") == 0)) {
@@ -2201,6 +2202,18 @@ void pair(FILE *fp, Data &data, char *style, int flag)
int tail_flag = read_int(fp);
int ncoultablebits = read_int(fp);
double tabinner = read_double(fp);
+ } else if (strcmp(style,"lj/cut/tip4p/cut") == 0) {
+ m = 0;
+ int typeO = read_int(fp);
+ int typeH = read_int(fp);
+ int typeB = read_int(fp);
+ int typeA = read_int(fp);
+ double qdist = read_double(fp);
+ double cut_lj_global = read_double(fp);
+ double cut_lj_coul = read_double(fp);
+ int offset_flag = read_int(fp);
+ int mix_flag = read_int(fp);
+ int tail_flag = read_int(fp);
} else if (strcmp(style,"lj/cut/tip4p/long") == 0) {
m = 0;
int typeO = read_int(fp);
@@ -3358,6 +3371,7 @@ void Data::write(FILE *fp, FILE *fp2, int write_coeffs, int write_vels)
(strcmp(pair_style,"lj/cut/coul/debye") == 0) ||
(strcmp(pair_style,"lj/cut/coul/long") == 0) ||
(strcmp(pair_style,"lj/cut/coul/long/proxy") == 0) ||
+ (strcmp(pair_style,"lj/cut/tip4p/cut") == 0) ||
(strcmp(pair_style,"lj/cut/tip4p/long") == 0) ||
(strcmp(pair_style,"lj/long/coul/long") == 0)) {
for (int i = 1; i <= ntypes; i++)