ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'lammps-ro/master' into lammps-icms

Browse Source 
Resolved Conflicts:
	doc/Manual.txt
	src/USER-OMP/msm_omp.h
This commit is contained in:
Axel Kohlmeyer
2013-03-16 10:18:47 +01:00
parent 1bf81d24a4 0a97cb8e9a
commit a3603ca40f
23 changed files with 1532 additions and 349 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/Manual.html
View File

@ -22,7 +22,7 @@
<CENTER><H3>LAMMPS-ICMS Documentation
</H3></CENTER>
<CENTER><H4>11 Mar 2013 version
<CENTER><H4>14 Mar 2013 version
</H4></CENTER>
<H4>Version info:
</H4>

2
doc/Manual.txt
View File

@ -18,7 +18,7 @@
<H1></H1>
LAMMPS-ICMS Documentation :c,h3
11 Mar 2013 version :c,h4
14 Mar 2013 version :c,h4
Version info: :h4

134
doc/Section_howto.html
View File

@ -445,41 +445,34 @@ set to 0.0, it corresponds to the original 1983 TIP3P model
<A HREF = "#Jorgensen">(Jorgensen)</A>.
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -0.834<BR>
H charge = 0.417 <BR>
</P>
<P>LJ epsilon of OO = 0.1521<BR>
H mass = 1.008<BR>
O charge = -0.834<BR>
H charge = 0.417<BR>
LJ epsilon of OO = 0.1521<BR>
LJ sigma of OO = 3.1507<BR>
LJ epsilon of HH = 0.0460<BR>
LJ sigma of HH = 0.4000<BR>
LJ epsilon of OH = 0.0836<BR>
LJ sigma of OH = 1.7753 <BR>
</P>
<P>K of OH bond = 450<BR>
r0 of OH bond = 0.9572 <BR>
</P>
<P>K of HOH angle = 55<BR>
LJ sigma of OH = 1.7753<BR>
K of OH bond = 450<BR>
r0 of OH bond = 0.9572<BR>
K of HOH angle = 55<BR>
theta of HOH angle = 104.52 <BR>
</P>
<P>These are the parameters to use for TIP3P with a long-range Coulombic
solver (Ewald or PPPM in LAMMPS), see <A HREF = "#Price">(Price)</A> for details:
solver (e.g. Ewald or PPPM in LAMMPS), see <A HREF = "#Price">(Price)</A> for
details:
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -0.830<BR>
H charge = 0.415 <BR>
</P>
<P>LJ epsilon of OO = 0.102<BR>
H mass = 1.008<BR>
O charge = -0.830<BR>
H charge = 0.415<BR>
LJ epsilon of OO = 0.102<BR>
LJ sigma of OO = 3.188<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
</P>
<P>K of OH bond = 450<BR>
r0 of OH bond = 0.9572 <BR>
</P>
<P>K of HOH angle = 55<BR>
LJ epsilon, sigma of OH, HH = 0.0<BR>
K of OH bond = 450<BR>
r0 of OH bond = 0.9572<BR>
K of HOH angle = 55<BR>
theta of HOH angle = 104.52 <BR>
</P>
<P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">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 <I>harmonic</I> and an
angle style of <I>harmonic</I> or <I>charmm</I> should also be used.
</P>
<P>A TIP4P model is run with LAMMPS using two commands:
<P>A TIP4P model is run with LAMMPS using either this command
for a cutoff model:
</P>
<P><A HREF = "pair_lj.html">pair_style lj/cut/tip4p/cut</A>
</P>
<P>or these two commands for a long-range model:
</P>
<UL><LI><A HREF = "pair_lj.html">pair_style lj/cut/tip4p/long</A>
<LI><A HREF = "kspace_style.html">kspace_style pppm/tip4p</A>
</UL>
<P>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 <A HREF = "fix_shake.html">fix shake</A> command.
<P>For both models, the bond lengths and bond angles should be held fixed
using the <A HREF = "fix_shake.html">fix shake</A> command.
</P>
<P>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 <A HREF = "pair_style.html">pair_style</A> command, not as part of the pair
coefficients.
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -1.040<BR>
H charge = 0.520 <BR>
</P>
<P>r0 of OH bond = 0.9572<BR>
H mass = 1.008<BR>
O charge = -1.040<BR>
H charge = 0.520<BR>
r0 of OH bond = 0.9572<BR>
theta of HOH angle = 104.52 <BR>
</P>
<P>OM distance = 0.15 <BR>
</P>
<P>LJ epsilon of O-O = 0.1550<BR>
OM distance = 0.15<BR>
LJ epsilon of O-O = 0.1550<BR>
LJ sigma of O-O = 3.1536<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
LJ epsilon, sigma of OH, HH = 0.0<BR>
Coulombic cutoff = 8.5 <BR>
</P>
<P>These are the parameters to use for TIP4P with a long-range Coulombic
solver (Ewald or PPPM in LAMMPS):
<P>For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
http://dx.doi.org/10.1063/1.1931662) these values can be used:
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
H mass = 1.008<BR>
O charge = -1.1794<BR>
H charge = 0.5897<BR>
r0 of OH bond = 0.9572<BR>
theta of HOH angle = 104.52<BR>
OM distance = 0.1577<BR>
LJ epsilon of O-O = 0.21084<BR>
LJ sigma of O-O = 3.1668<BR>
LJ epsilon, sigma of OH, HH = 0.0<BR>
Coulombic cutoff = 8.5 <BR>
</P>
<P>O charge = -1.0484<BR>
H charge = 0.5242 <BR>
<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:
</P>
<P>r0 of OH bond = 0.9572<BR>
theta of HOH angle = 104.52 <BR>
<P>O mass = 15.9994<BR>
H mass = 1.008<BR>
O charge = -1.1128<BR>
H charge = 0.5564<BR>
r0 of OH bond = 0.9572<BR>
theta of HOH angle = 104.52<BR>
OM distance = 0.1546<BR>
LJ epsilon of O-O = 0.1852<BR>
LJ sigma of O-O = 3.1589<BR>
LJ epsilon, sigma of OH, HH = 0.0<BR>
Coulombic cutoff = 8.5 <BR>
</P>
<P>OM distance = 0.1250 <BR>
<P>These are the parameters to use for TIP4P with a long-range Coulombic
solver (e.g. Ewald or PPPM in LAMMPS):
</P>
<P>LJ epsilon of O-O = 0.16275<BR>
<P>O mass = 15.9994<BR>
H mass = 1.008<BR>
O charge = -1.0484<BR>
H charge = 0.5242<BR>
r0 of OH bond = 0.9572<BR>
theta of HOH angle = 104.52<BR>
OM distance = 0.1250<BR>
LJ epsilon of O-O = 0.16275<BR>
LJ sigma of O-O = 3.16435<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
</P>
@ -574,16 +593,13 @@ used.
atoms and the water molecule to run a rigid SPC model.
</P>
<P>O mass = 15.9994<BR>
H mass = 1.008 <BR>
</P>
<P>O charge = -0.820<BR>
H charge = 0.410 <BR>
</P>
<P>LJ epsilon of OO = 0.1553<BR>
H mass = 1.008<BR>
O charge = -0.820<BR>
H charge = 0.410<BR>
LJ epsilon of OO = 0.1553<BR>
LJ sigma of OO = 3.166<BR>
LJ epsilon, sigma of OH, HH = 0.0 <BR>
</P>
<P>r0 of OH bond = 1.0<BR>
LJ epsilon, sigma of OH, HH = 0.0<BR>
r0 of OH bond = 1.0<BR>
theta of HOH angle = 109.47 <BR>
</P>
<P>Note that as originally proposed, the SPC model was run with a 9

