diff --git a/src/compute_centro_atom.cpp b/src/compute_centro_atom.cpp
index 6c35ba1e8244014cd49516be139948eca11d49a6..9c3bb35e96ad548f261bc467bd1b75fd86620a22 100644
--- a/src/compute_centro_atom.cpp
+++ b/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);
- nmax = atom->nmax;
- memory->create(centro,nmax,"centro/atom:centro");
- vector_atom = centro;
+ if (!axes_flag) {
+ memory->destroy(centro);
+ nmax = atom->nmax;
+ 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,45 +204,122 @@ void ComputeCentroAtom::compute_peratom()
}
}
- // if not nnn neighbors, centro = 0.0
-
- if (n < nnn) {
- centro[i] = 0.0;
- continue;
- }
-
- // store nnn nearest neighs in 1st nnn locations of distsq and nearest
-
- select2(nnn,n,distsq,nearest);
-
- // R = Ri + Rj for each of npairs i,j pairs among nnn neighbors
- // pairs = squared length of each R
-
- n = 0;
- 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;
- pairs[n++] = delx*delx + dely*dely + delz*delz;
- }
+ // check whether to include local crystal symmetry axes
+
+ if (!axes_flag) {
+
+ // if not nnn neighbors, centro = 0.0
+
+ if (n < nnn) {
+ centro[i] = 0.0;
+ continue;
+ }
+
+ // store nnn nearest neighs in 1st nnn locations of distsq and nearest
+
+ select2(nnn,n,distsq,nearest);
+
+ // R = Ri + Rj for each of npairs i,j pairs among nnn neighbors
+ // pairs = squared length of each R
+
+ n = 0;
+ 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;
+ 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
-
+
select(nhalf,npairs,pairs);
-
+
// centrosymmetry = sum of nhalf smallest squared values
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];
+ }
}
/* ----------------------------------------------------------------------
diff --git a/src/compute_centro_atom.h b/src/compute_centro_atom.h
index 6575e3fb5808d587ddd7525f0f76769ac166757d..ef3128ee68e5ac9ba89e8a25d833b5d73a1fb553 100644
--- a/src/compute_centro_atom.h
+++ b/src/compute_centro_atom.h
@@ -39,7 +39,8 @@ 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 *);
};
diff --git a/src/math_extra.h b/src/math_extra.h
index 4e6ee4c77c3886ba2e73552e83f61dd66c3b1170..0bfc7193a4e31622af2b58db0a0c3b2af97e6d64 100755
--- a/src/math_extra.h
+++ b/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
------------------------------------------------------------------------- */