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

src/KSPACE/pppm.cpp

@@ -700,10 +700,13 @@ void PPPM::compute(int eflag, int vflag)
 
   // per-atom energy/virial
   // energy includes self-energy correction
+  // notal accounts for TIP4P tallying eatom/vatom for ghost atoms
 
   if (evflag_atom) {
     double *q = atom->q;
     int nlocal = atom->nlocal;
+    int ntotal = nlocal;
+    if (tip4pflag) ntotal += atom->nghost;
 
     if (eflag_atom) {
       for (i = 0; i < nlocal; i++) {
@@ -712,10 +715,11 @@ void PPPM::compute(int eflag, int vflag)
           (g_ewald*g_ewald*volume);
         eatom[i] *= qscale;
       }
+      for (i = nlocal; i < ntotal; i++) eatom[i] *= 0.5*qscale;
     }
 
     if (vflag_atom) {
-      for (i = 0; i < nlocal; i++)
+      for (i = 0; i < ntotal; i++)
         for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*qscale;
     }
   }

src/KSPACE/pppm_disp_tip4p.cpp

@@ -443,43 +443,44 @@ void PPPMDispTIP4P::fieldforce_c_peratom()
       }
     }
 
+    const double qfactor = 0.5*force->qqrd2e * scale * q[i];
 
     if (eflag_atom) {
       if (type[i] != typeO) {
-        eatom[i] += q[i]*u_pa;
+        eatom[i] += qfactor*u_pa;
       } else {
-        eatom[i] += q[i]*u_pa*(1-alpha);
+        eatom[i] += qfactor*u_pa*(1-alpha);
-        eatom[iH1] += q[i]*u_pa*alpha*0.5;
+        eatom[iH1] += qfactor*u_pa*alpha*0.5;
-        eatom[iH2] += q[i]*u_pa*alpha*0.5;
+        eatom[iH2] += qfactor*u_pa*alpha*0.5;
       }
     }
     if (vflag_atom) {
       if (type[i] != typeO) {
-        vatom[i][0] += v0*q[i];
+        vatom[i][0] += v0*qfactor;
-        vatom[i][1] += v1*q[i];
+        vatom[i][1] += v1*qfactor;
-        vatom[i][2] += v2*q[i];
+        vatom[i][2] += v2*qfactor;
-        vatom[i][3] += v3*q[i];
+        vatom[i][3] += v3*qfactor;
-        vatom[i][4] += v4*q[i];
+        vatom[i][4] += v4*qfactor;
-        vatom[i][5] += v5*q[i];
+        vatom[i][5] += v5*qfactor;
       } else {
-        vatom[i][0] += v0*(1-alpha)*q[i];
+        vatom[i][0] += v0*(1-alpha)*qfactor;
-        vatom[i][1] += v1*(1-alpha)*q[i];
+        vatom[i][1] += v1*(1-alpha)*qfactor;
-        vatom[i][2] += v2*(1-alpha)*q[i];
+        vatom[i][2] += v2*(1-alpha)*qfactor;
-        vatom[i][3] += v3*(1-alpha)*q[i];
+        vatom[i][3] += v3*(1-alpha)*qfactor;
-        vatom[i][4] += v4*(1-alpha)*q[i];
+        vatom[i][4] += v4*(1-alpha)*qfactor;
-        vatom[i][5] += v5*(1-alpha)*q[i];
+        vatom[i][5] += v5*(1-alpha)*qfactor;
-        vatom[iH1][0] += v0*alpha*0.5*q[i];
+        vatom[iH1][0] += v0*alpha*0.5*qfactor;
-        vatom[iH1][1] += v1*alpha*0.5*q[i];
+        vatom[iH1][1] += v1*alpha*0.5*qfactor;
-        vatom[iH1][2] += v2*alpha*0.5*q[i];
+        vatom[iH1][2] += v2*alpha*0.5*qfactor;
-        vatom[iH1][3] += v3*alpha*0.5*q[i];
+        vatom[iH1][3] += v3*alpha*0.5*qfactor;
-        vatom[iH1][4] += v4*alpha*0.5*q[i];
+        vatom[iH1][4] += v4*alpha*0.5*qfactor;
-        vatom[iH1][5] += v5*alpha*0.5*q[i];
+        vatom[iH1][5] += v5*alpha*0.5*qfactor;
-        vatom[iH2][0] += v0*alpha*0.5*q[i];
+        vatom[iH2][0] += v0*alpha*0.5*qfactor;
-        vatom[iH2][1] += v1*alpha*0.5*q[i];
+        vatom[iH2][1] += v1*alpha*0.5*qfactor;
-        vatom[iH2][2] += v2*alpha*0.5*q[i];
+        vatom[iH2][2] += v2*alpha*0.5*qfactor;
-        vatom[iH2][3] += v3*alpha*0.5*q[i];
+        vatom[iH2][3] += v3*alpha*0.5*qfactor;
-        vatom[iH2][4] += v4*alpha*0.5*q[i];
+        vatom[iH2][4] += v4*alpha*0.5*qfactor;
-        vatom[iH2][5] += v5*alpha*0.5*q[i];
+        vatom[iH2][5] += v5*alpha*0.5*qfactor;
       }
     }
   }