108
doc/Section_howto.txt
View File

@ -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

5
doc/compute_temp_asphere.html
View File

@ -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 <A HREF = "pair_gayberne.html">GayBerne
model</A> 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.
</P>
<P>The translational kinetic energy is computed the same as is described
by the <A HREF = "compute_temp.html">compute temp</A> command. The rotational

5
doc/compute_temp_asphere.txt
View File

@ -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

6
doc/pair_buck_long.html
View File

@ -140,9 +140,9 @@ interaction, assuming <I>flag_buck</I> is <I>cut</I>.
shift option for the energy of the Buckingham portion of the pair
interaction.
</P>
<P>This pair style does not support the <A HREF = "pair_modify.html">pair_modify</A>
table option since a tabulation capability has not yet been added to
this potential.
<P>This pair style supports the <A HREF = "pair_modify.html">pair_modify</A> table and
table/disp options since they can tabulate the short-range portion of
the long-range Coulombic and dispersion interactions.
</P>
<P>This pair style write its information to <A HREF = "restart.html">binary restart
files</A>, so pair_style and pair_coeff commands do not need

6
doc/pair_buck_long.txt
View File

@ -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

48
doc/pair_lj.html
View File

@ -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
</PRE>
<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>
<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
@ -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.
</P>
<P>Style <I>lj/cut/tip4p/long</I> implements the TIP4P water model of
<A HREF = "#Jorgensen">(Jorgensen)</A>, 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
<I>lj/cut/tip4p/cut</I> is similar to <I>lj/cut/tip4p/long</I>, only it uses a
cutoff for Coulomb interactions.
<P>Styles <I>lj/cut/tip4p/cut</I> and <I>lj/cut/tip4p/long</I> implement the TIP4P
water model of <A HREF = "#Jorgensen">(Jorgensen)</A>, 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 <I>lj/cut/tip4p/cut</I> uses a cutoff for
Coulomb interactions; style <I>lj/cut/tip4p/long</I> is for use with a
long-range Coulombic solver (Ewald or PPPM).
</P>
<P>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.
</P>
<P>See the <A HREF = "Section_howto.html#howto_8">howto section</A> 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.
</P>
<P>For all of the <I>lj/cut</I> 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 <I>lj/cut</I>, since it
has no Coulombic terms.
</P>
<P>For <I>lj/cut/coul/long</I> and <I>lj/cut/coul/msm</I> and <I>lj/cut/tip4p/long</I>
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.
<P>For <I>lj/cut/coul/long</I> and <I>lj/cut/coul/msm</I> and <I>lj/cut/tip4p/cut</I>
and <I>lj/cut/tip4p/long</I> 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.
</P>
<HR>

48
doc/pair_lj.txt
View File

@ -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

163
src/KSPACE/pair_buck_long_coul_long.cpp
View File

@ -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 (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
qri = qqrd2e*q[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 = listinner->firstneigh[i])+listinner->numneigh[i];
jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
for (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
@ -697,15 +716,15 @@ void PairBuckLongCoulLong::compute_middle()
double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi;
vector xi, d;
ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum;
ineighn = (ineigh = list->ilist)+list->inum;
for (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
qri = qqrd2e*q[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 = listmiddle->firstneigh[i])+listmiddle->numneigh[i];
jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
for (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
@ -765,7 +784,7 @@ void PairBuckLongCoulLong::compute_middle()
void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
{
double evdwl,ecoul,fpair;
double evdwl,ecoul,fpair,fvirial;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = 0;
@ -796,7 +815,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;
ineighn = (ineigh = listouter->ilist)+listouter->inum;
ineighn = (ineigh = list->ilist)+list->inum;
for (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
@ -806,7 +825,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
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 = listouter->firstneigh[i])+listouter->numneigh[i];
jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
for (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
@ -822,9 +841,13 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
r2inv = 1.0/rsq;
r = sqrt(rsq);
if ((respa_flag = (rsq>cut_in_off_sq)&&(rsq<cut_in_on_sq))) {
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 = 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]);
}
}
}
}

6
src/KSPACE/pair_buck_long_coul_long.h
View File

@ -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;

12
src/KSPACE/pppm_disp.cpp
View File

@ -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;
}
}
}

6
src/MC/fix_gcmc.cpp
View File

@ -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);

63
src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

92
src/USER-OMP/msm_omp.h
View File

@ -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 <int, int, int> void direct_eval(int);
template <int> 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 <int, int, int> void direct_eval(int);
template <int> void direct_peratom(int);
};
}
#endif
#endif

1062
src/USER-OMP/pair_buck_long_coul_long_omp.cpp
View File

File diff suppressed because it is too large Load Diff

79
src/USER-OMP/pair_buck_long_coul_long_omp.h
View File

@ -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 <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void eval(int ifrom, int ito, ThrData * const thr);
template <const int EVFLAG, const int EFLAG,
const int NEWTON_PAIR, const int CTABLE, const int LJTABLE,
const int ORDER1, const int ORDER6 >
void eval(int, int, ThrData * const);
template <const int EVFLAG, const int EFLAG,
const int NEWTON_PAIR, const int CTABLE, const int LJTABLE,
const int ORDER1, const int ORDER6 >
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.
*/

12
src/USER-OMP/pppm_disp_omp.cpp
View File

@ -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;
}
}
}

3
src/pair_dpd_tstat.cpp
View File

@ -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;

3
src/pair_lj_cut_coul_cut.cpp
View File

@ -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;

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "11 Mar 2013"
#define LAMMPS_VERSION "14 Mar 2013"

14
tools/restart2data.cpp
View File

@ -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++)