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

src/USER-EFF/README

@@ -74,3 +74,18 @@ real units, upgraded post-processing and visualization tools, added
 support for "compute pair eff" to allow thermo prints with the
 different eFF energy components (eke, epauli, ecoul and errestrain),
 fixed radial scaling factors in the eff langevin thermostat.
+
+12/2011: Added support for "zero" option in fix langevin/eff (see doc), and
+adjusted fix_langevin_eff.cpp to correctly thermostat between nuclear and electronic dof 
+(required additional scaling of friction term in the Langevin equations of motion).
+Radial electron mass now scales as a function of system dimension.
+
+Bug fixes:
+(10-2011): Thanks to Christian Chenard-Lemire (U Montreal) for reporting a bug in the
+fixed pair_eff_cut.cpp fixedcore-pseudocore interactions (an incorrect index and a missing
+elec-core call to account for the 2 electrons from the fixed core.)     
+(12-2011): Corrected undefined natom variable in fix_langevin_eff (recent changes in 
+main fix_langevin class caused compilation error in user-eff).  
+Corrected thermostat in fix langevin/eff as described in version 12/2011.
+
+

src/USER-EFF/compute_ke_atom_eff.cpp

@@ -25,6 +25,7 @@
 #include "force.h"
 #include "memory.h"
 #include "error.h"
+#include "domain.h"
 
 using namespace LAMMPS_NS;
 
@@ -90,6 +91,7 @@ void ComputeKEAtomEff::compute_peratom()
   int *mask = atom->mask;
   int *type = atom->type;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
   
   if (mass)
     for (int i = 0; i < nlocal; i++) {
@@ -97,7 +99,7 @@ void ComputeKEAtomEff::compute_peratom()
         ke[i] = 0.5 * mvv2e * mass[type[i]] * 
           (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
         if (fabs(spin[i])==1)
-          ke[i] += 0.5 * mvv2e * mass[type[i]] * ervel[i]*ervel[i] * 0.75;
+          ke[i] += 0.5 * mvv2e * mefactor * mass[type[i]] * ervel[i]*ervel[i];
       } else ke[i] = 0.0;
     }
 }

src/USER-EFF/compute_ke_eff.cpp

@@ -63,6 +63,7 @@ double ComputeKEEff::compute_scalar()
   int *mask = atom->mask;
   int *type = atom->type;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
 
   double ke = 0.0;
 
@@ -71,7 +72,7 @@ double ComputeKEEff::compute_scalar()
       if (mask[i] & groupbit) {
         ke += mass[type[i]]*(v[i][0]*v[i][0] + v[i][1]*v[i][1] + 
                              v[i][2]*v[i][2]);
-        if (fabs(spin[i])==1) ke += 0.75*mass[type[i]]*ervel[i]*ervel[i];
+        if (fabs(spin[i])==1) ke += mefactor*mass[type[i]]*ervel[i]*ervel[i];
       }
   }
 

src/USER-EFF/compute_temp_deform_eff.cpp

@@ -113,6 +113,7 @@ void ComputeTempDeformEff::dof_compute()
   int nelectrons;
   MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
 
+  // Assume 3/2 k T per nucleus 
   dof -= domain->dimension * nelectrons;
   
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
@@ -135,6 +136,7 @@ double ComputeTempDeformEff::compute_scalar()
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
 
   // lamda = 0-1 triclinic lamda coords
   // vstream = streaming velocity = Hrate*lamda + Hratelo
@@ -159,7 +161,7 @@ double ComputeTempDeformEff::compute_scalar()
       if (mass) {
         t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + 
               vthermal[2]*vthermal[2])* mass[type[i]];
-        if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
+        if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i];
       }
     }
   
@@ -185,6 +187,7 @@ void ComputeTempDeformEff::compute_vector()
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
 
   double *h_rate = domain->h_rate;
   double *h_ratelo = domain->h_ratelo;
@@ -211,9 +214,9 @@ void ComputeTempDeformEff::compute_vector()
       t[4] += massone * vthermal[0]*vthermal[2];
       t[5] += massone * vthermal[1]*vthermal[2];
       if (fabs(spin[i])==1) {
-        t[0] += 0.75 * massone * ervel[i]*ervel[i];
+        t[0] += mefactor * massone * ervel[i]*ervel[i];
-        t[1] += 0.75 * massone * ervel[i]*ervel[i];
+        t[1] += mefactor * massone * ervel[i]*ervel[i];
-        t[2] += 0.75 * massone * ervel[i]*ervel[i];
+        t[2] += mefactor * massone * ervel[i]*ervel[i];
       }
     }
   