src/compute_pe_atom.cpp

@@ -97,14 +97,17 @@ void ComputePEAtom::compute_peratom()
   //   b/c some bonds/dihedrals call pair::ev_tally with pairwise info
   // nbond includes ghosts if newton_bond is set
   // ntotal includes ghosts if either newton flag is set
+  // KSpace includes ghosts if tip4pflag is set
 
   int nlocal = atom->nlocal;
   int npair = nlocal;
   int nbond = nlocal;
   int ntotal = nlocal;
+  int nkspace = nlocal;
   if (force->newton) npair += atom->nghost;
   if (force->newton_bond) nbond += atom->nghost;
   if (force->newton) ntotal += atom->nghost;
+  if (force->kspace->tip4pflag) nkspace += atom->nghost;
 
   // clear local energy array
 
@@ -137,17 +140,16 @@ void ComputePEAtom::compute_peratom()
     for (i = 0; i < nbond; i++) energy[i] += eatom[i];
   }
 
-  // communicate ghost energy between neighbor procs
-
-  if (force->newton) comm->reverse_comm_compute(this);
-
-  // KSpace contribution is already per local atom
-
   if (kspaceflag && force->kspace) {
     double *eatom = force->kspace->eatom;
-    for (i = 0; i < nlocal; i++) energy[i] += eatom[i];
+    for (i = 0; i < nkspace; i++) energy[i] += eatom[i];
   }
 
+  // communicate ghost energy between neighbor procs
+
+  if (force->newton || force->kspace->tip4pflag) 
+    comm->reverse_comm_compute(this);
+
   // zero energy of atoms not in group
   // only do this after comm since ghost contributions must be included
 

src/compute_stress_atom.cpp

@@ -111,14 +111,17 @@ void ComputeStressAtom::compute_peratom()
   //   b/c some bonds/dihedrals call pair::ev_tally with pairwise info
   // nbond includes ghosts if newton_bond is set
   // ntotal includes ghosts if either newton flag is set
+  // KSpace includes ghosts if tip4pflag is set
 
   int nlocal = atom->nlocal;
   int npair = nlocal;
   int nbond = nlocal;
   int ntotal = nlocal;
+  int nkspace = nlocal;
   if (force->newton) npair += atom->nghost;
   if (force->newton_bond) nbond += atom->nghost;
   if (force->newton) ntotal += atom->nghost;
+  if (force->kspace->tip4pflag) nkspace += atom->nghost;
 
   // clear local stress array
 
@@ -163,6 +166,13 @@ void ComputeStressAtom::compute_peratom()
         stress[i][j] += vatom[i][j];
   }
 
