Merge pull request #301 from akohlmey/corrections-and-bugfixes

Collected corrections and bugfixes

src/USER-DPD/fix_eos_cv.cpp

@@ -70,8 +70,8 @@ void FixEOScv::init()
       if (mask[i] & groupbit) {
         if(dpdTheta[i] <= 0.0)
           error->one(FLERR,"Internal temperature <= zero");
-        uCond[i] = 0.5*cvEOS*dpdTheta[i];
-        uMech[i] = 0.5*cvEOS*dpdTheta[i];
+        uCond[i] = 0.0;
+        uMech[i] = cvEOS*dpdTheta[i];
       }
   }
 }

src/USER-DPD/fix_eos_cv.h

@@ -50,6 +50,10 @@ 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.
 
+E:  FixEOScv requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory.
+
 E: EOS cv must be > 0.0
 
 The constant volume heat capacity must be larger than zero.

src/USER-DPD/fix_eos_table.cpp

@@ -122,8 +122,8 @@ void FixEOStable::init()
         if(dpdTheta[i] <= 0.0)
           error->one(FLERR,"Internal temperature <= zero");
         energy_lookup(dpdTheta[i],tmp);
-        uCond[i] = tmp / 2.0;
-        uMech[i] = tmp / 2.0;
+        uCond[i] = 0.0;
+        uMech[i] = tmp;
       }
   }
 }

src/USER-DPD/fix_eos_table.h

@@ -92,6 +92,10 @@ E: eos/table values are not increasing
 
 The EOS must be a monotonically, increasing function
 
+E:  FixEOStable requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory.
+
 E: Internal temperature < zero
 
 Self-explanatory.  EOS may not be valid under current simulation conditions.

src/USER-DPD/fix_eos_table_rx.cpp

@@ -162,13 +162,15 @@ void FixEOStableRX::setup(int vflag)
   double *uCG   = atom->uCG;
   double *uCGnew = atom->uCGnew;
 
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit){
-      duChem = uCG[i] - uCGnew[i];
-      uChem[i] += duChem;
-      uCG[i] = 0.0;
-      uCGnew[i] = 0.0;
-    }
+  if(!this->restart_reset){
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit){
+        duChem = uCG[i] - uCGnew[i];
+        uChem[i] += duChem;
+        uCG[i] = 0.0;
+        uCGnew[i] = 0.0;
+      }
+  }
 
   // Communicate the updated momenta and velocities to all nodes
   comm->forward_comm_fix(this);
@@ -200,8 +202,8 @@ void FixEOStableRX::init()
         if(dpdTheta[i] <= 0.0)
           error->one(FLERR,"Internal temperature <= zero");
         energy_lookup(i,dpdTheta[i],tmp);
-        uCond[i] = tmp / 2.0;
-        uMech[i] = tmp / 2.0;
+        uCond[i] = 0.0;
+        uMech[i] = tmp;
         uChem[i] = 0.0;
       }
   }

src/USER-DPD/fix_eos_table_rx.h

@@ -106,6 +106,10 @@ E:  eos/table/rx values are not increasing
 
 The equation-of-state must an increasing function
 
+E:  FixEOStableRX requires atom_style with internal temperature and energies (e.g. dpd)
+
+Self-explanatory.
+
 E:  Internal temperature <= zero.
 
 Self-explanatory.

src/USER-DPD/fix_rx.cpp

@@ -69,7 +69,6 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
   id_fix_species_old(NULL), fix_species(NULL), fix_species_old(NULL)
 {
   if (narg < 7 || narg > 12) error->all(FLERR,"Illegal fix rx command");
-  restart_peratom = 1;
   nevery = 1;
 
   nreactions = maxparam = 0;
@@ -366,6 +365,7 @@ void FixRX::post_constructor()
 
   modify->add_fix(nspecies+5,newarg);
   fix_species = (FixPropertyAtom *) modify->fix[modify->nfix-1];
+  restartFlag = modify->fix[modify->nfix-1]->restart_reset;
 
   modify->add_fix(nspecies+5,newarg2);
   fix_species_old = (FixPropertyAtom *) modify->fix[modify->nfix-1];
@@ -647,34 +647,38 @@ void FixRX::setup_pre_force(int vflag)
   double tmp;
   int ii;
 
-  if(localTempFlag){
-    int count = nlocal + (newton_pair ? nghost : 0);
-    dpdThetaLocal = new double[count];
-    memset(dpdThetaLocal, 0, sizeof(double)*count);
-    computeLocalTemperature();
+  if(restartFlag){
+    restartFlag = 0;
+  } else {
+    if(localTempFlag){
+      int count = nlocal + (newton_pair ? nghost : 0);
+      dpdThetaLocal = new double[count];
+      memset(dpdThetaLocal, 0, sizeof(double)*count);
+      computeLocalTemperature();
+    }
+  
+    for (int id = 0; id < nlocal; id++)
+      for (int ispecies=0; ispecies<nspecies; ispecies++){
+        tmp = atom->dvector[ispecies][id];
+        atom->dvector[ispecies+nspecies][id] = tmp;
+      }
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit){
+
+        // Set the reaction rate constants to zero:  no reactions occur at step 0
+        for(int irxn=0;irxn<nreactions;irxn++)
+	  kR[irxn] = 0.0;
+
+        if (odeIntegrationFlag == ODE_LAMMPS_RK4)
+	  rk4(i,NULL);
+        else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
+	  rkf45(i,NULL);
+      }
+  
+    // Communicate the updated momenta and velocities to all nodes
+    comm->forward_comm_fix(this);
+    if(localTempFlag) delete [] dpdThetaLocal;
   }
-
-  for (int id = 0; id < nlocal; id++)
-    for (int ispecies=0; ispecies<nspecies; ispecies++){
-      tmp = atom->dvector[ispecies][id];
-      atom->dvector[ispecies+nspecies][id] = tmp;
-    }
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit){
-
-      // Set the reaction rate constants to zero:  no reactions occur at step 0
-      for(int irxn=0;irxn<nreactions;irxn++)
-        kR[irxn] = 0.0;
-
-      if (odeIntegrationFlag == ODE_LAMMPS_RK4)
-        rk4(i,NULL);
-      else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
-        rkf45(i,NULL);
-    }
-
-  // Communicate the updated momenta and velocities to all nodes
-  comm->forward_comm_fix(this);
-  if(localTempFlag) delete [] dpdThetaLocal;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-DPD/fix_rx.h

