From 7fbf90ea90b9aae42ee2686b71c4279b4d05f120 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 13 Mar 2013 15:01:49 +0000 Subject: [PATCH 01/15] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9656 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Section_howto.html | 134 +++++++++++++++++++++++------------------ doc/Section_howto.txt | 108 +++++++++++++++++++-------------- 2 files changed, 137 insertions(+), 105 deletions(-) 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 From 9897d385421655e491823e07d2cf5b00ba84cf56 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 13 Mar 2013 15:26:18 +0000 Subject: [PATCH 02/15] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9657 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/pair_dpd_tstat.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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; From 0cd14bfb847ad71d08dac2387cd6280fa6b50411 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 13 Mar 2013 15:26:52 +0000 Subject: [PATCH 03/15] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9658 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/KSPACE/pppm_disp.cpp | 3 +- src/MOLECULE/pair_lj_cut_tip4p_cut.cpp | 41 +++++++++++++++++++++++--- 2 files changed, 38 insertions(+), 6 deletions(-) diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 58cecb039f..a4d385291d 100755 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -221,12 +221,11 @@ void PPPMDisp::init() } // free all arrays previously allocated + deallocate(); deallocate_peratom(); peratom_allocate_flag = 0; - - // set scale scale = 1.0; diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp index 4bd52d293e..d7ad345b77 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,8 @@ 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 +438,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 +449,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 +529,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() From a58779babf6220e6fd0566c20cff38c311bf47b5 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 13 Mar 2013 15:30:49 +0000 Subject: [PATCH 04/15] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9659 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/pair_lj.html | 48 ++++++++++++++++++++++++++++-------------------- doc/pair_lj.txt | 48 ++++++++++++++++++++++++++++-------------------- 2 files changed, 56 insertions(+), 40 deletions(-) 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
 

