git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9701 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/KSPACE/pair_buck_long_coul_long.cpp

@@ -236,7 +236,7 @@ void PairBuckLongCoulLong::init_style()
 
   int irequest;
 
-  if (update->whichflag == 0 && strstr(update->integrate_style,"respa")) {
+  if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
     int respa = 0;
     if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
     if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
@@ -627,14 +627,14 @@ void PairBuckLongCoulLong::compute_inner()
   double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi;
   vector xi, d;
 
-  ineighn = (ineigh = list->ilist) + list->inum;
+  ineighn = (ineigh = listinner->ilist) + listinner->inum;
   for (; ineigh<ineighn; ++ineigh) {                        // loop over my atoms
     i = *ineigh; 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];
+    jneighn = (jneigh = listinner->firstneigh[i])+listinner->numneigh[i];
 
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;
@@ -716,7 +716,7 @@ void PairBuckLongCoulLong::compute_middle()
   double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi;
   vector xi, d;
 
-  ineighn = (ineigh = list->ilist)+list->inum;
+  ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum;
 
   for (; ineigh<ineighn; ++ineigh) {                        // loop over my atoms
     i = *ineigh; fi = f0+3*i;
@@ -724,7 +724,7 @@ void PairBuckLongCoulLong::compute_middle()
     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];
+    jneighn = (jneigh = listmiddle->firstneigh[i])+listmiddle->numneigh[i];
 
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;
@@ -815,7 +815,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
   double cut_in_off_sq = cut_in_off*cut_in_off;
   double cut_in_on_sq = cut_in_on*cut_in_on;
 
-  ineighn = (ineigh = list->ilist)+list->inum;
+  ineighn = (ineigh = listouter->ilist)+listouter->inum;
 
   for (; ineigh<ineighn; ++ineigh) {                        // loop over my atoms
     i = *ineigh; fi = f0+3*i;
@@ -825,7 +825,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
     buckai = buck_a[typei]; buckci = buck_c[typei]; rhoinvi = rhoinv[typei];
     cutsqi = cutsq[typei]; cut_bucksqi = cut_bucksq[typei];
     memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
-    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+    jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i];
 
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;

src/KSPACE/pair_lj_long_coul_long.cpp

@@ -236,7 +236,7 @@ void PairLJLongCoulLong::init_style()
 
   int irequest;
 
-  if (update->whichflag == 0 && strstr(update->integrate_style,"respa")) {
+  if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
     int respa = 0;
     if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
     if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
@@ -633,13 +633,13 @@ void PairLJLongCoulLong::compute_inner()
   double qri, *cut_ljsqi, *lj1i, *lj2i;
   vector xi, d;
 
-  ineighn = (ineigh = list->ilist)+list->inum;
+  ineighn = (ineigh = listinner->ilist)+listinner->inum;
   for (; ineigh<ineighn; ++ineigh) {                        // loop over my atoms
     i = *ineigh; fi = f0+3*i;
     memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
     cut_ljsqi = cut_ljsq[typei = type[i]];
     lj1i = lj1[typei]; lj2i = lj2[typei];
-    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+    jneighn = (jneigh = listinner->firstneigh[i])+listinner->numneigh[i];
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;
       ni = sbmask(j);
@@ -720,7 +720,7 @@ void PairLJLongCoulLong::compute_middle()
   double qri, *cut_ljsqi, *lj1i, *lj2i;
   vector xi, d;
 
-  ineighn = (ineigh = list->ilist)+list->inum;
+  ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum;
 
   for (; ineigh<ineighn; ++ineigh) {                        // loop over my atoms
     i = *ineigh; fi = f0+3*i;
@@ -728,7 +728,7 @@ void PairLJLongCoulLong::compute_middle()
     memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
     cut_ljsqi = cut_ljsq[typei = type[i]];
     lj1i = lj1[typei]; lj2i = lj2[typei];
-    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+    jneighn = (jneigh = listmiddle->firstneigh[i])+listmiddle->numneigh[i];
 
     for (; jneigh<jneighn; ++jneigh) {
       j = *jneigh;
@@ -817,7 +817,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
   double cut_in_off_sq = cut_in_off*cut_in_off;
   double cut_in_on_sq = cut_in_on*cut_in_on;
   
-  ineighn = (ineigh = list->ilist)+list->inum;
+  ineighn = (ineigh = listouter->ilist)+listouter->inum;
   
   for (; ineigh<ineighn; ++ineigh) {                        // loop over my atoms
     i = *ineigh; fi = f0+3*i;
@@ -826,7 +826,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
     lj1i = lj1[typei]; lj2i = lj2[typei]; lj3i = lj3[typei]; lj4i = lj4[typei];
     cutsqi = cutsq[typei]; cut_ljsqi = cut_ljsq[typei];
     memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
-    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+    jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i];
     
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;

src/KSPACE/pair_lj_long_tip4p_long.cpp

@@ -53,7 +53,7 @@ PairLJLongTIP4PLong::PairLJLongTIP4PLong(LAMMPS *lmp) :
 {
   dispersionflag = tip4pflag = 1;
   single_enable = 0;
-  respa_enable = 0;
+  respa_enable = 1;
 
   nmax = 0;
   hneigh = NULL;
@@ -184,22 +184,40 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
       if (rsq < cut_ljsq[itype][jtype]) {			// lj
         r2inv = 1.0/rsq;
        	if (order6) {					// long-range lj
-	  register double rn = r2inv*r2inv*r2inv;
+          if (!ndisptablebits || rsq <= tabinnerdispsq) {
-	  register double x2 = g2*rsq, a2 = 1.0/x2;
+	    register double rn = r2inv*r2inv*r2inv;
-	  x2 = a2*exp(-x2)*lj4i[jtype];
+	    register double x2 = g2*rsq, a2 = 1.0/x2;
-	  if (ni == 0) {
+	    x2 = a2*exp(-x2)*lj4i[jtype];
-	    forcelj =
+	    if (ni == 0) {
-	      (rn*=rn)*lj1i[jtype]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
+	      forcelj =
-	    if (eflag)
+	        (rn*=rn)*lj1i[jtype]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
-	      evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
+	      if (eflag)
-	  }
+	        evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
-	  else {					// special case
+	    }
-	    register double f = special_lj[ni], t = rn*(1.0-f);
+	    else {					// special case
-	    forcelj = f*(rn *= rn)*lj1i[jtype]-
+	      register double f = special_lj[ni], t = rn*(1.0-f);
-	      g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype];
+	      forcelj = f*(rn *= rn)*lj1i[jtype]-
-	    if (eflag) 
+	        g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype];
-	      evdwl = f*rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[jtype];
+	      if (eflag) 
-	  }
+	        evdwl = f*rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[jtype];
+	    }
+          }
+          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) {
+              forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype];
+              if (eflag) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
+            }
+            else {					// special case
+              register double f = special_lj[ni], t = rn*(1.0-f);
+              forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype];
+              if (eflag) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
+            }
+          }
 	}
 	else {						// cut lj
 	  register double rn = r2inv*r2inv*r2inv;
@@ -437,6 +455,922 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
   }
 }
 