src/USER-EFF/compute_temp_eff.cpp

@@ -84,7 +84,7 @@ void ComputeTempEff::dof_compute()
   int nelectrons;
   MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
 
-  // average over nuclear dof only
+  // Assume 3/2 k T per nucleus 
 
   dof -= domain->dimension * nelectrons;
 
@@ -105,6 +105,7 @@ double ComputeTempEff::compute_scalar()
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
 
   double t = 0.0;
    
@@ -113,7 +114,7 @@ double ComputeTempEff::compute_scalar()
       if (mask[i] & groupbit) {
         t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * 
           mass[type[i]];
-        if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
+        if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i];
       }
     }
   }
@@ -139,6 +140,7 @@ void ComputeTempEff::compute_vector()
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
 
   double massone,t[6];
   for (i = 0; i < 6; i++) t[i] = 0.0;
@@ -153,9 +155,9 @@ void ComputeTempEff::compute_vector()
       t[4] += massone * v[i][0]*v[i][2];
       t[5] += massone * v[i][1]*v[i][2];
       if (fabs(spin[i])==1) {
-        t[0] += 0.75*massone*ervel[i]*ervel[i];
+        t[0] += mefactor*massone*ervel[i]*ervel[i];
-        t[1] += 0.75*massone*ervel[i]*ervel[i];
+        t[1] += mefactor*massone*ervel[i]*ervel[i];
-        t[2] += 0.75*massone*ervel[i]*ervel[i];
+        t[2] += mefactor*massone*ervel[i]*ervel[i];
       }
     }
 

src/USER-EFF/compute_temp_region_eff.cpp

@@ -104,9 +104,11 @@ double ComputeTempRegionEff::compute_scalar()
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
 
   Region *region = domain->regions[iregion];
   int count = 0;
+  int ecount = 0;
   double t = 0.0;
 
   if (mass) {
@@ -115,12 +117,16 @@ double ComputeTempRegionEff::compute_scalar()
         count++;
         t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * 
           mass[type[i]];
-        if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
+        if (fabs(spin[i])==1) {
+          t += mefactor*mass[type[i]]*ervel[i]*ervel[i];
+          ecount++;
+        }
       }
   }
 
   double tarray[2],tarray_all[2];
-  tarray[0] = count;
+  // Assume 3/2 k T per nucleus
+  tarray[0] = count-ecount;
   tarray[1] = t;
   MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world);
   dof = domain->dimension * tarray_all[0] - extra_dof;
@@ -130,12 +136,6 @@ double ComputeTempRegionEff::compute_scalar()
     if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
       if (fabs(spin[i])==1) one++;
     }
-  int nelectrons_region;
-  MPI_Allreduce(&one,&nelectrons_region,1,MPI_INT,MPI_SUM,world);
-
-  // average over nuclear dof only
-
-  dof -= domain->dimension * nelectrons_region ;
 
   if (dof > 0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz);
   else scalar = 0.0;
@@ -158,6 +158,7 @@ void ComputeTempRegionEff::compute_vector()
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
 
   Region *region = domain->regions[iregion];
   double massone,t[6];
@@ -175,9 +176,9 @@ void ComputeTempRegionEff::compute_vector()
       t[5] += massone * v[i][1]*v[i][2];
                 
       if (fabs(spin[i])==1) {
-        t[0] += 0.75 * massone * ervel[i]*ervel[i];
+        t[0] += mefactor * massone * ervel[i]*ervel[i];
-        t[1] += 0.75 * massone * ervel[i]*ervel[i];
+        t[1] += mefactor * massone * ervel[i]*ervel[i];
-        t[2] += 0.75 * massone * ervel[i]*ervel[i];
+        t[2] += mefactor * massone * ervel[i]*ervel[i];
       }
     }
 

src/USER-EFF/fix_langevin_eff.cpp

@@ -20,6 +20,7 @@
 #include "string.h"
 #include "stdlib.h"
 #include "fix_langevin_eff.h"
+#include "math_extra.h"
 #include "atom.h"
 #include "force.h"
 #include "update.h"
@@ -29,13 +30,20 @@
 #include "region.h"
 #include "respa.h"
 #include "comm.h"
+#include "input.h"
+#include "variable.h"
 #include "random_mars.h"
 #include "memory.h"
 #include "error.h"