From b78c4aa841671f2bf162f3de0f35d3c0388c3718 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Wed, 13 Mar 2013 15:36:13 +0000
Subject: [PATCH 05/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9660
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/pair_lj_cut_coul_cut.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

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;

From a200f0ae3814a0aa0054a2969fb8bf7cb4126fee Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Wed, 13 Mar 2013 15:43:17 +0000
Subject: [PATCH 06/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9661
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 tools/restart2data.cpp | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/tools/restart2data.cpp b/tools/restart2data.cpp
index 77852ed9a6..f5e465a488 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++)

From ba0c0a959c8e479a91358d3bb256b7d35ef9a103 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Wed, 13 Mar 2013 15:43:36 +0000
Subject: [PATCH 07/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9662
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/MOLECULE/pair_lj_cut_tip4p_cut.cpp | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
index d7ad345b77..72f809aec7 100644
--- a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
+++ b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
@@ -421,6 +421,7 @@ void PairLJCutTIP4PCut::settings(int narg, char **arg)
   cut_lj_global = force->numeric(arg[5]);
   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);
 
@@ -631,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);
 }
 
 /* ----------------------------------------------------------------------
@@ -649,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);
@@ -660,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);
 }
 
 /* ----------------------------------------------------------------------

From 788b657922e1fc853c50d01ba24fcd1ec9d76c57 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Wed, 13 Mar 2013 18:48:49 +0000
Subject: [PATCH 08/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9663
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/compute_temp_asphere.html | 5 +++++
 doc/compute_temp_asphere.txt  | 5 +++++
 2 files changed, 10 insertions(+)

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

From eeecd6980262145cb6aef6e1fbd4b4c5b03de9ad Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Thu, 14 Mar 2013 14:30:02 +0000
Subject: [PATCH 09/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9680
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/USER-OMP/msm_omp.cpp | 263 +++++++++++++++++++++++++++++++++++----
 src/USER-OMP/msm_omp.h   |  91 +++++++-------
 2 files changed, 286 insertions(+), 68 deletions(-)

diff --git a/src/USER-OMP/msm_omp.cpp b/src/USER-OMP/msm_omp.cpp
index 480e52ce2f..18e2086a3b 100755
--- a/src/USER-OMP/msm_omp.cpp
+++ b/src/USER-OMP/msm_omp.cpp
@@ -112,12 +112,16 @@ void MSMOMP::direct(int n)
         direct_eval<0,0,0>(n);
     }
   }
+
+  if (vflag_atom)
+    direct_peratom<1>(n);
+  else
+    direct_peratom<0>(n);
 }
 
 template 
 void MSMOMP::direct_eval(const int nn)
 {
-  double * _noalias const * _noalias const * _noalias const egridn  = egrid[nn];
   const double * _noalias const * _noalias const * _noalias const qgridn = qgrid[nn];
   const double * _noalias const g_directn = g_direct[nn];
   const double * _noalias const v0_directn = v0_direct[nn];
@@ -157,14 +161,14 @@ void MSMOMP::direct_eval(const int nn)
 #pragma omp parallel default(none) reduction(+:v0,v1,v2,v3,v4,v5,emsm)
 #endif
   {
-    double qtmp,esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
+    double esum,v0sum,v1sum,v2sum,v3sum,v4sum,v5sum;
     int i,ifrom,ito,tid,icx,icy,icz,ix,iy,iz,k;
 
     loop_setup_thr(ifrom, ito, tid, inum, comm->nthreads);
 
     for (i = ifrom; i < ito; ++i) {
 
-      // infer outer loop indices icx, icy, icz from master loop index
+      // infer outer loop indices icx, icy, icz from master loop index i
 
       icz = i/(numy*numx);
       icy = (i - icz*numy*numx) / numx;
@@ -173,41 +177,102 @@ void MSMOMP::direct_eval(const int nn)
       icy += nylo_inn;
       icx += nxlo_inn;
       
-      const int kmin = zper ? nzlo_direct : MAX(nzlo_direct,alphan - icz);
       const int kmax = zper ? nzhi_direct : MIN(nzhi_direct,betazn - icz);
       const int jmin = yper ? nylo_direct : MAX(nylo_direct,alphan - icy);
       const int jmax = yper ? nyhi_direct : MIN(nyhi_direct,betayn - icy);
       const int imin = xper ? nxlo_direct : MAX(nxlo_direct,alphan - icx);
       const int imax = xper ? nxhi_direct : MIN(nxhi_direct,betaxn - icx);
 
+      const double qtmp = qgridn[icz][icy][icx]; // charge on center grid point
+
       esum = 0.0;
       if (VFLAG_GLOBAL || VFLAG_ATOM)
         v0sum = v1sum = v2sum = v3sum = v4sum = v5sum = 0.0;
 
-      for (iz = kmin; iz <= kmax; iz++) {
+      // use hemisphere to avoid double computation of pair-wise
+      //   interactions in direct sum (no computations in -z direction)
+
+      for (iz = 1; iz <= kmax; iz++) {
         const int kk = icz+iz;
         const int zk = (iz + nzhi_direct)*ny;
         for (iy = jmin; iy <= jmax; iy++) {
           const int jj = icy+iy;
           const int zyk = (zk + iy + nyhi_direct)*nx;
-          const double * _noalias const qgridnkj = & qgridn[kk][jj][icx];
+          const double * _noalias const qgridnkj = &qgridn[kk][jj][icx];
           for (ix = imin; ix <= imax; ix++) {
-            qtmp = qgridnkj[ix];
+            const double qtmp2 = qgridnkj[ix];
             k = zyk + ix + nxhi_direct;
-            esum += g_directn[k] * qtmp;
+            const double gtmp = g_directn[k];
+            esum += gtmp * qtmp2;
 
             if (VFLAG_GLOBAL || VFLAG_ATOM) {
-              v0sum += v0_directn[k] * qtmp;
-              v1sum += v1_directn[k] * qtmp;
-              v2sum += v2_directn[k] * qtmp;
-              v3sum += v3_directn[k] * qtmp;
-              v4sum += v4_directn[k] * qtmp;
-              v5sum += v5_directn[k] * qtmp;
+              v0sum += v0_directn[k] * qtmp2;
+              v1sum += v1_directn[k] * qtmp2;
+              v2sum += v2_directn[k] * qtmp2;
+              v3sum += v3_directn[k] * qtmp2;
+              v4sum += v4_directn[k] * qtmp2;
+              v5sum += v5_directn[k] * qtmp2;
             }
           }
         }
       }
-      egridn[icz][icy][icx] = esum;
+
+      // iz=0
+
+      const int zk = nzhi_direct*ny;
+      for (iy = 1; iy <= jmax; iy++) {
+        const int jj = icy+iy;
+        const int zyk = (zk + iy + nyhi_direct)*nx;
+        const double * _noalias const qgridnkj = &qgridn[icz][jj][icx];
+        for (ix = imin; ix <= imax; ix++) {
+          const double qtmp2 = qgridnkj[ix];
+          k = zyk + ix + nxhi_direct;
+          const double gtmp = g_directn[k];
+          esum += gtmp * qtmp2;
+
+          if (VFLAG_GLOBAL || VFLAG_ATOM) {
+            v0sum += v0_directn[k] * qtmp2;
+            v1sum += v1_directn[k] * qtmp2;
+            v2sum += v2_directn[k] * qtmp2;
+            v3sum += v3_directn[k] * qtmp2;
+            v4sum += v4_directn[k] * qtmp2;
+            v5sum += v5_directn[k] * qtmp2;
+          }
+        }
+      }
+
+      // iz=0, iy=0
+
+      const int zyk = (zk + nyhi_direct)*nx;
+      const double * _noalias const qgridnkj = &qgridn[icz][icy][icx];
+      for (ix = 1; ix <= imax; ix++) {
+        const double qtmp2 = qgridnkj[ix];
+        k = zyk + ix + nxhi_direct;
+        const double gtmp = g_directn[k];
+        esum += gtmp * qtmp2;
+
+        if (VFLAG_GLOBAL || VFLAG_ATOM) {
+          v0sum += v0_directn[k] * qtmp2;
+          v1sum += v1_directn[k] * qtmp2;
+          v2sum += v2_directn[k] * qtmp2;
+          v3sum += v3_directn[k] * qtmp2;
+          v4sum += v4_directn[k] * qtmp2;
+          v5sum += v5_directn[k] * qtmp2;
+        }
+      }
+
+      // iz=0, iy=0, ix=0
+
+      const double qtmp2 = qgridnkj[0];
+      k = zyk + nxhi_direct;
+      const double gtmp = g_directn[k];
+      esum += 0.5 * gtmp * qtmp2;
+
+      // virial is zero for iz=0, iy=0, ix=0
+
+      // accumulate per-atom energy/virial
+
+      egrid[n][icz][icy][icx] = esum;
 
       if (VFLAG_ATOM) {
         v0grid[n][icz][icy][icx] = v0sum;
@@ -219,15 +284,15 @@ void MSMOMP::direct_eval(const int nn)
       }
 
       if (EFLAG_GLOBAL || VFLAG_GLOBAL) {
-        qtmp = qgridn[icz][icy][icx];
-        if (EFLAG_GLOBAL) emsm += esum * qtmp;
+        const double qtmp3 = qgridn[icz][icy][icx];
+        if (EFLAG_GLOBAL) emsm += 2.0 * esum * qtmp3;
         if (VFLAG_GLOBAL) {
-          v0 += v0sum * qtmp;
-          v1 += v1sum * qtmp;
-          v2 += v2sum * qtmp;
-          v3 += v3sum * qtmp;
-          v4 += v4sum * qtmp;
-          v5 += v5sum * qtmp;
+          v0 += 2.0 * v0sum * qtmp3;
+          v1 += 2.0 * v1sum * qtmp3;
+          v2 += 2.0 * v2sum * qtmp3;
+          v3 += 2.0 * v3sum * qtmp3;
+          v4 += 2.0 * v4sum * qtmp3;
+          v5 += 2.0 * v5sum * qtmp3;
         }
       }
     }
@@ -245,3 +310,155 @@ void MSMOMP::direct_eval(const int nn)
     }
   }
 }
+
+template 
+void MSMOMP::direct_peratom(const int nn)
+{
+  double * _noalias const * _noalias const * _noalias const egridn  = egrid[nn];
+  double * _noalias const * _noalias const * _noalias const v0gridn  = v0grid[nn];
+  double * _noalias const * _noalias const * _noalias const v1gridn  = v1grid[nn];
+  double * _noalias const * _noalias const * _noalias const v2gridn  = v2grid[nn];
+  double * _noalias const * _noalias const * _noalias const v3gridn  = v3grid[nn];
+  double * _noalias const * _noalias const * _noalias const v4gridn  = v4grid[nn];
+  double * _noalias const * _noalias const * _noalias const v5gridn  = v5grid[nn];
+  const double * _noalias const * _noalias const * _noalias const qgridn = qgrid[nn];
+  const double * _noalias const g_directn = g_direct[nn];
+  const double * _noalias const v0_directn = v0_direct[nn];
+  const double * _noalias const v1_directn = v1_direct[nn];
+  const double * _noalias const v2_directn = v2_direct[nn];
+  const double * _noalias const v3_directn = v3_direct[nn];
+  const double * _noalias const v4_directn = v4_direct[nn];
+  const double * _noalias const v5_directn = v5_direct[nn];
+
+
+  const int alphan = alpha[nn];
+  const int betaxn = betax[nn];
+  const int betayn = betay[nn];
+  const int betazn = betaz[nn];
+
+  const int nx = nxhi_direct - nxlo_direct + 1;
+  const int ny = nyhi_direct - nylo_direct + 1;
+
+  // merge three outer loops into one
+
+  const int nzlo_inn = nzlo_in[nn];
+  const int nylo_inn = nylo_in[nn];
+  const int nxlo_inn = nxlo_in[nn];
+  const int numz = nzhi_in[nn] - nzlo_inn + 1;
+  const int numy = nyhi_in[nn] - nylo_inn + 1;
+  const int numx = nxhi_in[nn] - nxlo_inn + 1;
+  const int inum = numz*numy*numx;
+
+  const int zper = domain->zperiodic;
+  const int yper = domain->yperiodic;
+  const int xper = domain->xperiodic;
+
+  const int n=nn;
+  int i,ifrom,ito,tid,icx,icy,icz,ix,iy,iz,k;
+
+
+  for (i = 0; i < inum; ++i) {
+
+    // infer outer loop indices icx, icy, icz from master loop index i
+
+    icz = i/(numy*numx);
+    icy = (i - icz*numy*numx) / numx;
+    icx = i - icz*numy*numx - icy*numx;
+    icz += nzlo_inn;
+    icy += nylo_inn;
+    icx += nxlo_inn;
+    
+    const int kmax = zper ? nzhi_direct : MIN(nzhi_direct,betazn - icz);
+    const int jmin = yper ? nylo_direct : MAX(nylo_direct,alphan - icy);
+    const int jmax = yper ? nyhi_direct : MIN(nyhi_direct,betayn - icy);
+    const int imin = xper ? nxlo_direct : MAX(nxlo_direct,alphan - icx);
+    const int imax = xper ? nxhi_direct : MIN(nxhi_direct,betaxn - icx);
+
+    const double qtmp = qgridn[icz][icy][icx]; // charge on center grid point
+
+
+    // use hemisphere to avoid double computation of pair-wise
+    //   interactions in direct sum (no computations in -z direction)
+
+    for (iz = 1; iz <= kmax; iz++) {
+      const int kk = icz+iz;
+      const int zk = (iz + nzhi_direct)*ny;
+      for (iy = jmin; iy <= jmax; iy++) {
+        const int jj = icy+iy;
+        const int zyk = (zk + iy + nyhi_direct)*nx;
+        double * _noalias const egridnkj = &egridn[kk][jj][icx];
+        for (ix = imin; ix <= imax; ix++) {
+          k = zyk + ix + nxhi_direct;
+          const int ii = icx+ix;
+          const double gtmp = g_directn[k];
+
+          egridnkj[ix] += gtmp * qtmp;
+
+          if (VFLAG_ATOM) {
+            v0gridn[kk][jj][ii] += v0_directn[k] * qtmp;
+            v1gridn[kk][jj][ii] += v1_directn[k] * qtmp;
+            v2gridn[kk][jj][ii] += v2_directn[k] * qtmp;
+            v3gridn[kk][jj][ii] += v3_directn[k] * qtmp;
+            v4gridn[kk][jj][ii] += v4_directn[k] * qtmp;
+            v5gridn[kk][jj][ii] += v5_directn[k] * qtmp;
+          }
+        }
+      }
+    }
+
+    // iz=0
+
+    const int zk = nzhi_direct*ny;
+    for (iy = 1; iy <= jmax; iy++) {
+      const int jj = icy+iy;
+      const int zyk = (zk + iy + nyhi_direct)*nx;
+      double * _noalias const egridnkj = &egridn[icz][jj][icx];
+      for (ix = imin; ix <= imax; ix++) {
+        k = zyk + ix + nxhi_direct;
+        const int ii = icx+ix;
+        const double gtmp = g_directn[k];
+
+        egridnkj[ix] += gtmp * qtmp;
+
+        if (VFLAG_ATOM) {
+          v0gridn[icz][jj][ii] += v0_directn[k] * qtmp;
+          v1gridn[icz][jj][ii] += v1_directn[k] * qtmp;
+          v2gridn[icz][jj][ii] += v2_directn[k] * qtmp;
+          v3gridn[icz][jj][ii] += v3_directn[k] * qtmp;
+          v4gridn[icz][jj][ii] += v4_directn[k] * qtmp;
+          v5gridn[icz][jj][ii] += v5_directn[k] * qtmp;
+        }
+      }
+    }
+
+    // iz=0, iy=0
+
+    const int zyk = (zk + nyhi_direct)*nx;
+    double * _noalias const egridnkj = &egridn[icz][icy][icx];
+    for (ix = 1; ix <= imax; ix++) {
+      k = zyk + ix + nxhi_direct;
+      const int ii = icx+ix;
+      const double gtmp = g_directn[k];
+
+      egridnkj[ix] += gtmp * qtmp;
+
+      if (VFLAG_ATOM) {
+        v0gridn[icz][icy][ii] += v0_directn[k] * qtmp;
+        v1gridn[icz][icy][ii] += v1_directn[k] * qtmp;
+        v2gridn[icz][icy][ii] += v2_directn[k] * qtmp;
+        v3gridn[icz][icy][ii] += v3_directn[k] * qtmp;
+        v4gridn[icz][icy][ii] += v4_directn[k] * qtmp;
+        v5gridn[icz][icy][ii] += v5_directn[k] * qtmp;
+      }
+    }
+
+    // iz=0, iy=0, ix=0
+
+    k = zyk + nxhi_direct;
+    const double gtmp = g_directn[k];
+    egridnkj[0] += 0.5 * gtmp * qtmp;
+
+    // virial is zero for iz=0, iy=0, ix=0
+
+  }
+}
diff --git a/src/USER-OMP/msm_omp.h b/src/USER-OMP/msm_omp.h
index 369e23be49..0a41c57ac3 100755
--- a/src/USER-OMP/msm_omp.h
+++ b/src/USER-OMP/msm_omp.h
@@ -1,45 +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);
-
-};
-
-}
-
-#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

From 040cedbe5aa1d9098e6f1e877ce2a47a54f48fda Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Thu, 14 Mar 2013 16:34:11 +0000
Subject: [PATCH 10/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9687
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/KSPACE/pair_buck_long_coul_long.cpp | 163 +++++++++++++++---------
 src/KSPACE/pair_buck_long_coul_long.h   |   6 +-
 src/KSPACE/pppm_disp.cpp                |  15 ++-
 3 files changed, 117 insertions(+), 67 deletions(-)

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 a4d385291d..62a22bf7a7 100755
--- a/src/KSPACE/pppm_disp.cpp
+++ b/src/KSPACE/pppm_disp.cpp
@@ -221,11 +221,12 @@ void PPPMDisp::init()
   }
 
   // free all arrays previously allocated
-
   deallocate();
   deallocate_peratom(); 
   peratom_allocate_flag = 0;
 
+
+
   // set scale
 
   scale = 1.0;
@@ -3354,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;
       }
     }
   }

From 28480c5221d53347bdc1187d704a6f2e47674ca1 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Thu, 14 Mar 2013 16:34:14 +0000
Subject: [PATCH 11/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9688
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/USER-OMP/pair_buck_long_coul_long_omp.cpp | 1062 ++++++++++++++++-
 src/USER-OMP/pair_buck_long_coul_long_omp.h   |   79 +-
 src/USER-OMP/pppm_disp_omp.cpp                |   12 +-
 3 files changed, 1083 insertions(+), 70 deletions(-)

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

From 6b884058df864e222b97beca99c08a29ceb1cdc6 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Thu, 14 Mar 2013 16:37:43 +0000
Subject: [PATCH 12/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9689
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/pair_buck_long.html | 6 +++---
 doc/pair_buck_long.txt  | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

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

From 32f9bf830620bf3135696bb5f2541ff41ebb3a5a Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Thu, 14 Mar 2013 16:46:01 +0000
Subject: [PATCH 13/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9690
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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"

From 505a0dd4f4798d151bf1bf48196377beee9c9ca9 Mon Sep 17 00:00:00 2001
From: sjplimp 
Date: Thu, 14 Mar 2013 16:46:01 +0000
Subject: [PATCH 14/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9691
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/Manual.html | 2 +-
 doc/Manual.txt  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/Manual.html b/doc/Manual.html
index 344f70c7eb..8c35c8dc96 100644
--- a/doc/Manual.html
+++ b/doc/Manual.html
@@ -22,7 +22,7 @@
 
 
LAMMPS Documentation 
 

-
11 Mar 2013 version 
+
14 Mar 2013 version 
 

 
Version info: 
 

diff --git a/doc/Manual.txt b/doc/Manual.txt
index ad85c3d217..9d678a242c 100644
--- a/doc/Manual.txt
+++ b/doc/Manual.txt
@@ -18,7 +18,7 @@
 

 
 LAMMPS Documentation :c,h3
-11 Mar 2013 version :c,h4
+14 Mar 2013 version :c,h4
 
 Version info: :h4
 

From 0a97cb8e9a6ad27e60ee526ac90a2c5ca576fcc9 Mon Sep 17 00:00:00 2001
From: pscrozi 
Date: Fri, 15 Mar 2013 00:14:36 +0000
Subject: [PATCH 15/15] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@9694
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/MC/fix_gcmc.cpp | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

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