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

src/GPU/pppm_gpu.cpp

@@ -284,21 +284,24 @@ void PPPMGPU::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum energy across procs and add in volume-dependent term
-  // reset qsum and qsqsum if atom count has changed
 
   if (eflag_global) {
     double energy_all;
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (atom->natoms != natoms_original) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     energy *= 0.5*volume;
-    energy -= g_ewald*qsqsum/1.772453851 +
+    energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qscale;
   }

src/KSPACE/ewald.cpp

@@ -431,19 +431,22 @@ void Ewald::compute(int eflag, int vflag)
     if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2];
   }
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum global energy across Kspace vevs and add in volume-dependent term
-  // reset qsum and qsqsum if atom count has changed
 
   if (eflag_global) {
     for (k = 0; k < kcount; k++)
       energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] +
                          sfacim_all[k]*sfacim_all[k]);
 
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qscale;

src/KSPACE/ewald_disp.cpp

@@ -246,6 +246,9 @@ void EwaldDisp::setup()
     error->all(FLERR,"KSpace accuracy too low");
   }
 
+  if (qsum_update_flag)
+    error->all(FLERR,"Kspace style ewald/disp does not support dynamic charges");
+
   bigint natoms = atom->natoms;
   double err;
   int kxmax = 1;

src/KSPACE/msm.cpp

@@ -566,8 +566,16 @@ void MSM::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum global energy across procs and add in self-energy term
-  // reset qsum and qsqsum if atom count has changed
 
   const double qscale = qqrd2e * scale;
 
@@ -576,11 +584,6 @@ void MSM::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     double e_self = qsqsum*gamma(0.0)/cutoff;
     energy -= e_self;
     energy *= 0.5*qscale;

src/KSPACE/msm_cg.cpp

@@ -248,7 +248,16 @@ void MSMCG::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // total long-range energy
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
+  // sum global energy across procs and add in self-energy term
 
   const double qscale = force->qqrd2e * scale;
 
@@ -257,7 +266,7 @@ void MSMCG::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    double e_self = qsqsum*gamma(0.0)/cutoff;  // Self-energy term
+    double e_self = qsqsum*gamma(0.0)/cutoff;
     energy -= e_self;
     energy *= 0.5*qscale;
   }

src/KSPACE/pppm.cpp

@@ -667,8 +667,16 @@ void PPPM::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum global energy across procs and add in volume-dependent term
-  // reset qsum and qsqsum if atom count has changed
 
   const double qscale = qqrd2e * scale;
 
@@ -677,11 +685,6 @@ void PPPM::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);

src/KSPACE/pppm_cg.cpp

@@ -208,8 +208,16 @@ void PPPMCG::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum global energy across procs and add in volume-dependent term
-  // reset qsum and qsqsum if atom count has changed
 
   const double qscale = qqrd2e * scale;
 
@@ -218,11 +226,6 @@ void PPPMCG::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);

src/KSPACE/pppm_disp.cpp

@@ -1121,8 +1121,16 @@ void PPPMDisp::compute(int eflag, int vflag)
     if (evflag_atom) fieldforce_none_peratom();
   }
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum energy across procs and add in volume-dependent term
-  // reset qsum and qsqsum if atom count has changed
 
   const double qscale = force->qqrd2e * scale;
   if (eflag_global) {
@@ -1135,11 +1143,6 @@ void PPPMDisp::compute(int eflag, int vflag)
     energy_1 *= 0.5*volume;
     energy_6 *= 0.5*volume;
     
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     energy_1 -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij +

src/KSPACE/pppm_stagger.cpp

@@ -202,8 +202,16 @@ void PPPMStagger::compute(int eflag, int vflag)
     stagger += 1.0/float(nstagger);
   }
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum global energy across procs and add in volume-dependent term
-  // reset qsum and qsqsum if atom count has changed
 
   const double qscale = qqrd2e * scale;
 
@@ -212,11 +220,6 @@ void PPPMStagger::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (qsum_update_flag || (atom->natoms != natoms_original)) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     energy *= 0.5*volume/float(nstagger);
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);

src/USER-CUDA/pppm_cuda.cpp

@@ -690,17 +690,18 @@ void PPPMCuda::compute(int eflag, int vflag)
   }
 
   // extend size of per-atom arrays if necessary
+  // and force update of device data, if charges can change
 
