Added axes keyword

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

src/compute_centro_atom.cpp

@@ -27,6 +27,7 @@
 #include "force.h"
 #include "pair.h"
 #include "comm.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "error.h"
 
@@ -37,23 +38,42 @@ using namespace LAMMPS_NS;
 ComputeCentroAtom::ComputeCentroAtom(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
-  if (narg != 4) error->all(FLERR,"Illegal compute centro/atom command");
+  if (narg < 4 || narg > 6) error->all(FLERR,"Illegal compute centro/atom command");
 
   if (strcmp(arg[3],"fcc") == 0) nnn = 12;
   else if (strcmp(arg[3],"bcc") == 0) nnn = 8;
   else nnn = force->inumeric(FLERR,arg[3]);
 
+  // default values
+
+  axes_flag = 0;
+  
+  // optional keywords
+  
+  int iarg = 4;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"axes") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute centro/atom command3");
+      if (strcmp(arg[iarg+1],"yes") == 0) axes_flag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) axes_flag = 0;
+      else error->all(FLERR,"Illegal compute centro/atom command2");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute centro/atom command1");
+  }
+
   if (nnn <= 0 || nnn % 2)
     error->all(FLERR,"Illegal neighbor value for compute centro/atom command");
 
   peratom_flag = 1;
-  size_peratom_cols = 0;
+  if (!axes_flag) size_peratom_cols = 0;
+  else size_peratom_cols = 10;
   
   nmax = 0;
   centro = NULL;
   maxneigh = 0;
   distsq = NULL;
   nearest = NULL;
+  array_atom = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -63,6 +83,7 @@ ComputeCentroAtom::~ComputeCentroAtom()
   memory->destroy(centro);
   memory->destroy(distsq);
   memory->destroy(nearest);
+  if (axes_flag) memory->destroy(array_atom);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -106,12 +127,21 @@ void ComputeCentroAtom::compute_peratom()
   invoked_peratom = update->ntimestep;
 
   // grow centro array if necessary
+  // grow array_atom array if axes_flag set
 
   if (atom->nmax > nmax) {
-    memory->destroy(centro);
+    if (!axes_flag) {
-    nmax = atom->nmax;
+      memory->destroy(centro);
-    memory->create(centro,nmax,"centro/atom:centro");
+      nmax = atom->nmax;
-    vector_atom = centro;
+      memory->create(centro,nmax,"centro/atom:centro");
+      vector_atom = centro;
+    } else {
+      memory->destroy(centro);
+      memory->destroy(array_atom);
+      nmax = atom->nmax;
+      memory->create(centro,nmax,"centro/atom:centro");
+      memory->create(array_atom,nmax,size_peratom_cols,"centro/atom:array_atom");
+    }
   }
 
   // invoke full neighbor list (will copy or build if necessary)
@@ -174,30 +204,92 @@ void ComputeCentroAtom::compute_peratom()
         }
       }
 
-      // if not nnn neighbors, centro = 0.0
+      // check whether to include local crystal symmetry axes
       