@@ -132,6 +132,7 @@ class FixRX : public Fix {
   char *kineticsFile;
   char *id_fix_species,*id_fix_species_old;
   class FixPropertyAtom *fix_species,*fix_species_old;
+  int restartFlag;
 
 };
 

src/USER-DPD/pair_exp6_rx.cpp

@@ -310,13 +310,13 @@ void PairExp6rx::compute(int eflag, int vflag)
 
             uin1 = buck1*(6.0*rin1exp - alphaOld12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
 
-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;
 
-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;
 
             uin1rep = aRep/powint(rin1,nRep);
 
-            forceExp6 = -double(nRep)*aRep/powint(r,nRep);
+            forceExp6 = double(nRep)*aRep/powint(r,nRep);
             fpairOldEXP6_12 = factor_lj*forceExp6*r2inv;
 
             evdwlOldEXP6_12 = uin1 - uin1rep + aRep/powint(r,nRep);
@@ -350,13 +350,13 @@ void PairExp6rx::compute(int eflag, int vflag)
 
             uin1 = buck1*(6.0*rin1exp - alphaOld21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
 
-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;
 
-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;
 
             uin1rep = aRep/powint(rin1,nRep);
 
-            forceExp6 = -double(nRep)*aRep/powint(r,nRep);
+            forceExp6 = double(nRep)*aRep/powint(r,nRep);
             fpairOldEXP6_21 = factor_lj*forceExp6*r2inv;
 
             evdwlOldEXP6_21 = uin1 - uin1rep + aRep/powint(r,nRep);
@@ -405,9 +405,9 @@ void PairExp6rx::compute(int eflag, int vflag)
 
             uin1 = buck1*(6.0*rin1exp - alpha12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
 
-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;
 
-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;
 
             uin1rep = aRep/powint(rin1,nRep);
 
@@ -437,9 +437,9 @@ void PairExp6rx::compute(int eflag, int vflag)
 
             uin1 = buck1*(6.0*rin1exp - alpha21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
 
-            win1 = -buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
+            win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) + rin1*durc;
 
-            aRep = -1.0*win1*powint(rin1,nRep)/nRep;
+            aRep = win1*powint(rin1,nRep)/nRep;
 
             uin1rep = aRep/powint(rin1,nRep);
 
@@ -1102,6 +1102,32 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
       rm2_old *= pow(nTotalOFA_old,fuchslinR);
     }
   }
+
+  // Check that no fractions are less than zero
+  if(fraction1 < 0.0){
+    if(fraction1 < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction1 = 0.0;
+  }
+  if(fraction2 < 0.0){
+    if(fraction2 < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction2 = 0.0;
+  }
+  if(fraction1_old < 0.0){
+    if(fraction1_old < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction1_old = 0.0;
+  }
+  if(fraction2_old < 0.0){
+    if(fraction2_old < -1.0e-10){
+      error->all(FLERR,"Computed fraction less than -1.0e-10");
+    }
+    fraction2_old = 0.0;
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-DPD/pair_multi_lucy.h

@@ -18,7 +18,7 @@ PairStyle(multi/lucy,PairMultiLucy)
 #else
 
 #ifndef LMP_PAIR_MULTI_LUCY_H
-#define LMP_PAIR_MUTLI_LUCY_H
+#define LMP_PAIR_MULTI_LUCY_H
 
 #include "pair.h"
 

src/USER-DPD/pair_multi_lucy_rx.h

@@ -18,7 +18,7 @@ PairStyle(multi/lucy/rx,PairMultiLucyRX)
 #else
 
 #ifndef LMP_PAIR_MULTI_LUCY_RX_H
-#define LMP_PAIR_MUTLI_LUCY_RX_H
+#define LMP_PAIR_MULTI_LUCY_RX_H
 
 #include "pair.h"
 

src/memory.cpp

@@ -26,6 +26,10 @@
 #endif
 #endif
 
+#if defined(LMP_USER_INTEL) && !defined(LAMMPS_MEMALIGN)
+#define LAMMPS_MEMALIGN 64
+#endif
+
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */