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

src/atom.cpp

@@ -76,6 +76,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   spin = NULL;
   eradius = ervel = erforce = NULL;
 
+  length = theta = NULL;
+
   maxspecial = 1;
   nspecial = NULL;
   special = NULL;
@@ -102,6 +104,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
   quat_flag = omega_flag = angmom_flag = torque_flag = 0;
   radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
   spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
+  line_flag = 0;
 
   // ntype-length arrays
 
@@ -185,6 +188,9 @@ Atom::~Atom()
   memory->sfree(ervel);
   memory->sfree(erforce);
 
+  memory->sfree(length);
+  memory->sfree(theta);
+
   memory->sfree(molecule);
 
   memory->destroy_2d_int_array(nspecial);
@@ -262,7 +268,8 @@ void Atom::create_avec(const char *style, int narg, char **arg)
   quat_flag = omega_flag = angmom_flag = torque_flag = 0;
   radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
   spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
-
+  line_flag = 0;
+    
   avec = new_avec(style,narg,arg);
   int n = strlen(style) + 1;
   atom_style = new char[n];

src/atom.h

@@ -55,6 +55,8 @@ class Atom : protected Pointers {
   int *spin;
   double *eradius,*ervel,*erforce;
 
+  double *length,*theta;
+
   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]
@@ -84,6 +86,7 @@ class Atom : protected Pointers {
   int quat_flag,omega_flag,angmom_flag,torque_flag;
   int radius_flag,density_flag,rmass_flag,vfrac_flag;
   int spin_flag,eradius_flag,ervel_flag,erforce_flag;
+  int line_flag;
 
   // extra peratom info in restart file destined for fix & diag 
 

src/compute_erotate_sphere.cpp

@@ -23,7 +23,8 @@
 
 using namespace LAMMPS_NS;
 
-#define INERTIA 0.4          // moment of inertia for sphere
+//#define INERTIA 0.4          // moment of inertia for sphere
+#define INERTIA (1.0/12.0)          // moment of inertia for sphere
 
 /* ---------------------------------------------------------------------- */
 
@@ -37,11 +38,13 @@ ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) :
 
   // error checks
 
+  /*
   if (!atom->omega_flag) 
     error->all("Compute erotate/sphere requires atom attribute omega");
   if (!atom->radius_flag && !atom->avec->shape_type)
     error->all("Compute erotate/sphere requires atom attribute "
 	       "radius or shape");
+  */
 }
 
 /* ---------------------------------------------------------------------- */
@@ -53,6 +56,7 @@ void ComputeERotateSphere::init()
   // if shape used, check that all particles are spherical
   // point particles are allowed
 
+  /*
   if (atom->radius == NULL) {
     double **shape = atom->shape;
     int *type = atom->type;
@@ -68,6 +72,7 @@ void ComputeERotateSphere::init()
 		     "spherical particle shapes");
       }
   }
+  */
 
   pfactor = 0.5 * force->mvv2e * INERTIA;
 }