+#include "group.h"
 
 using namespace LAMMPS_NS;
 
 enum{NOBIAS,BIAS};
+enum{CONSTANT,EQUAL,ATOM};
+
+#define SINERTIA 0.4          // moment of inertia prefactor for sphere
+#define EINERTIA 0.2          // moment of inertia prefactor for ellipsoid
 
 /* ---------------------------------------------------------------------- */
 
@@ -56,22 +64,51 @@ FixLangevinEff::~FixLangevinEff()
 
 void FixLangevinEff::post_force_no_tally()
 {
-  double gamma1,gamma2;
+  double gamma1,gamma2,t_target;
 
   double **v = atom->v;
   double **f = atom->f;
-  int *type = atom->type;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
   double *ervel = atom->ervel;
   double *erforce = atom->erforce;
   int *spin = atom->spin;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  double mefactor = domain->dimension/4.0;
+  double sqrtmefactor = sqrt(mefactor);
 
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
-  double t_target = t_start + delta * (t_stop-t_start);
+
-  double tsqrt = sqrt(t_target);
+  // set current t_target and t_sqrt
+  // if variable temp, evaluate variable, wrap with clear/add
+  // reallocate tforce array if necessary
+
+  if (tstyle == CONSTANT) {
+    t_target = t_start + delta * (t_stop-t_start);
+    tsqrt = sqrt(t_target);
+  } else {
+    modify->clearstep_compute();
+    if (tstyle == EQUAL) {
+      t_target = input->variable->compute_equal(tvar);
+      if (t_target < 0.0)
+	error->one(FLERR,"Fix langevin/eff variable returned negative temperature");
+      tsqrt = sqrt(t_target);
+    } else {
+      if (nlocal > maxatom2) {
+	maxatom2 = atom->nmax;
+	memory->destroy(tforce);
+	memory->create(tforce,maxatom2,"langevin/eff:tforce");
+      }
+      input->variable->compute_atom(tvar,igroup,tforce,1,0);
+      for (int i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit)
+	    if (tforce[i] < 0.0) 
+	      error->one(FLERR,
+			 "Fix langevin/eff variable returned negative temperature");
+    }
+    modify->addstep_compute(update->ntimestep + 1);
+  }
 
   // apply damping and thermostat to atoms in group
   // for BIAS:
@@ -79,72 +116,179 @@ void FixLangevinEff::post_force_no_tally()
   //   computed on current nlocal atoms to remove bias
   //   test v = 0 since some computes mask non-participating atoms via v = 0
   //   and added force has extra term not multiplied by v = 0
+  // for ZEROFLAG:
+  //   sum random force over all atoms in group
+  //   subtract sum/particles from each atom in group    
+
+  double fran[4],fsum[4],fsumall[4];
+  fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
+
+  double boltz = force->boltz;
+  double dt = update->dt;
+  double mvv2e = force->mvv2e;
+  double ftm2v = force->ftm2v;
+
+  int particles = group->count(igroup);
+  if (zeroflag) {
+    if (particles == 0)
+      error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons");
+  }
+
+  // find number of electrons in group
+  int dof,fix_dof;
+  dof = domain->dimension * particles;
+  fix_dof = 0;
+  for (int i = 0; i < modify->nfix; i++)
+    fix_dof += modify->fix[i]->dof(igroup);
+
+  // extra_dof = domain->dimension
+  dof -= domain->dimension + fix_dof;
+
+  int one = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (fabs(spin[i])==1) one++;
+    }
+  int nelectrons, dofelectrons, dofnuclei;
+  MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
+  dofelectrons = domain->dimension*nelectrons;
+  dofnuclei = dof-dofelectrons;
+
+  // thermal partitioning factor between nuclei and electrons 
+  // extra dof from electron size
+  double gfactor3=(double) (dof+nelectrons)/dofnuclei;
 
   if (which == NOBIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
-        gamma1 = gfactor1[type[i]];
+        if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
+	gamma1 = gfactor1[type[i]] * gfactor3;
 	gamma2 = gfactor2[type[i]] * tsqrt;
-	f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
+	fran[0] = gamma2*(random->uniform()-0.5);
-	f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
+	fran[1] = gamma2*(random->uniform()-0.5);
-	f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+	fran[2] = gamma2*(random->uniform()-0.5);
-        if (abs(spin[i])==1) 
+	f[i][0] += gamma1*v[i][0] + fran[0];
-	  erforce[i] += 0.75*gamma1*ervel[i] + 
+	f[i][1] += gamma1*v[i][1] + fran[1];
-	    0.866025404*gamma2*(random->uniform()-0.5);
+	f[i][2] += gamma1*v[i][2] + fran[2];
+	fsum[0] += fran[0];
+	fsum[1] += fran[1];
+	fsum[2] += fran[2];
+        if (fabs(spin[i])==1) {
+          fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5);
+          erforce[i] += mefactor*gamma1*ervel[i]+fran[3];
+          fsum[3] += fran[3];
+        }
       }
     }
   } else if (which == BIAS) {
-    double tmp = temperature->compute_scalar();
+    temperature->compute_scalar();
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
-	gamma1 = gfactor1[type[i]];
+	if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
+	gamma1 = gfactor1[type[i]] * gfactor3;
 	gamma2 = gfactor2[type[i]] * tsqrt;
 	temperature->remove_bias(i,v[i]);
+	fran[0] = gamma2*(random->uniform()-0.5);
+	fran[1] = gamma2*(random->uniform()-0.5);
+	fran[2] = gamma2*(random->uniform()-0.5);
 	if (v[i][0] != 0.0)
-	  f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
+	  f[i][0] += gamma1*v[i][0] + fran[0];
 	if (v[i][1] != 0.0)
-	  f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
+	  f[i][1] += gamma1*v[i][1] + fran[1];
 	if (v[i][2] != 0.0)
-	  f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+	  f[i][2] += gamma1*v[i][2] + fran[2];
-        if (abs(spin[i])==1 && ervel[i] != 0.0)
+	fsum[0] += fran[0];
-          erforce[i] += 0.75*gamma1*ervel[i] + 
+	fsum[1] += fran[1];
-	    0.866025404*gamma2*(random->uniform()-0.5);
+	fsum[2] += fran[2];
+        if (fabs(spin[i])==1) {
+          fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5);
+          if (ervel[i] != 0.0) erforce[i] += mefactor*gamma1*ervel[i]+fran[3];
+          fsum[3] += fran[3];
+        }
 	temperature->restore_bias(i,v[i]);
       }
     }
   }
+
+  // set total force to zero
+
+  if (zeroflag) {
+    MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
+    fsumall[0] /= particles;
+    fsumall[1] /= particles;
+    fsumall[2] /= particles;
+    fsumall[3] /= nelectrons;
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	f[i][0] -= fsumall[0];
+	f[i][1] -= fsumall[1];
+	f[i][2] -= fsumall[2];
+        if (fabs(spin[i])==1) erforce[i] -= fsumall[3];
+      }
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixLangevinEff::post_force_tally()
 {
-  double gamma1,gamma2;
+  double gamma1,gamma2,t_target;
 
-  // reallocate flangevin if necessary
+  // reallocate flangevin and erforcelangevin if necessary
 
-  if (atom->nmax > nmax) {
+  if (atom->nlocal > maxatom1) {
     memory->destroy(flangevin);
     memory->destroy(erforcelangevin);
-    nmax = atom->nmax;
+    maxatom1 = atom->nmax;
-    memory->create(flangevin,nmax,3,"langevin:flangevin");
+    memory->create(flangevin,maxatom1,3,"langevin/eff:flangevin");
-    memory->create(erforcelangevin,nmax,"langevin/eff:erforcelangevin");
+    memory->create(erforcelangevin,maxatom1,"langevin/eff:erforcelangevin");
   }
 
   double **v = atom->v;
   double **f = atom->f;
+  double *erforce = atom->erforce;
+  double *ervel = atom->ervel;
+  int *spin = atom->spin;
+  double mefactor = domain->dimension/4.0;
+  double sqrtmefactor = sqrt(mefactor);
+
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double *erforce = atom->erforce;
-  double *ervel = atom->ervel;
-  int *spin = atom->spin;
-
   double delta = update->ntimestep - update->beginstep;
   delta /= update->endstep - update->beginstep;
-  double t_target = t_start + delta * (t_stop-t_start);
+
-  double tsqrt = sqrt(t_target);
+  // set current t_target and t_sqrt
+  // if variable temp, evaluate variable, wrap with clear/add
+  // reallocate tforce array if necessary
+
+  if (tstyle == CONSTANT) {
+    t_target = t_start + delta * (t_stop-t_start);
+    tsqrt = sqrt(t_target);
+  } else {
+    modify->clearstep_compute();
+    if (tstyle == EQUAL) {
+      t_target = input->variable->compute_equal(tvar);
+      if (t_target < 0.0)
+	error->one(FLERR,"Fix langevin/eff variable returned negative temperature");
+      tsqrt = sqrt(t_target);
+    } else {
+      if (nlocal > maxatom2) {
+	maxatom2 = atom->nmax;
+	memory->destroy(tforce);
+	memory->create(tforce,maxatom2,"langevin/eff:tforce");
+      }
+      input->variable->compute_atom(tvar,igroup,tforce,1,0);
+      for (int i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit)
+	    if (tforce[i] < 0.0) 
+	      error->one(FLERR,
+			 "Fix langevin/eff variable returned negative temperature");
+    }
+    modify->addstep_compute(update->ntimestep + 1);
+  }
 
   // apply damping and thermostat to appropriate atoms
   // for BIAS:
@@ -153,42 +297,81 @@ void FixLangevinEff::post_force_tally()
   //   test v = 0 since some computes mask non-participating atoms via v = 0
   //   and added force has extra term not multiplied by v = 0
 
+  double boltz = force->boltz;
+  double dt = update->dt;
+  double mvv2e = force->mvv2e;
+  double ftm2v = force->ftm2v;
+
+  int particles = group->count(igroup);
+  if (zeroflag) {
+    if (particles == 0)
+      error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons");
+  }
+
+  // find number of electrons in group
+  int dof,fix_dof;
+  dof = domain->dimension * particles;
+  fix_dof = 0;
+  for (int i = 0; i < modify->nfix; i++)
+    fix_dof += modify->fix[i]->dof(igroup);
+
+  // extra_dof = domain->dimension
+  dof -= domain->dimension + fix_dof;
+
+  int one = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (fabs(spin[i])==1) one++;
+    }
+  int nelectrons, dofelectrons, dofnuclei;
+  MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
+  dofelectrons = domain->dimension*nelectrons;
+  dofnuclei = dof-dofelectrons;
+
+  // thermal partitioning factor between nuclei and electrons
+  // extra dof from electron size
+  double gfactor3=(double) (dof+nelectrons)/dofnuclei;
+
   if (which == NOBIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
-	gamma1 = gfactor1[type[i]];
+	if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
+	gamma1 = gfactor1[type[i]] * gfactor3;
 	gamma2 = gfactor2[type[i]] * tsqrt;
 	flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
 	flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
 	flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-        erforcelangevin[i] = 0.75*gamma1*ervel[i] + 
-	  0.866025404*gamma2*(random->uniform()-0.5);
 	f[i][0] += flangevin[i][0];
 	f[i][1] += flangevin[i][1];
 	f[i][2] += flangevin[i][2];
-        if (abs(spin[i])==1) erforce[i] += erforcelangevin[i];
+        if (fabs(spin[i])==1) {
+          erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5);
+          erforce[i] += erforcelangevin[i];
+        }
       }
     }
   } else if (which == BIAS) {
-    double tmp = temperature->compute_scalar();
+    temperature->compute_scalar();
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
-	gamma1 = gfactor1[type[i]];
+	if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
+	gamma1 = gfactor1[type[i]] * gfactor3;
 	gamma2 = gfactor2[type[i]] * tsqrt;
 	temperature->remove_bias(i,v[i]);
 	flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
 	flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
 	flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
-        erforcelangevin[i] = 0.75*gamma1*ervel[i] + 
-	  0.866025404*gamma2*(random->uniform()-0.5);
 	if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
 	else flangevin[i][0] = 0.0;
 	if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
 	else flangevin[i][1] = 0.0;
 	if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
 	else flangevin[i][2] = 0.0;
-        if (abs(spin[i])==1 && ervel[i] != 0.0)
+        if (fabs(spin[i])==1) {
-	  erforce[i] += erforcelangevin[i];
+          erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5);
+          if (ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
+          else erforcelangevin[i] = 0.0;
+        }
 	temperature->restore_bias(i,v[i]);
       }
     }