-  if ((cu_atom->update_nmax)||(old_nmax==0)) {
+  if ((cu_atom->update_nmax)||(old_nmax==0)||qsum_update_flag) {
     memory->destroy(part2grid);
     nmax = atom->nmax;
     memory->create(part2grid,nmax,3,"pppm:part2grid");
-         delete cu_part2grid;
-         delete [] adev_data_array;
-         adev_data_array=new dev_array[1];
-         cu_part2grid = new cCudaData<int  , int   , yx > ((int*)part2grid,adev_data_array, nmax,3);
+    delete cu_part2grid;
+    delete [] adev_data_array;
+    adev_data_array=new dev_array[1];
+    cu_part2grid = new cCudaData<int  , int   , yx > ((int*)part2grid,adev_data_array, nmax,3);
 
-          pppm_device_update(&cuda->shared_data,cu_part2grid->dev_data(),atom->nlocal,atom->nmax);
+    pppm_device_update(&cuda->shared_data,cu_part2grid->dev_data(),atom->nlocal,atom->nmax);
     old_nmax=nmax;
   }
   if(cu_atom->update_nlocal) {pppm_update_nlocal(cu_atom->nlocal);}

src/USER-CUDA/pppm_old.cpp

@@ -724,8 +724,16 @@ void PPPMOld::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // sum global energy across procs and add in volume-dependent term
-  // reset qsum and qsqsum if atom count has changed
 
   const double qscale = qqrd2e * scale;
 
@@ -734,11 +742,6 @@ void PPPMOld::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (atom->natoms != natoms_original) {
-      qsum_qsq(0);
-      natoms_original = atom->natoms;
-    }
-
     energy *= 0.5*volume;
     energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);

src/USER-OMP/ewald_omp.cpp

@@ -86,6 +86,16 @@ void EwaldOMP::compute(int eflag, int vflag)
   MPI_Allreduce(sfacrl,sfacrl_all,kcount,MPI_DOUBLE,MPI_SUM,world);
   MPI_Allreduce(sfacim,sfacim_all,kcount,MPI_DOUBLE,MPI_SUM,world);
 
+  // update qsum and qsqsum, if needed
+  // (n.b. needs to be done outside of the multi-threaded region)
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
   // K-space portion of electric field
   // double loop over K-vectors and local atoms
 

src/USER-OMP/msm_cg_omp.cpp

@@ -250,7 +250,16 @@ void MSMCGOMP::compute(int eflag, int vflag)
 
   if (evflag_atom) fieldforce_peratom();
 
-  // total long-range energy
+  // update qsum and qsqsum, if needed
+
+  if (eflag_global || eflag_atom) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
+      qsum_qsq(0);
+      natoms_original = atom->natoms;
+    }
+  }
+
+  // sum global energy across procs and add in self-energy term
 
   const double qscale = force->qqrd2e * scale;
 
@@ -259,7 +268,7 @@ void MSMCGOMP::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    double e_self = qsqsum*gamma(0.0)/cutoff;  // Self-energy term
+    double e_self = qsqsum*gamma(0.0)/cutoff;
     energy -= e_self;
     energy *= 0.5*qscale;
   }

src/force.cpp

@@ -906,8 +906,7 @@ FILE *Force::open_potential(const char *name)
   fp = fopen(name,"r");
   if (fp) {
     potential_date(fp,name);
-    fclose(fp);
-    fp = fopen(name,"r");
+    rewind(fp);
     return fp;
   }
 
@@ -934,9 +933,10 @@ FILE *Force::open_potential(const char *name)
   strcat(newpath,pot);
 
   fp = fopen(newpath,"r");
-  potential_date(fp,name);
-  fclose(fp);
-  fp = fopen(newpath,"r");
+  if (fp) {
+    potential_date(fp,name);
+    rewind(fp);
+  }
 
   delete [] newpath;
   return fp;

src/kspace.cpp

@@ -264,17 +264,20 @@ void KSpace::ev_setup(int eflag, int vflag)
 
 void KSpace::qsum_qsq(int flag)
 {
-  qsum = qsqsum = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) {
-    qsum += atom->q[i];
-    qsqsum += atom->q[i]*atom->q[i];
+  const double * const q = atom->q;
+  const int nlocal = atom->nlocal;
+  double qsum_local(0.0), qsqsum_local(0.0);
+
+#if defined(_OPENMP)
+#pragma omp parallel for default(none) reduction(+:qsum_local,qsqsum_local)
+#endif
+  for (int i = 0; i < nlocal; i++) {
+    qsum_local += q[i];
+    qsqsum_local += q[i]*q[i];
   }
 
-  double tmp;
-  MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsum = tmp;
-  MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  qsqsum = tmp;
+  MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
 
   if (qsqsum == 0.0)
     error->all(FLERR,"Cannot use kspace solver on system with no charge");

src/pair_hybrid.cpp

@@ -745,7 +745,7 @@ void *PairHybrid::extract(const char *str, int &dim)
 {
   void *cutptr = NULL;
   void *ptr;
-  double cutvalue;
+  double cutvalue = 0.0;
 
   for (int m = 0; m < nstyles; m++) {
     ptr = styles[m]->extract(str,dim);