-      if (n < nnn) {
+      if (!axes_flag) {
-        centro[i] = 0.0;
-        continue;
-      }
 	
-      // store nnn nearest neighs in 1st nnn locations of distsq and nearest
+	// if not nnn neighbors, centro = 0.0
 
-      select2(nnn,n,distsq,nearest);
+	if (n < nnn) {
+	  centro[i] = 0.0;
+	  continue;
+	}
 
-      // R = Ri + Rj for each of npairs i,j pairs among nnn neighbors
+	// store nnn nearest neighs in 1st nnn locations of distsq and nearest
-      // pairs = squared length of each R
 
-      n = 0;
+	select2(nnn,n,distsq,nearest);
-      for (j = 0; j < nnn; j++) {
+
-        jj = nearest[j];
+	// R = Ri + Rj for each of npairs i,j pairs among nnn neighbors
-        for (k = j+1; k < nnn; k++) {
+	// pairs = squared length of each R
-          kk = nearest[k];
+
-          delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
+	n = 0;
-          dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
+	for (j = 0; j < nnn; j++) {
-          delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
+	  jj = nearest[j];
-          pairs[n++] = delx*delx + dely*dely + delz*delz;
+	  for (k = j+1; k < nnn; k++) {
-        }
+	    kk = nearest[k];
+	    delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
+	    dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
+	    delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
+	    pairs[n++] = delx*delx + dely*dely + delz*delz;
+
+	  }
+	}
+	
+      } else {
+	
+	// calculate local crystal symmetry axes
+
+	// rsq1, rsq2 are two smallest values of R^2 
+	// R1, R2 are corresponding vectors Ri - Rj
+	// R3 is normal to R1, R2
+
+	double rsq1,rsq2;
+
+	double* r1 = &array_atom[i][1];
+	double* r2 = &array_atom[i][4];
+	double* r3 = &array_atom[i][7];
+	
+	if (n < nnn) {
+	  centro[i] = 0.0;
+	  MathExtra::zero3(r1);
+	  MathExtra::zero3(r2);
+	  MathExtra::zero3(r3);
+	  continue;
+	}
+
+	// store nnn nearest neighs in 1st nnn locations of distsq and nearest
+
+	select2(nnn,n,distsq,nearest);
+
+	n = 0;
+	rsq1 = rsq2 = cutsq;
+	for (j = 0; j < nnn; j++) {
+	  jj = nearest[j];
+	  for (k = j+1; k < nnn; k++) {
+	    kk = nearest[k];
+	    delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
+	    dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
+	    delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
+	    double rsq = delx*delx + dely*dely + delz*delz;
+	    pairs[n++] = rsq;
+	    
+	    if (rsq < rsq2)
+	      if (rsq < rsq1) {
+		rsq2 = rsq1;
+		MathExtra::copy3(r1, r2);
+		rsq1 = rsq;
+		MathExtra::sub3(x[jj],x[kk],r1);
+	      } else {
+		rsq2 = rsq;
+		MathExtra::sub3(x[jj],x[kk],r2);
+	      }
+	  }
+	}
+
+	MathExtra::cross3(r1,r2,r3);
+	MathExtra::norm3(r1);
+	MathExtra::norm3(r2);
+	MathExtra::norm3(r3);
       }
       
       // store nhalf smallest pair distances in 1st nhalf locations of pairs
@@ -209,10 +301,25 @@ void ComputeCentroAtom::compute_peratom()
       value = 0.0;
       for (j = 0; j < nhalf; j++) value += pairs[j];
       centro[i] = value;
-    } else centro[i] = 0.0;
+
+    } else {
+      centro[i] = 0.0;
+      if (axes_flag) {
+	MathExtra::zero3(&array_atom[i][1]);
+	MathExtra::zero3(&array_atom[i][4]);
+	MathExtra::zero3(&array_atom[i][7]);
+      }
+    }
   }
 
   delete [] pairs;
+
+  if (axes_flag)
+    for (ii = 0; ii < inum; ii++) {
+      i = ilist[ii];
+      if (mask[i] & groupbit)
+	array_atom[i][0] = centro[i];
+    }
 }
 
 /* ----------------------------------------------------------------------

src/compute_centro_atom.h

@@ -39,6 +39,7 @@ class ComputeCentroAtom : public Compute {
   int *nearest;
   class NeighList *list;
   double *centro;
+  int axes_flag;
   
   void select(int, int, double *);
   void select2(int, int, double *, int *);

src/math_extra.h

@@ -27,12 +27,15 @@ namespace MathExtra {
 
   // 3 vector operations
 
+  inline void copy3(const double *v, double *ans);
+  inline void zero3(double *v);
   inline void norm3(double *v);
   inline void normalize3(const double *v, double *ans);
   inline void snormalize3(const double, const double *v, double *ans);
   inline void negate3(double *v);
   inline void scale3(double s, double *v);
   inline void add3(const double *v1, const double *v2, double *ans);
+  inline void scaleadd3(double s, const double *v1, const double *v2, double *ans);
   inline void sub3(const double *v1, const double *v2, double *ans);
   inline double len3(const double *v);
   inline double lensq3(const double *v);
@@ -128,6 +131,28 @@ namespace MathExtra {
                         double *inertia);
 }
 
+/* ----------------------------------------------------------------------
+   copy a vector, return in ans
+------------------------------------------------------------------------- */
+
+inline void MathExtra::copy3(const double *v, double *ans)
+{
+  ans[0] = v[0];
+  ans[1] = v[1];
+  ans[2] = v[2];
+}
+
+/* ----------------------------------------------------------------------
+   set vector equal to zero
+------------------------------------------------------------------------- */
+
+inline void MathExtra::zero3(double *v)
+{
+  v[0] = 0.0;
+  v[1] = 0.0;
+  v[2] = 0.0;
+}
+
 /* ----------------------------------------------------------------------
    normalize a vector in place
 ------------------------------------------------------------------------- */
@@ -198,6 +223,17 @@ inline void MathExtra::add3(const double *v1, const double *v2, double *ans)
   ans[2] = v1[2] + v2[2];
 }
 
+/* ----------------------------------------------------------------------
+   ans = s*v1 + v2
+------------------------------------------------------------------------- */
+
+inline void MathExtra::scaleadd3(double s, const double *v1, const double *v2, double *ans)
+{
+  ans[0] = s*v1[0] + v2[0];
+  ans[1] = s*v1[1] + v2[1];
+  ans[2] = s*v1[2] + v2[2];
+}
+
 /* ----------------------------------------------------------------------
    ans = v1 - v2
 ------------------------------------------------------------------------- */