@@ -214,7 +397,7 @@ void FixLangevinEff::end_of_step()
     if (mask[i] & groupbit) {
       energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + 
           flangevin[i][2]*v[i][2];
-      if (abs(spin[i])==1) energy_onestep += erforcelangevin[i];
+      if (fabs(spin[i])==1) energy_onestep += erforcelangevin[i];
     }
   energy += energy_onestep*update->dt;
 }
@@ -223,7 +406,7 @@ void FixLangevinEff::end_of_step()
 
 double FixLangevinEff::compute_scalar()
 {
-  if (!tally) return 0.0;
+  if (!tally || flangevin == NULL || erforcelangevin == NULL) return 0.0;
 
   // capture the very first energy transfer to thermal reservoir
 
@@ -238,7 +421,7 @@ double FixLangevinEff::compute_scalar()
       if (mask[i] & groupbit) {
 	energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + 
 	  flangevin[i][2]*v[i][2];
-        if (abs(spin[i])==1) energy_onestep += erforcelangevin[i];
+        if (fabs(spin[i])==1) energy_onestep += erforcelangevin[i];
       }
     energy = 0.5*energy_onestep*update->dt;
   }
@@ -249,3 +432,4 @@ double FixLangevinEff::compute_scalar()
   MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);	
   return -energy_all;
 }