+/* --------------------------------------------------------------------- */
+
+void PairLJLongTIP4PLong::compute_inner()
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype,itable;
+  int iH1,iH2,jH1,jH2;
+  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
+  double r,r2inv,forcecoul,forcelj,cforce;
+  double factor_coul,factor_lj;
+  double grij,expm2,prefactor,t,erfc;
+  double xiM[3],xjM[3],fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];// f1[3];
+  double *x1,*x2;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double rsq, qri;
+
+  double cut_out_on = cut_respa[0];
+  double cut_out_off = cut_respa[1];
+
+  double cut_out_diff = cut_out_off - cut_out_on;
+  double cut_out_on_sq = cut_out_on*cut_out_on;
+  double cut_out_off_sq = cut_out_off*cut_out_off;
+
+  // reallocate hneigh & newsite if necessary
+  // initialize hneigh[0] to -1 on steps when reneighboring occurred
+  // initialize hneigh[2] to 0 every step
+
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+
+  if (atom->nmax > nmax) {
+    nmax = atom->nmax;
+    memory->destroy(hneigh);
+    memory->create(hneigh,nmax,3,"pair:hneigh");
+    memory->destroy(newsite);
+    memory->create(newsite,nmax,3,"pair:newsite");
+  }
+  if (neighbor->ago == 0)
+    for (i = 0; i < nall; i++) hneigh[i][0] = -1;
+  for (i = 0; i < nall; i++) hneigh[i][2] = 0;
+
+  double **f = atom->f;
+  double **x = atom->x;
+  double *q = atom->q;
+  int *type = atom->type;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+  double qqrd2e = force->qqrd2e;
+  double cut_coulsqplus = (cut_coul+2.0*qdist)*(cut_coul+2.0*qdist);
+
+  int order1 = ewald_order&(1<<1);
+  int ni;
+  double  *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti;
+
+  inum = listinner->inum;
+  ilist = listinner->ilist;
+  numneigh = listinner->numneigh;
+  firstneigh = listinner->firstneigh;
+  
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    qtmp = q[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    if (itype == typeO && order1) {
+      if (hneigh[i][0] < 0) {
+        hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1);
+        hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2);
+        hneigh[i][2] = 1;
+        if (iH1 == -1 || iH2 == -1)
+          error->one(FLERR,"TIP4P hydrogen is missing");
+        if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
+          error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+        compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+      } else {
+        iH1 = hneigh[i][0];
+        iH2 = hneigh[i][1];
+        if (hneigh[i][2] == 0) {
+          hneigh[i][2] = 1;
+          compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+        }
+      }
+      x1 = newsite[i];
+    } else x1 = x[i];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    offseti = offset[itype];
+    lj1i = lj1[itype]; lj2i = lj2[itype]; lj3i = lj3[itype]; lj4i = lj4[itype];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      ni = sbmask(j);
+      factor_lj = special_lj[sbmask(j)];
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+        
+      if (rsq < cut_ljsq[itype][jtype] && rsq < cut_out_off_sq ) {  // lj
+        r2inv = 1.0/rsq;				
+	register double rn = r2inv*r2inv*r2inv;
+	if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
+	else {					// special case
+	  register double f = special_lj[ni];
+	  forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
+	 }
+
+        if (rsq > cut_out_on_sq) {                        // switching
+          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+        }
+
+        forcelj *= r2inv;
+	f[i][0] += delx*forcelj;
+	f[i][1] += dely*forcelj;
+	f[i][2] += delz*forcelj;
+	f[j][0] -= delx*forcelj;
+	f[j][1] -= dely*forcelj;
+	f[j][2] -= delz*forcelj;
+      }
+
+      
+      // adjust rsq and delxyz for off-site O charge(s)
+      // ADDITIONAL REQEUST REQUIRED HERE!!!!!
+
+      if (rsq < cut_coulsqplus && order1) {
+        if (itype == typeO || jtype == typeO) { 
+	  if (jtype == typeO) {
+            if (hneigh[j][0] < 0) {
+              hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
+              hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
+              hneigh[j][2] = 1;
+              if (jH1 == -1 || jH2 == -1)
+                error->one(FLERR,"TIP4P hydrogen is missing");
+              if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
+                error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+              compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+            } else {
+              jH1 = hneigh[j][0];
+              jH2 = hneigh[j][1];
+              if (hneigh[j][2] == 0) {
+                hneigh[j][2] = 1;
+                compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+              }
+            }
+            x2 = newsite[j];
+	  } else x2 = x[j];
+	  delx = x1[0] - x2[0];
+	  dely = x1[1] - x2[1];
+	  delz = x1[2] - x2[2];
+	  rsq = delx*delx + dely*dely + delz*delz;
+        }
+
+	// test current rsq against cutoff and compute Coulombic force
+      
+        if (rsq < cut_coulsq && rsq < cut_out_off_sq) {
+	  r2inv = 1.0 / rsq;
+          qri = qqrd2e*qtmp;
+          if (ni == 0) forcecoul = qri*q[j]*sqrt(r2inv);
+          else {
+            forcecoul = qri*q[j]*sqrt(r2inv)*special_coul[ni];
+          }
+
+          if (rsq > cut_out_on_sq) {                        // switching
+            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+          }
+
+	  cforce = forcecoul * r2inv;
+
+	  //if (evflag) ev_tally(i,j,nlocal,newton_pair,
+          //		       evdwl,0.0,cforce,delx,dely,delz);
+
+	  // if i,j are not O atoms, force is applied directly
+	  // if i or j are O atoms, force is on fictitious atom & partitioned
+	  // force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999)
+	  // f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f
+	  // preserves total force and torque on water molecule
+	  // virial = sum(r x F) where each water's atoms are near xi and xj
+	  // vlist stores 2,4,6 atoms whose forces contribute to virial
+
+	  if (itype != typeO) {
+	    f[i][0] += delx * cforce;
+	    f[i][1] += dely * cforce;
+	    f[i][2] += delz * cforce;
+
+	  } else {
+            fd[0] = delx*cforce;
+            fd[1] = dely*cforce;
+            fd[2] = delz*cforce;
+
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
+
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2];
+
+            f[i][0] += fO[0];
+            f[i][1] += fO[1];
+            f[i][2] += fO[2];
+
+            f[iH1][0] += fH[0];
+            f[iH1][1] += fH[1];
+            f[iH1][2] += fH[2];
+
+            f[iH2][0] += fH[0];
+            f[iH2][1] += fH[1];
+            f[iH2][2] += fH[2];
+          }
+
+	  if (jtype != typeO) {
+	    f[j][0] -= delx * cforce;
+	    f[j][1] -= dely * cforce;
+	    f[j][2] -= delz * cforce;
+
+	  } else {
+	    fd[0] = -delx*cforce;
+	    fd[1] = -dely*cforce;
+	    fd[2] = -delz*cforce;
+
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
+
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2]; 
+
+	    f[j][0] += fO[0];
+	    f[j][1] += fO[1];
+	    f[j][2] += fO[2];
+
+	    f[jH1][0] += fH[0];
+	    f[jH1][1] += fH[1];
+	    f[jH1][2] += fH[2];
+
+	    f[jH2][0] += fH[0];
+	    f[jH2][1] += fH[1];
+	    f[jH2][2] += fH[2];
+          }
+	}
+      }
+    }
+  }
+}
+
+/* --------------------------------------------------------------------- */
+
+void PairLJLongTIP4PLong::compute_middle()
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype,itable;
+  int iH1,iH2,jH1,jH2;
+  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
+  double r,r2inv,forcecoul,forcelj,cforce;
+  double factor_coul,factor_lj;
+  double grij,expm2,prefactor,t,erfc;
+  double xiM[3],xjM[3],fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];// f1[3];
+  double *x1,*x2;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double rsq,qri;
+
+  double cut_in_off = cut_respa[0];
+  double cut_in_on = cut_respa[1];
+  double cut_out_on = cut_respa[2];
+  double cut_out_off = cut_respa[3];
+
+  double cut_in_diff = cut_in_on - cut_in_off;
+  double cut_out_diff = cut_out_off - cut_out_on;
+  double cut_in_off_sq = cut_in_off*cut_in_off;
+  double cut_in_on_sq = cut_in_on*cut_in_on;
+  double cut_out_on_sq = cut_out_on*cut_out_on;
+  double cut_out_off_sq = cut_out_off*cut_out_off;
+
+  // reallocate hneigh & newsite if necessary
+  // initialize hneigh[0] to -1 on steps when reneighboring occurred
+  // initialize hneigh[2] to 0 every step
+
+  int nlocal = atom->nlocal;
+
+  double **f = atom->f;
+  double **x = atom->x;
+  double *q = atom->q;
+  int *type = atom->type;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+  double qqrd2e = force->qqrd2e;
+  double cut_coulsqplus = (cut_coul+2.0*qdist)*(cut_coul+2.0*qdist);
+
+  int order1 = ewald_order&(1<<1);
+  int ni;
+  double  *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti;
+  double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2;
+
+  inum = listmiddle->inum;
+  ilist = listmiddle->ilist;
+  numneigh = listmiddle->numneigh;
+  firstneigh = listmiddle->firstneigh;
+  
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    qtmp = q[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    if (itype == typeO && order1) {
+      if (hneigh[i][0] < 0) {
+        hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1);
+        hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2);
+        hneigh[i][2] = 1;
+        if (iH1 == -1 || iH2 == -1)
+          error->one(FLERR,"TIP4P hydrogen is missing");
+        if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
+          error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+        compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+      } else {
+        iH1 = hneigh[i][0];
+        iH2 = hneigh[i][1];
+        if (hneigh[i][2] == 0) {
+          hneigh[i][2] = 1;
+          compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+        }
+      }
+      x1 = newsite[i];
+    } else x1 = x[i];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    offseti = offset[itype];
+    lj1i = lj1[itype]; lj2i = lj2[itype]; lj3i = lj3[itype]; lj4i = lj4[itype];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      ni = sbmask(j);
+      factor_lj = special_lj[sbmask(j)];
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+        
+      if (rsq < cut_ljsq[itype][jtype] && rsq >= cut_in_off_sq && rsq <= cut_out_off_sq ) {  // lj
+        r2inv = 1.0/rsq;				
+	register double rn = r2inv*r2inv*r2inv;
+	if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
+	else {					// special case
+	  register double f = special_lj[ni];
+	  forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
+	 }
+
+        if (rsq < cut_in_on_sq) {                                // switching
+          register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+          forcelj  *= rsw*rsw*(3.0 - 2.0*rsw);
+        }
+        if (rsq > cut_out_on_sq) {
+          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+        }
+
+        forcelj *= r2inv;
+	f[i][0] += delx*forcelj;
+	f[i][1] += dely*forcelj;
+	f[i][2] += delz*forcelj;
+	f[j][0] -= delx*forcelj;
+	f[j][1] -= dely*forcelj;
+	f[j][2] -= delz*forcelj;
+      }
+
+      
+      // adjust rsq and delxyz for off-site O charge(s)
+      // ADDITIONAL REQEUST REQUIRED HERE!!!!!
+
+      if (rsq < cut_coulsqplus && order1) {
+        if (itype == typeO || jtype == typeO) { 
+	  if (jtype == typeO) {
+            if (hneigh[j][0] < 0) {
+              hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
+              hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
+              hneigh[j][2] = 1;
+              if (jH1 == -1 || jH2 == -1)
+                error->one(FLERR,"TIP4P hydrogen is missing");
+              if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
+                error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+              compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+            } else {
+              jH1 = hneigh[j][0];
+              jH2 = hneigh[j][1];
+              if (hneigh[j][2] == 0) {
+                hneigh[j][2] = 1;
+                compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+              }
+            }
+            x2 = newsite[j];
+	  } else x2 = x[j];
+	  delx = x1[0] - x2[0];
+	  dely = x1[1] - x2[1];
+	  delz = x1[2] - x2[2];
+	  rsq = delx*delx + dely*dely + delz*delz;
+        }
+
+	// test current rsq against cutoff and compute Coulombic force
+      
+        if (rsq < cut_coulsq &&  rsq >= cut_in_off_sq && rsq <= cut_out_off_sq) {
+	  r2inv = 1.0 / rsq;
+          qri = qqrd2e*qtmp;
+          if (ni == 0) forcecoul = qri*q[j]*sqrt(r2inv);
+          else {
+            forcecoul = qri*q[j]*sqrt(r2inv)*special_coul[ni];
+          }
+
+          if (rsq < cut_in_on_sq) {                                // switching
+            register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+            forcecoul  *= rsw*rsw*(3.0 - 2.0*rsw);
+          }
+          if (rsq > cut_out_on_sq) {
+            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
+          }
+
+	  cforce = forcecoul * r2inv;
+
+	  //if (evflag) ev_tally(i,j,nlocal,newton_pair,
+          //		       evdwl,0.0,cforce,delx,dely,delz);
+
+	  // if i,j are not O atoms, force is applied directly
+	  // if i or j are O atoms, force is on fictitious atom & partitioned
+	  // force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999)
+	  // f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f
+	  // preserves total force and torque on water molecule
+	  // virial = sum(r x F) where each water's atoms are near xi and xj
+	  // vlist stores 2,4,6 atoms whose forces contribute to virial
+
+	  if (itype != typeO) {
+	    f[i][0] += delx * cforce;
+	    f[i][1] += dely * cforce;
+	    f[i][2] += delz * cforce;
+
+	  } else {
+            fd[0] = delx*cforce;
+            fd[1] = dely*cforce;
+            fd[2] = delz*cforce;
+
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
+
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2];
+
+            f[i][0] += fO[0];
+            f[i][1] += fO[1];
+            f[i][2] += fO[2];
+
+            f[iH1][0] += fH[0];
+            f[iH1][1] += fH[1];
+            f[iH1][2] += fH[2];
+
+            f[iH2][0] += fH[0];
+            f[iH2][1] += fH[1];
+            f[iH2][2] += fH[2];
+          }
+
+	  if (jtype != typeO) {
+	    f[j][0] -= delx * cforce;
+	    f[j][1] -= dely * cforce;
+	    f[j][2] -= delz * cforce;
+
+	  } else {
+	    fd[0] = -delx*cforce;
+	    fd[1] = -dely*cforce;
+	    fd[2] = -delz*cforce;
+
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
+
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2]; 
+
+	    f[j][0] += fO[0];
+	    f[j][1] += fO[1];
+	    f[j][2] += fO[2];
+
+	    f[jH1][0] += fH[0];
+	    f[jH1][1] += fH[1];
+	    f[jH1][2] += fH[2];
+
+	    f[jH2][0] += fH[0];
+	    f[jH2][1] += fH[1];
+	    f[jH2][2] += fH[2];
+          }
+	}
+      }
+    }
+  }
+}
+
+/* --------------------------------------------------------------------- */
+
+void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype,itable;
+  int n,vlist[6];
+  int key;
+  int iH1,iH2,jH1,jH2;
+  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
+  double fraction,table;
+  double r,r2inv,forcecoul,forcelj,cforce, respa_coul, respa_lj, frespa,fvirial;
+  double factor_coul,factor_lj;
+  double grij,expm2,prefactor,t,erfc;
+  double xiM[3],xjM[3],fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];// f1[3];
+  double *x1,*x2;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double rsq,qri;
+  int respa_flag;
+ 
+  evdwl = ecoul = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  // reallocate hneigh & newsite if necessary
+  // initialize hneigh[0] to -1 on steps when reneighboring occurred
+  // initialize hneigh[2] to 0 every step
+
+  int nlocal = atom->nlocal;
+
+  double **f = atom->f;
+  double **x = atom->x;
+  double *q = atom->q;
+  int *type = atom->type;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+  double qqrd2e = force->qqrd2e;
+  double cut_coulsqplus = (cut_coul+2.0*qdist)*(cut_coul+2.0*qdist);
+
+  int order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6);
+  int ni;
+  double  *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti;
+  double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2;
+
+  double cut_in_off = cut_respa[2];
+  double cut_in_on = cut_respa[3];
+  
+  double cut_in_diff = cut_in_on - cut_in_off;
+  double cut_in_off_sq = cut_in_off*cut_in_off;
+  double cut_in_on_sq = cut_in_on*cut_in_on;
+
+  inum = listouter->inum;
+  ilist = listouter->ilist;
+  numneigh = listouter->numneigh;
+  firstneigh = listouter->firstneigh;
+  
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    qtmp = q[i];
+    qri = qtmp*qqrd2e;
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    if (itype == typeO) {
+      if (hneigh[i][0] < 0) {
+        hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1);
+        hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2);
+        hneigh[i][2] = 1;
+        if (iH1 == -1 || iH2 == -1)
+          error->one(FLERR,"TIP4P hydrogen is missing");
+        if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
+          error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+        compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+      } else {
+        iH1 = hneigh[i][0];
+        iH2 = hneigh[i][1];
+        if (hneigh[i][2] == 0) {
+          hneigh[i][2] = 1;
+          compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+        }
+      }
+      x1 = newsite[i];
+    } else x1 = x[i];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    offseti = offset[itype];
+    lj1i = lj1[itype]; lj2i = lj2[itype]; lj3i = lj3[itype]; lj4i = lj4[itype];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      ni = sbmask(j);
+      factor_lj = special_lj[sbmask(j)];
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+        
+      respa_coul = 0;
+      respa_lj = 0;
+      if (rsq < cut_ljsq[itype][jtype]) {			// lj
+        frespa = 1.0;                                       // check whether and how to compute respa corrections
+        respa_flag = rsq < cut_in_on_sq ? 1 : 0;
+        if (respa_flag && (rsq > cut_in_off_sq)) {
+          register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+          frespa = 1-rsw*rsw*(3.0-2.0*rsw);
+        }
+
+        r2inv = 1.0/rsq;
+        register double rn = r2inv*r2inv*r2inv;
+        if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
+            frespa*rn*(rn*lj1i[jtype]-lj2i[jtype]) :
+            frespa*rn*(rn*lj1i[jtype]-lj2i[jtype])*special_lj[ni];
+        if (order6) {                                        // long-range form
+          if (!ndisptablebits || rsq <= tabinnerdispsq) {
+            register double x2 = g2*rsq, a2 = 1.0/x2;
+            x2 = a2*exp(-x2)*lj4i[jtype];
+            if (ni == 0) {
+              forcelj =
+                (rn*=rn)*lj1i[jtype]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq-respa_lj;
+              if (eflag) evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
+            }
+            else {                                        // correct for special
+              register double f = special_lj[ni], t = rn*(1.0-f);
+              forcelj = f*(rn *= rn)*lj1i[jtype]-
+                g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype]-respa_lj;
+              if (eflag)
+                evdwl = f*rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[jtype];
+            }
+          }
+          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) {
+              forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]-respa_lj;
+              if (eflag) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
+            }
+            else {					// special case
+              register double f = special_lj[ni], t = rn*(1.0-f);
+              forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype]-respa_lj;
+              if (eflag) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
+            }
+          }
+        }
+        else {                                                // cut form
+          if (ni == 0) {
+            forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype])-respa_lj;
+            if (eflag) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
+          }
+          else {                                        // correct for special
+            register double f = special_lj[ni];
+            forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype])-respa_lj;
+            if (eflag)
+              evdwl = f*(rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
+          }
+        }
+
+        forcelj *= r2inv;
+	f[i][0] += delx*forcelj;
+	f[i][1] += dely*forcelj;
+	f[i][2] += delz*forcelj;
+	f[j][0] -= delx*forcelj;
+	f[j][1] -= dely*forcelj;
+	f[j][2] -= delz*forcelj;
+      
+        if (evflag) {
+          fvirial = forcelj + respa_lj*r2inv;
+          ev_tally(i,j,nlocal,newton_pair,
+	           evdwl,0.0,fvirial,delx,dely,delz);
+        }
+      }
+
+      
+      // adjust rsq and delxyz for off-site O charge(s)
+      // ADDITIONAL REQEUST REQUIRED HERE!!!!!
+
+      if (rsq < cut_coulsqplus) {
+        if (itype == typeO || jtype == typeO) { 
+	  if (jtype == typeO) {
+            if (hneigh[j][0] < 0) {
+              hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
+              hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
+              hneigh[j][2] = 1;
+              if (jH1 == -1 || jH2 == -1)
+                error->one(FLERR,"TIP4P hydrogen is missing");
+              if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
+                error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+              compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+            } else {
+              jH1 = hneigh[j][0];
+              jH2 = hneigh[j][1];
+              if (hneigh[j][2] == 0) {
+                hneigh[j][2] = 1;
+                compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+              }
+            }
+            x2 = newsite[j];
+	  } else x2 = x[j];
+	  delx = x1[0] - x2[0];
+	  dely = x1[1] - x2[1];
+	  delz = x1[2] - x2[2];
+	  rsq = delx*delx + dely*dely + delz*delz;
+        }
+
+	// test current rsq against cutoff and compute Coulombic force
+        if ((rsq < cut_coulsq) && order1) {
+
+          frespa = 1.0;                                       // check whether and how to compute respa corrections
+          respa_flag = rsq < cut_in_on_sq ? 1 : 0;
+          if (respa_flag && (rsq > cut_in_off_sq)) {
+            register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+            frespa = 1-rsw*rsw*(3.0-2.0*rsw);
+          }
+
+          r2inv = 1.0 / rsq;
+          if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
+            register double r = sqrt(rsq), 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);
+              forcecoul = (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);
+              forcecoul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r-respa_coul;
+              if (eflag) ecoul = t-r;
+            }
+          }                                                // table real space
+          else {
+            if (respa_flag) {
+              register double r = sqrt(rsq), 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 = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
+            if (ni == 0) {
+              forcecoul = 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]);
+              forcecoul = 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);
+              }
+            }
+          } 
+
+          cforce = forcecoul * r2inv;
+          fvirial = (forcecoul + respa_coul) * r2inv;
+          double fvcf = fvirial/cforce;
+
+          // if i,j are not O atoms, force is applied directly
+	  // if i or j are O atoms, force is on fictitious atom & partitioned
+	  // force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999)
+	  // f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f
+	  // preserves total force and torque on water molecule
+	  // virial = sum(r x F) where each water's atoms are near xi and xj
+	  // vlist stores 2,4,6 atoms whose forces contribute to virial
+
+	  n = 0;
+          key = 0;
+
+	  if (itype != typeO) {
+	    f[i][0] += delx * cforce;
+            f[i][1] += dely * cforce;
+	    f[i][2] += delz * cforce;
+
+            if (vflag) {
+              v[0] = x[i][0] * delx * fvirial;
+              v[1] = x[i][1] * dely * fvirial;
+              v[2] = x[i][2] * delz * fvirial;
+              v[3] = x[i][0] * dely * fvirial;
+              v[4] = x[i][0] * delz * fvirial;
+              v[5] = x[i][1] * delz * fvirial;
+            }
+          vlist[n++] = i;
+
+	  } else {
+            key += 1;
+            fd[0] = delx*cforce;
+            fd[1] = dely*cforce;
+            fd[2] = delz*cforce;
+
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
+
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2];
+
+            f[i][0] += fO[0];
+            f[i][1] += fO[1];
+            f[i][2] += fO[2];
+
+            f[iH1][0] += fH[0];
+            f[iH1][1] += fH[1];
+            f[iH1][2] += fH[2];
+
+            f[iH2][0] += fH[0];
+            f[iH2][1] += fH[1];
+            f[iH2][2] += fH[2];
+
+	    if (vflag) {
+	      domain->closest_image(x[i],x[iH1],xH1);
+	      domain->closest_image(x[i],x[iH2],xH2);
+
+	      v[0] = (x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0])*fvcf;
+	      v[1] = (x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1])*fvcf;
+	      v[2] = (x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2])*fvcf;
+	      v[3] = (x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1])*fvcf;
+	      v[4] = (x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2])*fvcf;
+	      v[5] = (x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2])*fvcf;
+	    }
+	    vlist[n++] = i;
+	    vlist[n++] = iH1;
+	    vlist[n++] = iH2;
+  	  }
+
+	  if (jtype != typeO) {
+	    f[j][0] -= delx * cforce;
+	    f[j][1] -= dely * cforce;
+	    f[j][2] -= delz * cforce;
+
+	    if (vflag) {
+	      v[0] -= x[j][0] * delx * fvirial;
+	      v[1] -= x[j][1] * dely * fvirial;
+	      v[2] -= x[j][2] * delz * fvirial;
+	      v[3] -= x[j][0] * dely * fvirial;
+	      v[4] -= x[j][0] * delz * fvirial;
+	      v[5] -= x[j][1] * delz * fvirial;
+            }
+	    vlist[n++] = j;
+
+	  } else {
+            key += 2;
+
+	    fd[0] = -delx*cforce;
+	    fd[1] = -dely*cforce;
+	    fd[2] = -delz*cforce;
+
+            fO[0] = fd[0]*(1 - alpha);
+            fO[1] = fd[1]*(1 - alpha);
+            fO[2] = fd[2]*(1 - alpha);
+
+            fH[0] = 0.5 * alpha * fd[0];
+            fH[1] = 0.5 * alpha * fd[1];
+            fH[2] = 0.5 * alpha * fd[2]; 
+
+	    f[j][0] += fO[0];
+	    f[j][1] += fO[1];
+	    f[j][2] += fO[2];
+
+	    f[jH1][0] += fH[0];
+	    f[jH1][1] += fH[1];
+	    f[jH1][2] += fH[2];
+
+	    f[jH2][0] += fH[0];
+	    f[jH2][1] += fH[1];
+	    f[jH2][2] += fH[2];
+
+	    if (vflag) {
+	      domain->closest_image(x[j],x[jH1],xH1);
+	      domain->closest_image(x[j],x[jH2],xH2);
+
+	      v[0] += (x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0])*fvcf;
+	      v[1] += (x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1])*fvcf;
+	      v[2] += (x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2])*fvcf;
+	      v[3] += (x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1])*fvcf;
+	      v[4] += (x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2])*fvcf;
+	      v[5] += (x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2])*fvcf;
+            }
+      	    vlist[n++] = j;
+	    vlist[n++] = jH1;
+	    vlist[n++] = jH2;
+	  }
+ 
+          if (evflag) ev_tally_tip4p(key,vlist,v,ecoul,alpha);
+        }
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    global settings
 ------------------------------------------------------------------------- */

