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

src/CORESHELL/pair_lj_cut_coul_long_cs.cpp

@@ -556,6 +556,7 @@ void PairLJCutCoulLongCS::compute_outer(int eflag, int vflag)
             r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
           } else if (rsq <= cut_in_on_sq)
+            r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
 
           fpair = (forcecoul + factor_lj*forcelj) * r2inv;

src/KSPACE/fix_tune_kspace.cpp

@@ -212,6 +212,7 @@ void FixTuneKspace::store_old_kspace_settings()
   old_differentiation_flag = force->kspace->differentiation_flag;
   old_slabflag = force->kspace->slabflag;
   old_slab_volfactor = force->kspace->slab_volfactor;
+  delete[] old_kspace_style;
 }
 
 /* ----------------------------------------------------------------------

src/KSPACE/msm.cpp

@@ -1712,7 +1712,7 @@ void MSM::direct(int n)
 
         egridn[icz][icy][icx] += esum;
 
-        if (vflag_atom) {
+        if (vflag_atom && !scalar_pressure_flag) {
           v0gridn[icz][icy][icx] += v0sum;
           v1gridn[icz][icy][icx] += v1sum;
           v2gridn[icz][icy][icx] += v2sum;
@@ -2075,7 +2075,7 @@ void MSM::direct_top(int n)
 
         egridn[icz][icy][icx] += esum;
 
-        if (vflag_atom) {
+        if (vflag_atom && !scalar_pressure_flag) {
           v0gridn[icz][icy][icx] += v0sum;
           v1gridn[icz][icy][icx] += v1sum;
           v2gridn[icz][icy][icx] += v2sum;
@@ -2912,7 +2912,7 @@ void MSM::compute_phis_and_dphis(const double &dx, const double &dy,
 
 inline double MSM::compute_phi(const double &xi)
 {
-  double phi;
+  double phi = 0.0;
   double abs_xi = fabs(xi);
   double xi2 = xi*xi;
 
@@ -2991,7 +2991,7 @@ inline double MSM::compute_phi(const double &xi)
 
 inline double MSM::compute_dphi(const double &xi)
 {
-  double dphi;
+  double dphi = 0.0;
   double abs_xi = fabs(xi);
 
   if (order == 4) {

src/KSPACE/pair_born_coul_msm.cpp

@@ -126,9 +126,9 @@ void PairBornCoulMSM::compute(int eflag, int vflag)
 
       if (rsq < cutsq[itype][jtype]) {
         r2inv = 1.0/rsq;
+        r = sqrt(rsq);
 
         if (rsq < cut_coulsq) {
-          r = sqrt(rsq);
           prefactor = qqrd2e * qtmp*q[j]/r;
           egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
           fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);

src/KSPACE/pair_buck_coul_msm.cpp

@@ -123,9 +123,9 @@ void PairBuckCoulMSM::compute(int eflag, int vflag)
 
       if (rsq < cutsq[itype][jtype]) {
         r2inv = 1.0/rsq;
+        r = sqrt(rsq);
 
         if (rsq < cut_coulsq) {
-          r = sqrt(rsq);
           prefactor = qqrd2e * qtmp*q[j]/r;
           egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
           fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);

src/KSPACE/pair_lj_charmm_coul_long.cpp

@@ -597,6 +597,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
               forcelj = forcelj*switch1 + philj*switch2;
             }
           } else if (rsq <= cut_in_on_sq) {
+            r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
             if (rsq > cut_lj_innersq) {
               switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *

src/KSPACE/pair_lj_charmm_coul_msm.cpp

@@ -437,6 +437,7 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
               forcelj = forcelj*switch1 + philj*switch2;
             }
           } else if (rsq <= cut_in_on_sq) {
+            r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
             if (rsq > cut_lj_innersq) {
               switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *

src/KSPACE/pair_lj_cut_coul_long.cpp

@@ -551,6 +551,7 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
             r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
           } else if (rsq <= cut_in_on_sq)
+            r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
 
           fpair = (forcecoul + factor_lj*forcelj) * r2inv;

src/KSPACE/pair_lj_cut_coul_msm.cpp

@@ -404,6 +404,7 @@ void PairLJCutCoulMSM::compute_outer(int eflag, int vflag)
             r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
           } else if (rsq <= cut_in_on_sq)
+            r6inv = r2inv*r2inv*r2inv;
             forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
 
           fpair = (forcecoul + factor_lj*forcelj) * r2inv;

src/KSPACE/pppm_disp.cpp

@@ -3543,7 +3543,7 @@ void PPPMDisp::set_init_g6()
   double acc_rspace = accuracy;
   if (accuracy_real_6 > 0) acc_rspace = accuracy_real_6;
 
-  g_ewald_6 = 1.0/cutoff_lj;
+  g_ewald_old = g_ewald_6 = 1.0/cutoff_lj;
   df_real = lj_rspace_error() - acc_rspace;
   int counter = 0;
   if (df_real > 0) {

src/REPLICA/fix_neb.cpp

@@ -138,7 +138,7 @@ void FixNEB::min_post_force(int vflag)
   // veng = PE of this replica
   // vprev,vnext = PEs of adjacent replicas
 
-  veng = pe->compute_scalar();
+  vprev = vnext = veng = pe->compute_scalar();
 
   if (ireplica < nreplica-1) MPI_Send(&veng,1,MPI_DOUBLE,procnext,0,uworld);
   if (ireplica > 0) MPI_Recv(&vprev,1,MPI_DOUBLE,procprev,0,uworld,MPI_STATUS_IGNORE);

src/SHOCK/fix_append_atoms.cpp

@@ -50,6 +50,7 @@ FixAppendAtoms::FixAppendAtoms(LAMMPS *lmp, int narg, char **arg) :
 
   scaleflag = 1;
   spatflag=0;
+  spatialid = NULL;
   xloflag = xhiflag = yloflag = yhiflag = zloflag = zhiflag = 0;
 
   tempflag = 0;
@@ -210,6 +211,7 @@ FixAppendAtoms::~FixAppendAtoms()
   delete [] basistype;
 
   if (ranflag) delete randomx;
+  if (spatflag) delete[] spatialid;
   if (tempflag) {
     delete randomt;
     delete [] gfactor1;

src/USER-DRUDE/pair_thole.cpp

@@ -392,21 +392,18 @@ double PairThole::single(int i, int j, int itype, int jtype,
   }
 
   r2inv = 1.0/rsq;
+  fforce = phicoul = 0.0;
   if (rsq < cutsq[itype][jtype]) {
     rinv = sqrt(r2inv);
     r = sqrt(rsq);
     asr = ascreen[itype][jtype] * r;
     exp_asr = exp(-asr);
-    dcoul = force->qqrd2e * qi * qj * rinv;
+    dcoul = force->qqrd2e * qi * qj * scale[itype][jtype] * rinv;
     factor_f = 0.5*(2. + (exp_asr * (-2. - asr * (2. + asr)))) - factor_coul;
-    forcecoul += factor_f * dcoul;
+    fforce = factor_f * dcoul * r2inv;
     factor_e = 0.5*(2. - (exp_asr * (2. + asr))) - factor_coul; 
-  } else forcecoul= 0.0;
-  fforce = factor_f*forcecoul * r2inv;
-
-  if (rsq < cutsq[itype][jtype]) {
     phicoul = factor_e * dcoul;
-  } else phicoul = 0.0;
+  }
 
   return phicoul;
 }

src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp

@@ -565,6 +565,8 @@ void PairLJCharmmCoulLongSoft::compute_outer(int eflag, int vflag)
               forcelj = forcelj*switch1 + philj*switch2;
             }
           } else if (rsq <= cut_in_on_sq) {
+            r4sig6 = rsq*rsq / lj2[itype][jtype];
+            denlj = lj3[itype][jtype] + rsq*r4sig6;
             forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * 
               (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
             if (rsq > cut_lj_innersq) {

src/USER-FEP/pair_lj_cut_coul_long_soft.cpp

@@ -516,6 +516,8 @@ void PairLJCutCoulLongSoft::compute_outer(int eflag, int vflag)
             forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * 
               (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
           } else if (rsq < cut_in_on_sq) {
+            r4sig6 = rsq*rsq / lj2[itype][jtype];
+            denlj = lj3[itype][jtype] + rsq*r4sig6;
             forcelj = lj1[itype][jtype] * epsilon[itype][jtype] * 
               (48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
           }

src/USER-FEP/pair_lj_cut_soft.cpp

@@ -393,6 +393,7 @@ void PairLJCutSoft::compute_outer(int eflag, int vflag)
         }
 
         if (eflag) {
+          denlj = lj3[itype][jtype] + rsq*r4sig6;
           evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] * 
             (1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
           evdwl *= factor_lj;

src/USER-MISC/fix_smd.cpp

@@ -210,6 +210,10 @@ void FixSMD::smd_tether()
       fy = k_smd*dy*dr/r;
       fz = k_smd*dz*dr/r;
       pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt;
+    } else {
+      fx = 0.0;
+      fy = 0.0;
+      fz = 0.0;
     }
   } else {
     r_old = r;
@@ -306,8 +310,11 @@ void FixSMD::smd_couple()
       fy = k_smd*dy*dr/r;
       fz = k_smd*dz*dr/r;
       pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * update->dt;
+    } else {
+      fx = 0.0;
+      fy = 0.0;
+      fz = 0.0;
     }
-
   } else {
     dx = xcm2[0] - xcm[0];
     dy = xcm2[1] - xcm[1];

src/USER-MISC/pair_srp.cpp

@@ -153,6 +153,7 @@ void PairSRP::compute(int eflag, int vflag)
     double wd, lever0, lever1, evdwl, fpair;
     double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1;
     double fx, fy, fz;
+    evdwl = 0.0;
 
     // mapping global to local for atoms inside bond particles 
     // exclude 1-2 neighs if requested 

src/compute_bond_local.cpp

@@ -172,11 +172,9 @@ int ComputeBondLocal::compute_bonds(int flag)
         domain->minimum_image(delx,dely,delz);
         rsq = delx*delx + dely*dely + delz*delz;
 
-        if (singleflag) {
+        if (singleflag && (btype > 0))
-          if (btype > 0)
+          eng = bond->single(btype,rsq,atom1,atom2,fbond);
-            eng = bond->single(btype,rsq,atom1,atom2,fbond);
+        else eng = fbond = 0.0;
-          else eng = fbond = 0.0;
-        }
 
         if (nvalues == 1) ptr = &vector[m];
         else ptr = array[m];

src/compute_pair_local.cpp

@@ -221,6 +221,7 @@ int ComputePairLocal::compute_pairs(int flag)
       if (flag) {
         if (singleflag)
           eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
+        else eng = fpair = 0.0;
 
         if (nvalues == 1) ptr = &vector[m];
         else ptr = array[m];