+

src/USER-EFF/fix_nh_eff.cpp

@@ -21,6 +21,7 @@
 #include "atom_vec.h"
 #include "group.h"
 #include "error.h"
+#include "domain.h"
 
 using namespace LAMMPS_NS;
 
@@ -48,6 +49,7 @@ void FixNHEff::nve_v()
   double *ervel = atom->ervel;
   double *mass = atom->mass;
   int *spin = atom->spin;
+  double mefactor = domain->dimension/4.0;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -60,7 +62,7 @@ void FixNHEff::nve_v()
     if (mask[i] & groupbit) {
       if (fabs(spin[i])==1) {
 	dtfm = dtf / mass[type[i]];
-	ervel[i] = dtfm * erforce[i] / 0.75;
+	ervel[i] = dtfm * erforce[i] / mefactor;
       }
     }
   }

src/USER-EFF/fix_nve_eff.cpp

@@ -25,6 +25,7 @@
 #include "respa.h"
 #include "error.h"
 #include "math.h"
+#include "domain.h"
 
 using namespace LAMMPS_NS;
 
@@ -80,6 +81,7 @@ void FixNVEEff::initial_integrate(int vflag)
   double *erforce = atom->erforce;
   double *mass = atom->mass;
   int *spin = atom->spin;
+  double mefactor = domain->dimension/4.0;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -98,7 +100,7 @@ void FixNVEEff::initial_integrate(int vflag)
         x[i][1] += dtv * v[i][1];
         x[i][2] += dtv * v[i][2];
         if (fabs(spin[i])==1) {
-          ervel[i] += dtfm * erforce[i] / 0.75;
+          ervel[i] += dtfm * erforce[i] / mefactor;
           eradius[i] += dtv * ervel[i];
         }
       }
@@ -118,6 +120,7 @@ void FixNVEEff::final_integrate()
   double **f = atom->f;
   double *mass = atom->mass;
   int *spin = atom->spin;
+  double mefactor = domain->dimension/4.0;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -133,7 +136,7 @@ void FixNVEEff::final_integrate()
         v[i][1] += dtfm * f[i][1];
         v[i][2] += dtfm * f[i][2];
         if (fabs(spin[i])==1)
-          ervel[i] += dtfm * erforce[i] / 0.75;
+          ervel[i] += dtfm * erforce[i] / mefactor;
       }
     }
   }

src/USER-EFF/pair_eff_cut.cpp

@@ -402,6 +402,11 @@ void PairEffCut::compute(int eflag, int vflag)
             ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair);
             ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
                         &fpair,&e2rforce);
+            ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
+                        &fpair,&e2rforce);
+            PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
+                        &ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B,
+                        PAULI_CORE_C);
             PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
                         &ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B,
                         PAULI_CORE_C);
@@ -427,7 +432,7 @@ void PairEffCut::compute(int eflag, int vflag)
         else if ((abs(spin[i]) == 1 || spin[i] == 2) && spin[j] == 3) {
           e1rforce = ecp_e1rforce = 0.0;
 
-          if (abs(spin[j]) == 1) {
+          if (abs(spin[i]) == 1) {
             ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
                         &fpair,&e1rforce);
             PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
@@ -438,6 +443,11 @@ void PairEffCut::compute(int eflag, int vflag)
             ElecCoreNuc(qxq,rc,eradius[i],&ecoul,&fpair);
             ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
                         &fpair,&e1rforce);
+            ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
+                        &fpair,&e1rforce);
+            PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
+                        &ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B,
+                        PAULI_CORE_C);
             PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
                         &ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B,
                         PAULI_CORE_C);