@@ -112,6 +117,13 @@ double ComputeERotateSphere::compute_scalar()
     
   } else {
     if (rmass) {
+      for (i = 0; i < nlocal; i++) 
+	if (mask[i] & groupbit) {
+	  erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
+		      omega[i][2]*omega[i][2]) * 
+	    atom->length[i]*atom->length[i]*rmass[i];
+	}
+    } else if (rmass) {
       for (i = 0; i < nlocal; i++) 
 	if (mask[i] & groupbit) {
 	  itype = type[i];

src/compute_property_atom.cpp

@@ -11,6 +11,7 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "string.h"
 #include "compute_property_atom.h"
 #include "atom.h"
@@ -218,6 +219,37 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
 		   "atom property that isn't allocated");
       pack_choice[i] = &ComputePropertyAtom::pack_erforce;
 
+    } else if (strcmp(arg[iarg],"end1x") == 0) {
+      if (!atom->line_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_end1x;
+    } else if (strcmp(arg[iarg],"end1y") == 0) {
+      if (!atom->line_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_end1y;
+    } else if (strcmp(arg[iarg],"end1z") == 0) {
+      if (!atom->line_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_end1z;
+    } else if (strcmp(arg[iarg],"end2x") == 0) {
+      if (!atom->line_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_end2x;
+    } else if (strcmp(arg[iarg],"end2y") == 0) {
+      if (!atom->line_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_end2y;
+    } else if (strcmp(arg[iarg],"end2z") == 0) {
+      if (!atom->line_flag)
+	error->all("Compute property/atom for "
+		   "atom property that isn't allocated");
+      pack_choice[i] = &ComputePropertyAtom::pack_end2z;
+
     } else error->all("Invalid keyword in compute property/atom command");
   }
 
@@ -1100,3 +1132,102 @@ void ComputePropertyAtom::pack_erforce(int n)
     n += nvalues;
   }
 }
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end1x(int n)
+{
+  double **x = atom->x;
+  double *length = atom->length;
+  double *theta = atom->theta;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = x[i][0] - 0.5*length[i]*cos(theta[i]);
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end1y(int n)
+{
+  double **x = atom->x;
+  double *length = atom->length;
+  double *theta = atom->theta;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = x[i][1] - 0.5*length[i]*sin(theta[i]);
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end1z(int n)
+{
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = x[i][2];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end2x(int n)
+{
+  double **x = atom->x;
+  double *length = atom->length;
+  double *theta = atom->theta;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = x[i][0] + 0.5*length[i]*cos(theta[i]);
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end2y(int n)
+{
+  double **x = atom->x;
+  double *length = atom->length;
+  double *theta = atom->theta;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = x[i][1] + 0.5*length[i]*sin(theta[i]);
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end2z(int n)
+{
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = x[i][2];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+

src/compute_property_atom.h

@@ -94,6 +94,12 @@ class ComputePropertyAtom : public Compute {
   void pack_eradius(int);
   void pack_ervel(int);
   void pack_erforce(int);
+  void pack_end1x(int);
+  void pack_end1y(int);
+  void pack_end1z(int);
+  void pack_end2x(int);
+  void pack_end2y(int);
+  void pack_end2z(int);
 };
 
 }

src/input.cpp

@@ -333,13 +333,14 @@ void Input::parse()
 }
 
 /* ----------------------------------------------------------------------
-   substitute for $ variables in str
+   substitute for $ variables in str and return it
+   str assumed to be long enough to hold expanded version
    print updated string if flag is set and not searching for label
 ------------------------------------------------------------------------- */
 
 void Input::substitute(char *str, int flag)
 {
-  // use work[] as scratch space to expand str
+  // use work[] as scratch space to expand str, then copy back to str
   // do not replace $ inside single/double quotes
   // var = pts at variable name, ended by NULL
   //   if $ is followed by '{', trailing '}' becomes NULL
@@ -528,12 +529,15 @@ void Input::ifthenelse()
 
   // substitute for variables in Boolean expression for "if"
   // in case expression was enclosed in quotes
+  // must substitute on copy of arg else will step on subsequent args
 
-  substitute(arg[0],0);
+  char *scopy = new char[MAXLINE];
+  strcpy(scopy,arg[0]);
+  substitute(scopy,0);
 
   // evaluate Boolean expression for "if"
 
-  double btest = variable->evaluate_boolean(arg[0]);
+  double btest = variable->evaluate_boolean(scopy);
 
   // bound "then" commands
 
@@ -546,6 +550,8 @@ void Input::ifthenelse()
     iarg++;
   int last = iarg-1;
 
+  printf("FIRSTLAST %d %d\n",first,last);
+
   // execute "then" commands
   // make copies of all arg string commands
   // required because re-parsing a command via one() will wipe out args
@@ -569,23 +575,29 @@ void Input::ifthenelse()
     
     for (int i = 0; i < ncommands; i++) delete [] commands[i];
     delete [] commands;
+    delete [] scopy;
 
     return;
   }
 
   // done if no "elif" or "else"
 
-  if (iarg == narg) return;
+  if (iarg == narg) {
+    delete [] scopy;
+    return;
+  }
 
   // check "elif" or "else" until find commands to execute
   // substitute for variables and evaluate Boolean expression for "elif"
+  // must substitute on copy of arg else will step on subsequent args
   // bound and execute "elif" or "else" commands
 
   while (1) {
     if (iarg+2 > narg) error->all("Illegal if command");
     if (strcmp(arg[iarg],"elif") == 0) {
-      substitute(arg[iarg+1],0);
+      strcpy(scopy,arg[iarg+1]);
-      btest = variable->evaluate_boolean(arg[iarg+1]);
+      substitute(scopy,0);
+      btest = variable->evaluate_boolean(scopy);
       first = iarg+2;
     } else {
       btest = 1.0;
@@ -622,6 +634,7 @@ void Input::ifthenelse()
 
     for (int i = 0; i < ncommands; i++) delete [] commands[i];
     delete [] commands;
+    delete [] scopy;
 
     return;
   }