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

src/MANYBODY/fix_qeq_comb.cpp

@@ -110,7 +110,7 @@ int FixQEQComb::setmask()
 void FixQEQComb::init()
 {
   if (!atom->q_flag)
-    error->all("Atoms must have charge to use fix qeq/comb");
+    error->all("Fix qeq/comb requires atom attribute q");
 
   comb = (PairComb *) force->pair_match("comb",1);
   if (comb == NULL) error->all("Must use pair_style comb with fix qeq/comb");

src/MANYBODY/pair_comb.cpp

@@ -14,7 +14,7 @@
 /* ----------------------------------------------------------------------
    Contributing author: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
    LAMMPS implementation of the Charge-optimized many-body (COMB) potential 
-   based on HELL MD program (Prof Simon R. Phillpot, UF, sphil@mse.ufl.edu)
+   based on the HELL MD program (Prof Simon Phillpot, UF, sphil@mse.ufl.edu)
    and Aidan Thompson's Tersoff code in LAMMPS
 ------------------------------------------------------------------------- */
 

src/Makefile

@@ -16,8 +16,8 @@ OBJ = 	$(SRC:.cpp=.o)
 PACKAGE = asphere class2 colloid dipole dsmc gpu granular \
 	  kspace manybody meam molecule opt peri poems prd reax shock xtc
 
-PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \
-	   user-imd user-smd
+PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm \
+	   user-ewaldn user-imd user-smd
 
 PACKALL = $(PACKAGE) $(PACKUSER)
 

src/atom.cpp

@@ -73,6 +73,9 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   vfrac = s0 = NULL;
   x0 = NULL;
 
+  spin = NULL;
+  eradius = evel = eforce = NULL;
+
   maxspecial = 1;
   nspecial = NULL;
   special = NULL;
@@ -98,6 +101,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   q_flag = mu_flag = 0;
   quat_flag = omega_flag = angmom_flag = torque_flag = 0;
   radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
+  spin_flag = eradius_flag = evel_flag = eforce_flag = 0;
 
   // ntype-length arrays
 
@@ -176,6 +180,11 @@ Atom::~Atom()
   memory->sfree(s0);
   memory->destroy_2d_double_array(x0);
 
+  memory->sfree(spin);
+  memory->sfree(eradius);
+  memory->sfree(evel);
+  memory->sfree(eforce);
+
   memory->sfree(molecule);
 
   memory->destroy_2d_int_array(nspecial);

src/atom.h

@@ -52,6 +52,9 @@ class Atom : protected Pointers {
   double *radius,*density,*rmass,*vfrac,*s0;
   double **x0;
 
+  int *spin;
+  double *eradius,*evel,*eforce;
+
   int **nspecial;               // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs
   int **special;                // IDs of 1-2,1-3,1-4 neighs of each atom
   int maxspecial;               // special[nlocal][maxspecial]
@@ -80,6 +83,7 @@ class Atom : protected Pointers {
   int q_flag,mu_flag;
   int quat_flag,omega_flag,angmom_flag,torque_flag;
   int radius_flag,density_flag,rmass_flag,vfrac_flag;
+  int spin_flag,eradius_flag,evel_flag,eforce_flag;
 
   // extra peratom info in restart file destined for fix & diag 
 

src/compute_temp.h

@@ -27,16 +27,16 @@ namespace LAMMPS_NS {
 class ComputeTemp : public Compute {
  public:
   ComputeTemp(class LAMMPS *, int, char **);
-  ~ComputeTemp();
+  virtual ~ComputeTemp();
   void init();
   double compute_scalar();
   void compute_vector();
 
- private:
+ protected:
   int fix_dof;
   double tfactor;
 
-  void dof_compute();
+  virtual void dof_compute();
 };
 
 }

src/compute_temp_deform.h

@@ -27,10 +27,10 @@ namespace LAMMPS_NS {
 class ComputeTempDeform : public Compute {
  public:
   ComputeTempDeform(class LAMMPS *, int, char **);
-  ~ComputeTempDeform();
+  virtual ~ComputeTempDeform();
   void init();
-  double compute_scalar();
-  void compute_vector();
+  virtual double compute_scalar();
+  virtual void compute_vector();
 
   void remove_bias(int, double *);
   void remove_bias_all();
@@ -38,14 +38,14 @@ class ComputeTempDeform : public Compute {
   void restore_bias_all();
   double memory_usage();
 
- private:
+ protected:
   int fix_dof;
   double tfactor;
   double vbias[3];    // stored velocity bias for one atom
   double **vbiasall;  // stored velocity bias for all atoms
   int maxbias;        // size of vbiasall array
 
-  void dof_compute();
+  virtual void dof_compute();
 };
 
 }

src/compute_temp_region.h

@@ -27,10 +27,10 @@ namespace LAMMPS_NS {
 class ComputeTempRegion : public Compute {
  public:
   ComputeTempRegion(class LAMMPS *, int, char **);
-  ~ComputeTempRegion();
+  virtual ~ComputeTempRegion();
   void init();
-  double compute_scalar();
-  void compute_vector();
+  virtual double compute_scalar();
+  virtual void compute_vector();
 
   int dof_remove(int);
   void remove_bias(int, double *);
@@ -39,7 +39,7 @@ class ComputeTempRegion : public Compute {
   void restore_bias_all();
   double memory_usage();
 
- private:
+ protected:
   int iregion;
   char *idregion;
 };

src/dump_custom.cpp

@@ -11,6 +11,7 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "stdlib.h"
 #include "string.h"
 #include "dump_custom.h"
@@ -37,7 +38,7 @@ enum{ID,MOL,TYPE,MASS,
        X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ,
        VX,VY,VZ,FX,FY,FZ,
        Q,MUX,MUY,MUZ,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
-       QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,
+       QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,EVEL,EFORCE,
        COMPUTE,FIX,VARIABLE};
 enum{LT,LE,GT,GE,EQ,NEQ};
 enum{INT,DOUBLE};
@@ -586,6 +587,7 @@ int DumpCustom::count()
       } else if (thresh_array[ithresh] == FZ) {
 	ptr = &atom->f[0][2];
 	nstride = 3;
+
       } else if (thresh_array[ithresh] == Q) {
 	if (!atom->q_flag)
 	  error->all("Threshhold for an atom property that isn't allocated");
@@ -606,6 +608,7 @@ int DumpCustom::count()
 	  error->all("Threshhold for an atom property that isn't allocated");
 	ptr = &atom->mu[0][2];
 	nstride = 3;
+
       } else if (thresh_array[ithresh] == RADIUS) {
 	if (!atom->radius_flag)
 	  error->all("Threshhold for an atom property that isn't allocated");
@@ -641,6 +644,7 @@ int DumpCustom::count()
 	  error->all("Threshhold for an atom property that isn't allocated");
 	ptr = &atom->angmom[0][2];
 	nstride = 3;
+
       } else if (thresh_array[ithresh] == QUATW) {
 	if (!atom->quat_flag)
 	  error->all("Threshhold for an atom property that isn't allocated");
@@ -677,6 +681,29 @@ int DumpCustom::count()
 	ptr = &atom->torque[0][2];
 	nstride = 3;
 
+      } else if (thresh_array[ithresh] == SPIN) {
+	if (!atom->spin_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+        int *spin = atom->spin;
+        for (i = 0; i < nlocal; i++) dchoose[i] = spin[i];
+        ptr = dchoose;
+        nstride = 1;
+      } else if (thresh_array[ithresh] == ERADIUS) {
+	if (!atom->eradius_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+        ptr = atom->eradius;
+        nstride = 1;
+      } else if (thresh_array[ithresh] == EVEL) {
+	if (!atom->evel_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+        ptr = atom->evel;
+        nstride = 1;
+      } else if (thresh_array[ithresh] == EFORCE) {
+	if (!atom->eforce_flag)
+	  error->all("Threshhold for an atom property that isn't allocated");
+        ptr = atom->eforce;
+        nstride = 1;				
+
       } else if (thresh_array[ithresh] == COMPUTE) {
 	i = nfield + ithresh;
 	if (argindex[i] == 0) {
@@ -971,6 +998,27 @@ void DumpCustom::parse_fields(int narg, char **arg)
       pack_choice[i] = &DumpCustom::pack_tqz;
       vtype[i] = DOUBLE;
 
+    } else if (strcmp(arg[iarg],"spin") == 0) {
+      if (!atom->spin_flag)
+        error->all("Dumping an atom quantity that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_spin;
+      vtype[i] = INT;
+    } else if (strcmp(arg[iarg],"eradius") == 0) {
+      if (!atom->eradius_flag)
+        error->all("Dumping an atom quantity that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_eradius;
+      vtype[i] = DOUBLE;
+    } else if (strcmp(arg[iarg],"evel") == 0) {
+      if (!atom->evel_flag)
+        error->all("Dumping an atom quantity that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_evel;
+      vtype[i] = DOUBLE;
+    } else if (strcmp(arg[iarg],"eforce") == 0) {
+      if (!atom->evel_flag)
+        error->all("Dumping an atom quantity that isn't allocated");
+      pack_choice[i] = &DumpCustom::pack_eforce;
+      vtype[i] = DOUBLE;
+
     // compute value = c_ID
     // if no trailing [], then arg is set to 0, else arg is int between []
 
@@ -1234,10 +1282,12 @@ int DumpCustom::modify_param(int narg, char **arg)
     else if (strcmp(arg[1],"fx") == 0) thresh_array[nthresh] = FX;
     else if (strcmp(arg[1],"fy") == 0) thresh_array[nthresh] = FY;
     else if (strcmp(arg[1],"fz") == 0) thresh_array[nthresh] = FZ;
+
     else if (strcmp(arg[1],"q") == 0) thresh_array[nthresh] = Q;
     else if (strcmp(arg[1],"mux") == 0) thresh_array[nthresh] = MUX;
     else if (strcmp(arg[1],"muy") == 0) thresh_array[nthresh] = MUY;
     else if (strcmp(arg[1],"muz") == 0) thresh_array[nthresh] = MUZ;
+
     else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS;
     else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX;
     else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY;
@@ -1245,6 +1295,7 @@ int DumpCustom::modify_param(int narg, char **arg)
     else if (strcmp(arg[1],"angmomx") == 0) thresh_array[nthresh] = ANGMOMX;
     else if (strcmp(arg[1],"angmomy") == 0) thresh_array[nthresh] = ANGMOMY;
     else if (strcmp(arg[1],"angmomz") == 0) thresh_array[nthresh] = ANGMOMZ;
+
     else if (strcmp(arg[1],"quatw") == 0) thresh_array[nthresh] = QUATW;
     else if (strcmp(arg[1],"quati") == 0) thresh_array[nthresh] = QUATI;
     else if (strcmp(arg[1],"quatj") == 0) thresh_array[nthresh] = QUATJ;
@@ -1253,6 +1304,11 @@ int DumpCustom::modify_param(int narg, char **arg)
     else if (strcmp(arg[1],"tqy") == 0) thresh_array[nthresh] = TQY;
     else if (strcmp(arg[1],"tqz") == 0) thresh_array[nthresh] = TQZ;
 
+    else if (strcmp(arg[1],"spin") == 0) thresh_array[nthresh] = SPIN;
+    else if (strcmp(arg[1],"eradius") == 0) thresh_array[nthresh] = ERADIUS;
+    else if (strcmp(arg[1],"evel") == 0) thresh_array[nthresh] = EVEL;
+    else if (strcmp(arg[1],"eforce") == 0) thresh_array[nthresh] = EFORCE;
+
     // compute value = c_ID
     // if no trailing [], then arg is set to 0, else arg is between []
     // must grow field2index and argindex arrays, since access is beyond nfield
@@ -2175,3 +2231,86 @@ void DumpCustom::pack_tqz(int n)
       n += size_one;
     }
 }
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_spin(int n)
+{
+  int *spin = atom->spin;
+  int nlocal = atom->nlocal;
+  
+  for (int i = 0; i < nlocal; i++)
+    if (choose[i]) {
+      buf[n] = spin[i];
+      n += size_one;
+    }
+}
+
+/* ----------------------------------------------------------------------
+   different interpretation of electron radius
+   depending on dynamics vs minimization
+------------------------------------------------------------------------- */
+
+void DumpCustom::pack_eradius(int n)
+{
+  double *eradius = atom->eradius;
+  int *spin = atom->spin;
+  int nlocal = atom->nlocal;
+
+  if (update->whichflag == 1) {
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+	buf[n] = eradius[i];
+	n += size_one;
+      }
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+        if (spin[i]) buf[n] = exp(eradius[i]);
+	else buf[n] = 0.0;
+	n += size_one;
+      }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_evel(int n)
+{
+  double *evel = atom->evel;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (choose[i]) {
+      buf[n] = evel[i];
+      n += size_one;
+    }
+}
+
+/* ----------------------------------------------------------------------
+   different interpretation of electron radial force
+   depending on dynamics vs minimization
+------------------------------------------------------------------------- */
+
+void DumpCustom::pack_eforce(int n)
+{
+  double *eradius = atom->eradius;
+  double *eforce = atom->eforce;
+  int *spin = atom->spin;
+  int nlocal = atom->nlocal;
+  
+  if (update->whichflag == 1) {
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+	buf[n] = eforce[i];
+	n += size_one;
+      }
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+        if (spin[i]) buf[n] = eforce[i]/exp(eradius[i]);
+        else buf[n] = 0.0;
+	n += size_one;
+      }
+  }
+}

src/dump_custom.h

@@ -148,6 +148,10 @@ class DumpCustom : public Dump {
   void pack_tqx(int);
   void pack_tqy(int);
   void pack_tqz(int);
+  void pack_spin(int);
+  void pack_eradius(int);
+  void pack_evel(int);
+  void pack_eforce(int);
 
   void pack_compute(int);
   void pack_fix(int);

src/min.cpp

@@ -135,10 +135,12 @@ void Min::init()
   ev_setup();
 
   // set flags for what arrays to clear in force_clear()
-  // clear torques if array exists
+  // need to clear torques,eforce if arrays exists
 
   torqueflag = 0;
-  if (atom->torque) torqueflag = 1;
+  if (atom->torque_flag) torqueflag = 1;
+  eforceflag = 0;
+  if (atom->eforce_flag) eforceflag = 1;
 
   // orthogonal vs triclinic simulation box
 
@@ -496,6 +498,8 @@ double Min::energy_force(int resetflag)
 
 void Min::force_clear()
 {
+  int i;
+
   // clear global force array
   // nall includes ghosts only if either newton flag is set
 
@@ -504,7 +508,7 @@ void Min::force_clear()
   else nall = atom->nlocal;
 
   double **f = atom->f;
-  for (int i = 0; i < nall; i++) {
+  for (i = 0; i < nall; i++) {
     f[i][0] = 0.0;
     f[i][1] = 0.0;
     f[i][2] = 0.0;
@@ -512,12 +516,18 @@ void Min::force_clear()
 
   if (torqueflag) {
     double **torque = atom->torque;
-    for (int i = 0; i < nall; i++) {
+    for (i = 0; i < nall; i++) {
       torque[i][0] = 0.0;
       torque[i][1] = 0.0;
       torque[i][2] = 0.0;
     }
   }
+
+  if (eforceflag) {
+    double *eforce = atom->eforce;
+    for (i = 0; i < nall; i++)
+      eforce[i] = 0.0;
+  }
 }
 
 /* ----------------------------------------------------------------------

src/min.h

@@ -57,7 +57,8 @@ class Min : protected Pointers {
   class Compute **vlist_global;
   class Compute **vlist_atom;
 
-  int pairflag,torqueflag;
+  int pairflag;
+  int torqueflag,eforceflag;
   int triclinic;              // 0 if domain is orthog, 1 if triclinic
 
   int narray;                       // # of arrays stored by fix_minimize

src/update.cpp

@@ -166,6 +166,19 @@ void Update::set_units(const char *style)
     dt = 1.0e-8;
     neighbor->skin = 0.1;
 
+  } else if (strcmp(style,"electron") == 0) {
+    force->boltz = 3.16679e-6;
+    force->mvv2e = 1.0;
+    force->ftm2v = 1 / 1.0327499;	
+    force->mv2d = 1.0;
+    force->nktv2p = 2.9421e13;
+    force->qqr2e = 1.0;
+    force->qe2f = 1.0;
+    force->vxmu2f = 0.6241509647;
+    force->xxt2kmu = 1.0e-4;
+    dt = 0.001;
+    neighbor->skin = 2.0;
+
   } else error->all("Illegal units command");
 
   delete [] unit_style;

src/verlet.cpp

@@ -63,10 +63,12 @@ void Verlet::init()
   ev_setup();
 
   // set flags for what arrays to clear in force_clear()
-  // need to clear torques if array exists
+  // need to clear torques,eforce if arrays exists
 
   torqueflag = 0;
   if (atom->torque_flag) torqueflag = 1;
+  eforceflag = 0;
+  if (atom->eforce_flag) eforceflag = 1;
 
   // orthogonal vs triclinic simulation box
 
@@ -317,6 +319,12 @@ void Verlet::force_clear()
       }
     }
 
+    if (eforceflag) {
+      double *eforce = atom->eforce;
+      for (i = 0; i < nall; i++)
+	eforce[i] = 0.0;
+    }
+
   // neighbor includegroup flag is set
   // clear force only on initial nfirst particles
   // if either newton flag is set, also include ghosts
@@ -340,6 +348,12 @@ void Verlet::force_clear()
       }
     }
 
+    if (eforceflag) {
+      double *eforce = atom->eforce;
+      for (i = 0; i < nall; i++)
+	eforce[i] = 0.0;
+    }
+
     if (force->newton) {
       nall = atom->nlocal + atom->nghost;
 
@@ -357,6 +371,12 @@ void Verlet::force_clear()
 	  torque[i][2] = 0.0;
 	}
       }
+
+      if (eforceflag) {
+	double *eforce = atom->eforce;
+	for (i = atom->nlocal; i < nall; i++)
+	  eforce[i] = 0.0;
+      }
     }
   }
 }

src/verlet.h

@@ -35,7 +35,8 @@ class Verlet : public Integrate {
 
  private:
   int triclinic;                    // 0 if domain is orthog, 1 if triclinic
-  int torqueflag;                   // zero out array every step
+  int torqueflag;                   // zero out arrays every step
+  int eforceflag;
 
   void force_clear();
 };