From 92d15d4a8931286cb53ea321f4bcd4cabc579e3d Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 10 Jan 2017 12:52:37 -0500
Subject: [PATCH] replace string compare with enums, fix memory leak,
formatting cleanup
---
src/compute_coord_atom.cpp | 127 +++++++++++++-------------
src/compute_coord_atom.h | 5 +-
src/compute_orientorder_atom.cpp | 150 ++++++++++++++++---------------
src/compute_orientorder_atom.h | 2 +-
4 files changed, 146 insertions(+), 138 deletions(-)
diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp
index 1ad553ee37..36f0b63504 100644
--- a/src/compute_coord_atom.cpp
+++ b/src/compute_coord_atom.cpp
@@ -36,16 +36,15 @@ using namespace LAMMPS_NS;
ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
- cstyle(NULL), id_orientorder(NULL), typelo(NULL), typehi(NULL), cvec(NULL), carray(NULL)
+ typelo(NULL), typehi(NULL), cvec(NULL), carray(NULL),
+ id_orientorder(NULL), normv(NULL)
{
if (narg < 5) error->all(FLERR,"Illegal compute coord/atom command");
- int n = strlen(arg[3]) + 1;
- cstyle = new char[n];
- strcpy(cstyle,arg[3]);
-
- if (strcmp(cstyle,"cutoff") == 0) {
+ cstyle = NONE;
+ if (strcmp(arg[3],"cutoff") == 0) {
+ cstyle = CUTOFF;
double cutoff = force->numeric(FLERR,arg[4]);
cutsq = cutoff*cutoff;
@@ -68,34 +67,33 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
ncol++;
iarg++;
}
-
}
- } else if (strcmp(cstyle,"orientorder") == 0) {
+ } else if (strcmp(arg[3],"orientorder") == 0) {
+ cstyle = ORIENT;
+ if (narg != 6) error->all(FLERR,"Illegal compute coord/atom command");
- if (narg != 6) error->all(FLERR,"Illegal compute coord/atom command");
+ int n = strlen(arg[4]) + 1;
+ id_orientorder = new char[n];
+ strcpy(id_orientorder,arg[4]);
- n = strlen(arg[4]) + 1;
- id_orientorder = new char[n];
- strcpy(id_orientorder,arg[4]);
-
- int iorientorder = modify->find_compute(id_orientorder);
- if (iorientorder < 0)
- error->all(FLERR,"Could not find compute coord/atom compute ID");
- if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0)
- error->all(FLERR,"Compute coord/atom compute ID does not compute orientorder/atom");
+ int iorientorder = modify->find_compute(id_orientorder);
+ if (iorientorder < 0)
+ error->all(FLERR,"Could not find compute coord/atom compute ID");
+ if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0)
+ error->all(FLERR,"Compute coord/atom compute ID does not compute orientorder/atom");
- threshold = force->numeric(FLERR,arg[5]);
- if (threshold <= -1.0 || threshold >= 1.0)
- error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1");
+ threshold = force->numeric(FLERR,arg[5]);
+ if (threshold <= -1.0 || threshold >= 1.0)
+ error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1");
- ncol = 1;
- typelo = new int[ncol];
- typehi = new int[ncol];
- typelo[0] = 1;
- typehi[0] = atom->ntypes;
+ ncol = 1;
+ typelo = new int[ncol];
+ typehi = new int[ncol];
+ typelo[0] = 1;
+ typehi[0] = atom->ntypes;
- } else error->all(FLERR,"Invalid cstyle in compute coord/atom");
+ } else error->all(FLERR,"Invalid cstyle in compute coord/atom");
peratom_flag = 1;
if (ncol == 1) size_peratom_cols = 0;
@@ -112,21 +110,23 @@ ComputeCoordAtom::~ComputeCoordAtom()
delete [] typehi;
memory->destroy(cvec);
memory->destroy(carray);
+ delete [] id_orientorder;
}
/* ---------------------------------------------------------------------- */
void ComputeCoordAtom::init()
{
- if (strcmp(cstyle,"orientorder") == 0) {
+ if (cstyle == ORIENT) {
int iorientorder = modify->find_compute(id_orientorder);
c_orientorder = (ComputeOrientOrderAtom*)(modify->compute[iorientorder]);
cutsq = c_orientorder->cutsq;
l = c_orientorder->qlcomp;
- // communicate real and imaginary 2*l+1 components of the normalized vector
+ // communicate real and imaginary 2*l+1 components of the normalized vector
comm_forward = 2*(2*l+1);
if (c_orientorder->iqlcomp < 0)
- error->all(FLERR,"Compute coord/atom requires components option in compute orientorder/atom be defined");
+ error->all(FLERR,"Compute coord/atom requires components "
+ "option in compute orientorder/atom be defined");
}
if (force->pair == NULL)
@@ -188,7 +188,7 @@ void ComputeCoordAtom::compute_peratom()
}
}
- if (strcmp(cstyle,"orientorder") == 0) {
+ if (cstyle == ORIENT) {
if (!(c_orientorder->invoked_flag & INVOKED_PERATOM)) {
c_orientorder->compute_peratom();
c_orientorder->invoked_flag |= INVOKED_PERATOM;
@@ -217,7 +217,7 @@ void ComputeCoordAtom::compute_peratom()
int *type = atom->type;
int *mask = atom->mask;
- if (strcmp(cstyle,"cutoff") == 0) {
+ if (cstyle == CUTOFF) {
if (ncol == 1) {
@@ -241,7 +241,7 @@ void ComputeCoordAtom::compute_peratom()
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0])
- n++;
+ n++;
}
cvec[i] = n;
@@ -281,37 +281,36 @@ void ComputeCoordAtom::compute_peratom()
}
}
- } else if (strcmp(cstyle,"orientorder") == 0) {
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- if (mask[i] & groupbit) {
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- n = 0;
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
-
- delx = xtmp - x[j][0];
- dely = ytmp - x[j][1];
- delz = ztmp - x[j][2];
- rsq = delx*delx + dely*dely + delz*delz;
- if (rsq < cutsq) {
- double dot_product = 0.0;
- for (int m=0; m < 2*(2*l+1); m++) {
- dot_product += normv[i][nqlist+m]*normv[j][nqlist+m];
- }
- if (dot_product > threshold) n++;
- }
- }
- cvec[i] = n;
- } else cvec[i] = 0.0;
- }
+ } else if (cstyle == ORIENT) {
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ if (mask[i] & groupbit) {
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ n = 0;
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ if (rsq < cutsq) {
+ double dot_product = 0.0;
+ for (int m=0; m < 2*(2*l+1); m++) {
+ dot_product += normv[i][nqlist+m]*normv[j][nqlist+m];
+ }
+ if (dot_product > threshold) n++;
+ }
+ }
+ cvec[i] = n;
+ } else cvec[i] = 0.0;
+ }
}
}
diff --git a/src/compute_coord_atom.h b/src/compute_coord_atom.h
index 942882179f..2ad46fa854 100644
--- a/src/compute_coord_atom.h
+++ b/src/compute_coord_atom.h
@@ -34,6 +34,7 @@ class ComputeCoordAtom : public Compute {
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
double memory_usage();
+ enum {NONE,CUTOFF,ORIENT};
private:
int nmax,ncol;
@@ -45,10 +46,10 @@ class ComputeCoordAtom : public Compute {
double **carray;
class ComputeOrientOrderAtom *c_orientorder;
- char *cstyle,*id_orientorder;
+ char *id_orientorder;
double threshold;
double **normv;
- int nqlist,l;
+ int cstyle,nqlist,l;
};
}
diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp
index 42186f03b8..5f78b33b61 100644
--- a/src/compute_orientorder_atom.cpp
+++ b/src/compute_orientorder_atom.cpp
@@ -46,7 +46,8 @@ using namespace std;
ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
- distsq(NULL), nearest(NULL), rlist(NULL), qlist(NULL), qnarray(NULL), qnm_r(NULL), qnm_i(NULL)
+ qlist(NULL), distsq(NULL), nearest(NULL), rlist(NULL),
+ qnarray(NULL), qnm_r(NULL), qnm_i(NULL)
{
if (narg < 3 ) error->all(FLERR,"Illegal compute orientorder/atom command");
@@ -57,7 +58,7 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
qlcompflag = 0;
// specify which orders to request
-
+
nqlist = 5;
memory->create(qlist,nqlist,"orientorder/atom:qlist");
qlist[0] = 4;
@@ -73,48 +74,55 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
while (iarg < narg) {
if (strcmp(arg[iarg],"nnn") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
- if (strcmp(arg[iarg+1],"NULL") == 0)
- nnn = 0;
- else {
- nnn = force->numeric(FLERR,arg[iarg+1]);
- if (nnn <= 0)
- error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (strcmp(arg[iarg+1],"NULL") == 0) {
+ nnn = 0;
+ } else {
+ nnn = force->numeric(FLERR,arg[iarg+1]);
+ if (nnn <= 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
}
iarg += 2;
} else if (strcmp(arg[iarg],"degrees") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (iarg+2 > narg)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
nqlist = force->numeric(FLERR,arg[iarg+1]);
- if (nqlist <= 0) error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (nqlist <= 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
memory->destroy(qlist);
memory->create(qlist,nqlist,"orientorder/atom:qlist");
iarg += 2;
if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
qmax = 0;
for (int iw = 0; iw < nqlist; iw++) {
- qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
- if (qlist[iw] < 0)
- error->all(FLERR,"Illegal compute orientorder/atom command");
- if (qlist[iw] > qmax) qmax = qlist[iw];
+ qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
+ if (qlist[iw] < 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (qlist[iw] > qmax) qmax = qlist[iw];
}
iarg += nqlist;
if (strcmp(arg[iarg],"components") == 0) {
- qlcompflag = 1;
- if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+ qlcompflag = 1;
+ if (iarg+2 > narg)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
qlcomp = force->numeric(FLERR,arg[iarg+1]);
- if (qlcomp <= 0) error->all(FLERR,"Illegal compute orientorder/atom command");
- iqlcomp = -1;
+ if (qlcomp <= 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
+ iqlcomp = -1;
for (int iw = 0; iw < nqlist; iw++)
if (qlcomp == qlist[iw]) {
- iqlcomp = iw;
- break;
- }
- if (iqlcomp < 0) error->all(FLERR,"Illegal compute orientorder/atom command");
+ iqlcomp = iw;
+ break;
+ }
+ if (iqlcomp < 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
iarg += 2;
}
} else if (strcmp(arg[iarg],"cutoff") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (iarg+2 > narg)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
double cutoff = force->numeric(FLERR,arg[iarg+1]);
- if (cutoff <= 0.0) error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (cutoff <= 0.0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
cutsq = cutoff*cutoff;
iarg += 2;
} else error->all(FLERR,"Illegal compute orientorder/atom command");
@@ -141,7 +149,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
memory->destroy(qlist);
memory->destroy(qnm_r);
memory->destroy(qnm_i);
-
+
}
/* ---------------------------------------------------------------------- */
@@ -224,7 +232,7 @@ void ComputeOrientOrderAtom::compute_peratom()
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
-
+
// insure distsq and nearest arrays are long enough
if (jnum > maxneigh) {
@@ -253,9 +261,9 @@ void ComputeOrientOrderAtom::compute_peratom()
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq) {
distsq[ncount] = rsq;
- rlist[ncount][0] = delx;
- rlist[ncount][1] = dely;
- rlist[ncount][2] = delz;
+ rlist[ncount][0] = delx;
+ rlist[ncount][1] = dely;
+ rlist[ncount][2] = delz;
nearest[ncount++] = j;
}
}
@@ -263,16 +271,16 @@ void ComputeOrientOrderAtom::compute_peratom()
// if not nnn neighbors, order parameter = 0;
if ((ncount == 0) || (ncount < nnn)) {
- for (int iw = 0; iw < nqlist; iw++)
- qn[iw] = 0.0;
+ for (int iw = 0; iw < nqlist; iw++)
+ qn[iw] = 0.0;
continue;
}
// if nnn > 0, use only nearest nnn neighbors
if (nnn > 0) {
- select3(nnn,ncount,distsq,nearest,rlist);
- ncount = nnn;
+ select3(nnn,ncount,distsq,nearest,rlist);
+ ncount = nnn;
}
calc_boop(rlist, ncount, qn, qlist, nqlist);
@@ -287,8 +295,8 @@ void ComputeOrientOrderAtom::compute_peratom()
double ComputeOrientOrderAtom::memory_usage()
{
double bytes = ncol*nmax * sizeof(double);
- bytes += (qmax*(2*qmax+1)+maxneigh*4) * sizeof(double);
- bytes += (nqlist+maxneigh) * sizeof(int);
+ bytes += (qmax*(2*qmax+1)+maxneigh*4) * sizeof(double);
+ bytes += (nqlist+maxneigh) * sizeof(int);
return bytes;
}
@@ -300,18 +308,18 @@ double ComputeOrientOrderAtom::memory_usage()
// Use no-op do while to create single statement
-#define SWAP(a,b) do { \
- tmp = a; a = b; b = tmp; \
+#define SWAP(a,b) do { \
+ tmp = a; a = b; b = tmp; \
} while(0)
-#define ISWAP(a,b) do { \
- itmp = a; a = b; b = itmp; \
+#define ISWAP(a,b) do { \
+ itmp = a; a = b; b = itmp; \
} while(0)
-#define SWAP3(a,b) do { \
- tmp = a[0]; a[0] = b[0]; b[0] = tmp; \
- tmp = a[1]; a[1] = b[1]; b[1] = tmp; \
- tmp = a[2]; a[2] = b[2]; b[2] = tmp; \
+#define SWAP3(a,b) do { \
+ tmp = a[0]; a[0] = b[0]; b[0] = tmp; \
+ tmp = a[1]; a[1] = b[1]; b[1] = tmp; \
+ tmp = a[2]; a[2] = b[2]; b[2] = tmp; \
} while(0)
/* ---------------------------------------------------------------------- */
@@ -330,7 +338,7 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir]);
- ISWAP(iarr[l],iarr[ir]);
+ ISWAP(iarr[l],iarr[ir]);
SWAP3(arr3[l],arr3[ir]);
}
return;
@@ -342,17 +350,17 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir]);
ISWAP(iarr[l],iarr[ir]);
- SWAP3(arr3[l],arr3[ir]);
+ SWAP3(arr3[l],arr3[ir]);
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir]);
ISWAP(iarr[l+1],iarr[ir]);
- SWAP3(arr3[l+1],arr3[ir]);
+ SWAP3(arr3[l+1],arr3[ir]);
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1]);
ISWAP(iarr[l],iarr[l+1]);
- SWAP3(arr3[l],arr3[l+1]);
+ SWAP3(arr3[l],arr3[l+1]);
}
i = l+1;
j = ir;
@@ -367,7 +375,7 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (j < i) break;
SWAP(arr[i],arr[j]);
ISWAP(iarr[i],iarr[j]);
- SWAP3(arr3[i],arr3[j]);
+ SWAP3(arr3[i],arr3[j]);
}
arr[l+1] = arr[j];
arr[j] = a;
@@ -389,9 +397,9 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
calculate the bond orientational order parameters
------------------------------------------------------------------------- */
-void ComputeOrientOrderAtom::calc_boop(double **rlist,
- int ncount, double qn[],
- int qlist[], int nqlist) {
+void ComputeOrientOrderAtom::calc_boop(double **rlist,
+ int ncount, double qn[],
+ int qlist[], int nqlist) {
for (int iw = 0; iw < nqlist; iw++) {
int n = qlist[iw];
@@ -429,22 +437,22 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
double expphim_r = expphi_r;
double expphim_i = expphi_i;
for(int m = 1; m <= +n; m++) {
- double prefactor = polar_prefactor(n, m, costheta);
- double c_r = prefactor * expphim_r;
- double c_i = prefactor * expphim_i;
- qnm_r[iw][m+n] += c_r;
- qnm_i[iw][m+n] += c_i;
- if(m & 1) {
- qnm_r[iw][-m+n] -= c_r;
- qnm_i[iw][-m+n] += c_i;
- } else {
- qnm_r[iw][-m+n] += c_r;
- qnm_i[iw][-m+n] -= c_i;
- }
- double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
- double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
- expphim_r = tmp_r;
- expphim_i = tmp_i;
+ double prefactor = polar_prefactor(n, m, costheta);
+ double c_r = prefactor * expphim_r;
+ double c_i = prefactor * expphim_i;
+ qnm_r[iw][m+n] += c_r;
+ qnm_i[iw][m+n] += c_i;
+ if(m & 1) {
+ qnm_r[iw][-m+n] -= c_r;
+ qnm_i[iw][-m+n] += c_i;
+ } else {
+ qnm_r[iw][-m+n] += c_r;
+ qnm_i[iw][-m+n] -= c_i;
+ }
+ double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
+ double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
+ expphim_r = tmp_r;
+ expphim_i = tmp_i;
}
}
@@ -458,15 +466,15 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
for(int m = 0; m < 2*n+1; m++) {
qm_sum += qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m];
// printf("Ylm^2 = %d %d %g\n",n,m,
- // qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m]);
+ // qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m]);
}
qn[iw] = fac * sqrt(qm_sum / (2*n+1));
if (qlcompflag && iqlcomp == iw) normfac = 1.0/sqrt(qm_sum);
}
-
+
// output of the complex vector
-
+
if (qlcompflag) {
int j = nqlist;
for(int m = 0; m < 2*qlcomp+1; m++) {
@@ -485,7 +493,7 @@ double ComputeOrientOrderAtom::dist(const double r[]) {
}
/* ----------------------------------------------------------------------
- polar prefactor for spherical harmonic Y_l^m, where
+ polar prefactor for spherical harmonic Y_l^m, where
Y_l^m (theta, phi) = prefactor(l, m, cos(theta)) * exp(i*m*phi)
------------------------------------------------------------------------- */
diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h
index d6f32b7727..81b08dbddc 100644
--- a/src/compute_orientorder_atom.h
+++ b/src/compute_orientorder_atom.h
@@ -49,7 +49,7 @@ class ComputeOrientOrderAtom : public Compute {
double **qnm_i;
void select3(int, int, double *, int *, double **);
- void calc_boop(double **rlist, int numNeighbors,
+ void calc_boop(double **rlist, int numNeighbors,
double qn[], int nlist[], int nnlist);
double dist(const double r[]);
--
GitLab