src/KSPACE/pair_lj_long_tip4p_long.h

@@ -29,6 +29,9 @@ class PairLJLongTIP4PLong : public PairLJLongCoulLong {
   PairLJLongTIP4PLong(class LAMMPS *);
   ~PairLJLongTIP4PLong();
   virtual void compute(int, int);
+  virtual void compute_inner();
+  virtual void compute_middle();
+  virtual void compute_outer(int, int);
   void settings(int, char **);
   void init_style();
   double init_one(int, int);

src/KSPACE/pppm_disp.cpp

@@ -168,6 +168,8 @@ PPPMDisp::PPPMDisp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
   cg_peratom = NULL;
   cg_6 = NULL;
   cg_peratom_6 = NULL;
+
+  memset(function, 0, EWALD_FUNCS*sizeof(int));
 }
 
 /* ----------------------------------------------------------------------

src/OPT/pair_lj_long_coul_long_opt.cpp

@@ -726,7 +726,7 @@ void PairLJLongCoulLongOpt::eval_outer()
   double cut_in_off_sq = cut_in_off*cut_in_off;
   double cut_in_on_sq = cut_in_on*cut_in_on;
   
-  ineighn = (ineigh = list->ilist)+list->inum;
+  ineighn = (ineigh = listouter->ilist)+listouter->inum;
   
   for (; ineigh<ineighn; ++ineigh) {                        // loop over my atoms
     i = *ineigh; fi = f0+3*i;
@@ -735,7 +735,7 @@ void PairLJLongCoulLongOpt::eval_outer()
     lj1i = lj1[typei]; lj2i = lj2[typei]; lj3i = lj3[typei]; lj4i = lj4[typei];
     cutsqi = cutsq[typei]; cut_ljsqi = cut_ljsq[typei];
     memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
-    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
+    jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i];
     
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;

src/USER-OMP/neigh_respa_omp.cpp

@@ -100,37 +100,33 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list)
 #if defined(_OPENMP)
 #pragma omp critical
 #endif
+   {
     if (pgsize - npnt < oneatom) {
       npnt = 0;
       npage += nthreads;
-      if (npage == list->maxpage) list->add_pages(nthreads);
+      if (npage >= list->maxpage) list->add_pages(nthreads);
     }
     neighptr = &(list->pages[npage][npnt]);
     n = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (pgsize - npnt_inner < oneatom) {
       npnt_inner = 0;
       npage_inner += nthreads;
-      if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads);
+      if (npage_inner >= listinner->maxpage) listinner->add_pages(nthreads);
     }
     neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]);
     n_inner = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (respamiddle) {
       if (pgsize - npnt_middle < oneatom) {
         npnt_middle = 0;
         npage_middle += nthreads;
-        if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads);
+        if (npage_middle >= listmiddle->maxpage) listmiddle->add_pages(nthreads);
       }
       neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]);
       n_middle = 0;
     }
+   }
 
     itype = type[i];
     xtmp = x[i][0];
@@ -281,37 +277,33 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list)
 #if defined(_OPENMP)
 #pragma omp critical
 #endif
+   {
     if (pgsize - npnt < oneatom) {
       npnt = 0;
       npage += nthreads;
-      if (npage == list->maxpage) list->add_pages(nthreads);
+      if (npage >= list->maxpage) list->add_pages(nthreads);
     }
     neighptr = &(list->pages[npage][npnt]);
     n = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (pgsize - npnt_inner < oneatom) {
       npnt_inner = 0;
       npage_inner += nthreads;
-      if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads);
+      if (npage_inner >= listinner->maxpage) listinner->add_pages(nthreads);
     }
     neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]);
     n_inner = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (respamiddle) {
       if (pgsize - npnt_middle < oneatom) {
         npnt_middle = 0;
         npage_middle += nthreads;
-        if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads);
+        if (npage_middle >= listmiddle->maxpage) listmiddle->add_pages(nthreads);
       }
       neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]);
       n_middle = 0;
     }
+   }
 
     itag = tag[i];
     itype = type[i];
@@ -484,37 +476,33 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list)
 #if defined(_OPENMP)
 #pragma omp critical
 #endif
+   {
     if (pgsize - npnt < oneatom) {
       npnt = 0;
       npage += nthreads;
-      if (npage == list->maxpage) list->add_pages(nthreads);
+      if (npage >= list->maxpage) list->add_pages(nthreads);
     }
     neighptr = &(list->pages[npage][npnt]);
     n = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (pgsize - npnt_inner < oneatom) {
       npnt_inner = 0;
       npage_inner += nthreads;
-      if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads);
+      if (npage_inner >= listinner->maxpage) listinner->add_pages(nthreads);
     }
     neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]);
     n_inner = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (respamiddle) {
       if (pgsize - npnt_middle < oneatom) {
         npnt_middle = 0;
         npage_middle += nthreads;
-        if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads);
+        if (npage_middle >= listmiddle->maxpage) listmiddle->add_pages(nthreads);
       }
       neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]);
       n_middle = 0;
     }
+   }
 
     itype = type[i];
     xtmp = x[i][0];
@@ -677,37 +665,33 @@ void Neighbor::respa_bin_newton_omp(NeighList *list)
 #if defined(_OPENMP)
 #pragma omp critical
 #endif
+   {
     if (pgsize - npnt < oneatom) {
       npnt = 0;
       npage += nthreads;
-      if (npage == list->maxpage) list->add_pages(nthreads);
+      if (npage >= list->maxpage) list->add_pages(nthreads);
     }
     neighptr = &(list->pages[npage][npnt]);
     n = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (pgsize - npnt_inner < oneatom) {
       npnt_inner = 0;
       npage_inner += nthreads;
-      if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads);
+      if (npage_inner >= listinner->maxpage) listinner->add_pages(nthreads);
     }
     neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]);
     n_inner = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (respamiddle) {
       if (pgsize - npnt_middle < oneatom) {
         npnt_middle = 0;
         npage_middle += nthreads;
-        if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads);
+        if (npage_middle >= listmiddle->maxpage) listmiddle->add_pages(nthreads);
       }
       neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]);
       n_middle = 0;
     }
+   }
 
     itype = type[i];
     xtmp = x[i][0];
@@ -911,37 +895,33 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list)
 #if defined(_OPENMP)
 #pragma omp critical
 #endif
+   {
     if (pgsize - npnt < oneatom) {
       npnt = 0;
       npage += nthreads;
-      if (npage == list->maxpage) list->add_pages(nthreads);
+      if (npage >= list->maxpage) list->add_pages(nthreads);
     }
     neighptr = &(list->pages[npage][npnt]);
     n = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (pgsize - npnt_inner < oneatom) {
       npnt_inner = 0;
       npage_inner += nthreads;
-      if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads);
+      if (npage_inner >= listinner->maxpage) listinner->add_pages(nthreads);
     }
     neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]);
     n_inner = 0;
 
-#if defined(_OPENMP)
-#pragma omp critical
-#endif
     if (respamiddle) {
       if (pgsize - npnt_middle < oneatom) {
         npnt_middle = 0;
         npage_middle += nthreads;
-        if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads);
+        if (npage_middle >= listmiddle->maxpage) listmiddle->add_pages(nthreads);
       }
       neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]);
       n_middle = 0;
     }
+   }
 
     itype = type[i];
     xtmp = x[i][0];

src/USER-OMP/pair_buck_long_coul_long_omp.cpp

@@ -317,7 +317,7 @@ void PairBuckLongCoulLongOMP::compute_inner()
 
   const int nall = atom->nlocal + atom->nghost;
   const int nthreads = comm->nthreads;
-  const int inum = list->inum;
+  const int inum = listinner->inum;
 #if defined(_OPENMP)
 #pragma omp parallel default(none)
 #endif
@@ -339,7 +339,7 @@ void PairBuckLongCoulLongOMP::compute_middle()
 
   const int nall = atom->nlocal + atom->nghost;
   const int nthreads = comm->nthreads;
-  const int inum = list->inum;
+  const int inum = listmiddle->inum;
 
 #if defined(_OPENMP)
 #pragma omp parallel default(none)
@@ -367,7 +367,7 @@ void PairBuckLongCoulLongOMP::compute_outer(int eflag, int vflag)
 
   const int nall = atom->nlocal + atom->nghost;
   const int nthreads = comm->nthreads;
-  const int inum = list->inum;
+  const int inum = listouter->inum;
 
 #if defined(_OPENMP)
 #pragma omp parallel default(none) shared(eflag,vflag)
@@ -803,7 +803,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
   const double *x0 = x[0];
   double *f0 = f[0], *fi = 0;
 
-  int *ilist = list->ilist;
+  int *ilist = listinner->ilist;
 
   const int newton_pair = force->newton_pair;
 
@@ -827,7 +827,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
     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];
+    jneighn = (jneigh = listinner->firstneigh[i])+listinner->numneigh[i];
 
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;
@@ -896,7 +896,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
   const double *x0 = x[0];
   double *f0 = f[0], *fi = 0;
 
-  int *ilist = list->ilist;
+  int *ilist = listmiddle->ilist;
 
   const int newton_pair = force->newton_pair;
 
@@ -924,7 +924,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
     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];
+    jneighn = (jneigh = listmiddle->firstneigh[i])+listmiddle->numneigh[i];
 
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;
@@ -1001,7 +1001,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
   const double *x0 = x[0];
   double *f0 = f[0], *fi = f0;
 
-  int *ilist = list->ilist;
+  int *ilist = listouter->ilist;
   
   int i, j, ii;
   int *jneigh, *jneighn, typei, typej, ni, respa_flag;
@@ -1027,7 +1027,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
     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];
+    jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i];
 
     for (; jneigh<jneighn; ++jneigh) {                        // loop over neighbors
       j = *jneigh;

src/USER-OMP/pair_lj_long_coul_long_omp.h

@@ -35,11 +35,28 @@ class PairLJLongCoulLongOMP : public PairLJLongCoulLong, public ThrOMP {
   PairLJLongCoulLongOMP(class LAMMPS *);
 
   virtual void compute(int, int);
+  virtual void compute_inner();
+  virtual void compute_middle();
+  virtual void compute_outer(int, int);
   virtual double memory_usage();
 
  private:
-  template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
+  template <const int EVFLAG, const int EFLAG,
-  void eval(int ifrom, int ito, ThrData * const thr);
+    const int NEWTON_PAIR, const int CTABLE, const int LJTABLE,
+    const int ORDER1, const int ORDER6 >
+  void eval(int, int, ThrData * const);
+
+  template <const int EVFLAG, const int EFLAG,
+    const int NEWTON_PAIR, const int CTABLE, const int LJTABLE,
+    const int ORDER1, const int ORDER6 >
+  void eval_outer(int, int, ThrData * const);
+
+
+  void eval_inner(int, int, ThrData *const);
+  void eval_middle(int, int, ThrData *const);
+
+  
+
 };
 
 }