+  if (kspaceflag && force->kspace) {
+    double **vatom = force->kspace->vatom;
+    for (i = 0; i < nkspace; i++)
+      for (j = 0; j < 6; j++)
+        stress[i][j] += vatom[i][j];
+  }
+
   // add in per-atom contributions from relevant fixes
   // skip if vatom = NULL
   // possible during setup phase if fix has not initialized its vatom yet
@@ -180,18 +190,10 @@ void ComputeStressAtom::compute_peratom()
       }
   }
 
-  // communicate ghost atom virials between neighbor procs
+  // communicate ghost virials between neighbor procs
 
-  if (force->newton) comm->reverse_comm_compute(this);
+  if (force->newton || force->kspace->tip4pflag) 
-
+    comm->reverse_comm_compute(this);
-  // KSpace contribution is already per local atom
-
-  if (kspaceflag && force->kspace) {
-    double **vatom = force->kspace->vatom;
-    for (i = 0; i < nlocal; i++)
-      for (j = 0; j < 6; j++)
-        stress[i][j] += vatom[i][j];
-  }
 
   // zero virial of atoms not in group
   // only do this after comm since ghost contributions must be included

src/force.cpp

@@ -669,7 +669,6 @@ void Force::bounds(char *str, int nmax, int &nlo, int &nhi)
 /* ----------------------------------------------------------------------
    read a floating point value from a string
    generate an error if not a legitimate floating point value
-   called by force fields to check validity of their arguments
 ------------------------------------------------------------------------- */
 
 double Force::numeric(char *str)
@@ -689,7 +688,6 @@ double Force::numeric(char *str)
 /* ----------------------------------------------------------------------
    read an integer value from a string
    generate an error if not a legitimate integer value
-   called by force fields to check validity of their arguments
 ------------------------------------------------------------------------- */
 
 int Force::inumeric(char *str)

src/kspace.cpp

@@ -188,12 +188,12 @@ void KSpace::ev_setup(int eflag, int vflag)
 
   // reallocate per-atom arrays if necessary
 
-  if (eflag_atom && atom->nlocal > maxeatom) {
+  if (eflag_atom && atom->nmax > maxeatom) {
     maxeatom = atom->nmax;
     memory->destroy(eatom);
     memory->create(eatom,maxeatom,"kspace:eatom");
   }
-  if (vflag_atom && atom->nlocal > maxvatom) {
+  if (vflag_atom && atom->nmax > maxvatom) {
     maxvatom = atom->nmax;
     memory->destroy(vatom);
     memory->create(vatom,maxvatom,6,"kspace:vatom");
@@ -205,10 +205,12 @@ void KSpace::ev_setup(int eflag, int vflag)
   if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
   if (eflag_atom) {
     n = atom->nlocal;
+    if (tip4pflag) n += atom->nghost;
     for (i = 0; i < n; i++) eatom[i] = 0.0;
   }
   if (vflag_atom) {
     n = atom->nlocal;
+    if (tip4pflag) n += atom->nghost;
     for (i = 0; i < n; i++) {
       vatom[i][0] = 0.0;
       vatom[i][1] = 0.0;

src/lattice.cpp

@@ -18,6 +18,7 @@
 #include "update.h"
 #include "domain.h"
 #include "comm.h"
+#include "force.h"
 #include "memory.h"
 #include "error.h"
 
@@ -161,9 +162,9 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
       if (dim == 0) orient = orientx;
       else if (dim == 1) orient = orienty;
       else if (dim == 2) orient = orientz;
-      orient[0] = atoi(arg[iarg+2]);
+      orient[0] = force->inumeric(arg[iarg+2]);
-      orient[1] = atoi(arg[iarg+3]);
+      orient[1] = force->inumeric(arg[iarg+3]);
-      orient[2] = atoi(arg[iarg+4]);
+      orient[2] = force->inumeric(arg[iarg+4]);
       iarg += 5;
 
     } else if (strcmp(arg[iarg],"spacing") == 0) {