jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

replace string compare with enums, fix memory leak, formatting cleanup

This commit is contained in:
Axel Kohlmeyer 2017-01-10 12:52:37 -05:00
parent 95706ac846
commit 92d15d4a89
4 changed files with 145 additions and 137 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

125
src/compute_coord_atom.cpp
View File

@ -36,16 +36,15 @@ using namespace LAMMPS_NS;
ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, 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"); if (narg < 5) error->all(FLERR,"Illegal compute coord/atom command");
int n = strlen(arg[3]) + 1; cstyle = NONE;
cstyle = new char[n];
strcpy(cstyle,arg[3]);
if (strcmp(cstyle,"cutoff") == 0) {
if (strcmp(arg[3],"cutoff") == 0) {
cstyle = CUTOFF;
double cutoff = force->numeric(FLERR,arg[4]); double cutoff = force->numeric(FLERR,arg[4]);
cutsq = cutoff*cutoff; cutsq = cutoff*cutoff;
@ -68,34 +67,33 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
ncol++; ncol++;
iarg++; 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; int iorientorder = modify->find_compute(id_orientorder);
id_orientorder = new char[n]; if (iorientorder < 0)
strcpy(id_orientorder,arg[4]); 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); threshold = force->numeric(FLERR,arg[5]);
if (iorientorder < 0) if (threshold <= -1.0 || threshold >= 1.0)
error->all(FLERR,"Could not find compute coord/atom compute ID"); error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1");
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]); ncol = 1;
if (threshold <= -1.0 || threshold >= 1.0) typelo = new int[ncol];
error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1"); typehi = new int[ncol];
typelo[0] = 1;
typehi[0] = atom->ntypes;
ncol = 1; } else error->all(FLERR,"Invalid cstyle in compute coord/atom");
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");
peratom_flag = 1; peratom_flag = 1;
if (ncol == 1) size_peratom_cols = 0; if (ncol == 1) size_peratom_cols = 0;
@ -112,21 +110,23 @@ ComputeCoordAtom::~ComputeCoordAtom()
delete [] typehi; delete [] typehi;
memory->destroy(cvec); memory->destroy(cvec);
memory->destroy(carray); memory->destroy(carray);
delete [] id_orientorder;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputeCoordAtom::init() void ComputeCoordAtom::init()
{ {
if (strcmp(cstyle,"orientorder") == 0) { if (cstyle == ORIENT) {
int iorientorder = modify->find_compute(id_orientorder); int iorientorder = modify->find_compute(id_orientorder);
c_orientorder = (ComputeOrientOrderAtom*)(modify->compute[iorientorder]); c_orientorder = (ComputeOrientOrderAtom*)(modify->compute[iorientorder]);
cutsq = c_orientorder->cutsq; cutsq = c_orientorder->cutsq;
l = c_orientorder->qlcomp; 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); comm_forward = 2*(2*l+1);
if (c_orientorder->iqlcomp < 0) 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) 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)) { if (!(c_orientorder->invoked_flag & INVOKED_PERATOM)) {
c_orientorder->compute_peratom(); c_orientorder->compute_peratom();
c_orientorder->invoked_flag |= INVOKED_PERATOM; c_orientorder->invoked_flag |= INVOKED_PERATOM;
@ -217,7 +217,7 @@ void ComputeCoordAtom::compute_peratom()
int *type = atom->type; int *type = atom->type;
int *mask = atom->mask; int *mask = atom->mask;
if (strcmp(cstyle,"cutoff") == 0) { if (cstyle == CUTOFF) {
if (ncol == 1) { if (ncol == 1) {
@ -241,7 +241,7 @@ void ComputeCoordAtom::compute_peratom()
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0]) if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0])
n++; n++;
} }
cvec[i] = n; cvec[i] = n;
@ -281,37 +281,36 @@ void ComputeCoordAtom::compute_peratom()
} }
} }
} else if (strcmp(cstyle,"orientorder") == 0) { } else if (cstyle == ORIENT) {
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
i = ilist[ii]; i = ilist[ii];
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xtmp = x[i][0]; xtmp = x[i][0];
ytmp = x[i][1]; ytmp = x[i][1];
ztmp = x[i][2]; ztmp = x[i][2];
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
n = 0; n = 0;
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; j = jlist[jj];
j &= NEIGHMASK; j &= NEIGHMASK;
delx = xtmp - x[j][0];
delx = xtmp - x[j][0]; dely = ytmp - x[j][1];
dely = ytmp - x[j][1]; delz = ztmp - x[j][2];
delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz;
rsq = delx*delx + dely*dely + delz*delz; if (rsq < cutsq) {
if (rsq < cutsq) { double dot_product = 0.0;
double dot_product = 0.0; for (int m=0; m < 2*(2*l+1); m++) {
for (int m=0; m < 2*(2*l+1); m++) { dot_product += normv[i][nqlist+m]*normv[j][nqlist+m];
dot_product += normv[i][nqlist+m]*normv[j][nqlist+m]; }
} if (dot_product > threshold) n++;
if (dot_product > threshold) n++; }
} }
} cvec[i] = n;
cvec[i] = n; } else cvec[i] = 0.0;
} else cvec[i] = 0.0; }
}
} }
} }

5
src/compute_coord_atom.h
View File

@ -34,6 +34,7 @@ class ComputeCoordAtom : public Compute {
int pack_forward_comm(int, int *, double *, int, int *); int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *); void unpack_forward_comm(int, int, double *);
double memory_usage(); double memory_usage();
enum {NONE,CUTOFF,ORIENT};
private: private:
int nmax,ncol; int nmax,ncol;
@ -45,10 +46,10 @@ class ComputeCoordAtom : public Compute {
double **carray; double **carray;
class ComputeOrientOrderAtom *c_orientorder; class ComputeOrientOrderAtom *c_orientorder;
char *cstyle,*id_orientorder; char *id_orientorder;
double threshold; double threshold;
double **normv; double **normv;
int nqlist,l; int cstyle,nqlist,l;
}; };
} }

150
src/compute_orientorder_atom.cpp
View File

@ -46,7 +46,8 @@ using namespace std;
ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) : ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, 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"); 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; qlcompflag = 0;
// specify which orders to request // specify which orders to request
nqlist = 5; nqlist = 5;
memory->create(qlist,nqlist,"orientorder/atom:qlist"); memory->create(qlist,nqlist,"orientorder/atom:qlist");
qlist[0] = 4; qlist[0] = 4;
@ -73,48 +74,55 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"nnn") == 0) { if (strcmp(arg[iarg],"nnn") == 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");
if (strcmp(arg[iarg+1],"NULL") == 0) if (strcmp(arg[iarg+1],"NULL") == 0) {
nnn = 0; nnn = 0;
else { } else {
nnn = force->numeric(FLERR,arg[iarg+1]); nnn = force->numeric(FLERR,arg[iarg+1]);
if (nnn <= 0) if (nnn <= 0)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
} }
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"degrees") == 0) { } 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]); 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->destroy(qlist);
memory->create(qlist,nqlist,"orientorder/atom:qlist"); memory->create(qlist,nqlist,"orientorder/atom:qlist");
iarg += 2; iarg += 2;
if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
qmax = 0; qmax = 0;
for (int iw = 0; iw < nqlist; iw++) { for (int iw = 0; iw < nqlist; iw++) {
qlist[iw] = force->numeric(FLERR,arg[iarg+iw]); qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
if (qlist[iw] < 0) if (qlist[iw] < 0)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
if (qlist[iw] > qmax) qmax = qlist[iw]; if (qlist[iw] > qmax) qmax = qlist[iw];
} }
iarg += nqlist; iarg += nqlist;
if (strcmp(arg[iarg],"components") == 0) { if (strcmp(arg[iarg],"components") == 0) {
qlcompflag = 1; qlcompflag = 1;
if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); if (iarg+2 > narg)
error->all(FLERR,"Illegal compute orientorder/atom command");
qlcomp = force->numeric(FLERR,arg[iarg+1]); qlcomp = force->numeric(FLERR,arg[iarg+1]);
if (qlcomp <= 0) error->all(FLERR,"Illegal compute orientorder/atom command"); if (qlcomp <= 0)
iqlcomp = -1; error->all(FLERR,"Illegal compute orientorder/atom command");
iqlcomp = -1;
for (int iw = 0; iw < nqlist; iw++) for (int iw = 0; iw < nqlist; iw++)
if (qlcomp == qlist[iw]) { if (qlcomp == qlist[iw]) {
iqlcomp = iw; iqlcomp = iw;
break; break;
} }
if (iqlcomp < 0) error->all(FLERR,"Illegal compute orientorder/atom command"); if (iqlcomp < 0)
error->all(FLERR,"Illegal compute orientorder/atom command");
iarg += 2; iarg += 2;
} }
} else if (strcmp(arg[iarg],"cutoff") == 0) { } 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]); 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; cutsq = cutoff*cutoff;
iarg += 2; iarg += 2;
} else error->all(FLERR,"Illegal compute orientorder/atom command"); } else error->all(FLERR,"Illegal compute orientorder/atom command");
@ -141,7 +149,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
memory->destroy(qlist); memory->destroy(qlist);
memory->destroy(qnm_r); memory->destroy(qnm_r);
memory->destroy(qnm_i); memory->destroy(qnm_i);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -224,7 +232,7 @@ void ComputeOrientOrderAtom::compute_peratom()
ztmp = x[i][2]; ztmp = x[i][2];
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
// insure distsq and nearest arrays are long enough // insure distsq and nearest arrays are long enough
if (jnum > maxneigh) { if (jnum > maxneigh) {
@ -253,9 +261,9 @@ void ComputeOrientOrderAtom::compute_peratom()
rsq = delx*delx + dely*dely + delz*delz; rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq) { if (rsq < cutsq) {
distsq[ncount] = rsq; distsq[ncount] = rsq;
rlist[ncount][0] = delx; rlist[ncount][0] = delx;
rlist[ncount][1] = dely; rlist[ncount][1] = dely;
rlist[ncount][2] = delz; rlist[ncount][2] = delz;
nearest[ncount++] = j; nearest[ncount++] = j;
} }
} }
@ -263,16 +271,16 @@ void ComputeOrientOrderAtom::compute_peratom()
// if not nnn neighbors, order parameter = 0; // if not nnn neighbors, order parameter = 0;
if ((ncount == 0) || (ncount < nnn)) { if ((ncount == 0) || (ncount < nnn)) {
for (int iw = 0; iw < nqlist; iw++) for (int iw = 0; iw < nqlist; iw++)
qn[iw] = 0.0; qn[iw] = 0.0;
continue; continue;
} }
// if nnn > 0, use only nearest nnn neighbors // if nnn > 0, use only nearest nnn neighbors
if (nnn > 0) { if (nnn > 0) {
select3(nnn,ncount,distsq,nearest,rlist); select3(nnn,ncount,distsq,nearest,rlist);
ncount = nnn; ncount = nnn;
} }
calc_boop(rlist, ncount, qn, qlist, nqlist); calc_boop(rlist, ncount, qn, qlist, nqlist);
@ -287,8 +295,8 @@ void ComputeOrientOrderAtom::compute_peratom()
double ComputeOrientOrderAtom::memory_usage() double ComputeOrientOrderAtom::memory_usage()
{ {
double bytes = ncol*nmax * sizeof(double); double bytes = ncol*nmax * sizeof(double);
bytes += (qmax*(2*qmax+1)+maxneigh*4) * sizeof(double); bytes += (qmax*(2*qmax+1)+maxneigh*4) * sizeof(double);
bytes += (nqlist+maxneigh) * sizeof(int); bytes += (nqlist+maxneigh) * sizeof(int);
return bytes; return bytes;
} }
@ -300,18 +308,18 @@ double ComputeOrientOrderAtom::memory_usage()
// Use no-op do while to create single statement // Use no-op do while to create single statement
#define SWAP(a,b) do { \ #define SWAP(a,b) do { \
tmp = a; a = b; b = tmp; \ tmp = a; a = b; b = tmp; \
} while(0) } while(0)
#define ISWAP(a,b) do { \ #define ISWAP(a,b) do { \
itmp = a; a = b; b = itmp; \ itmp = a; a = b; b = itmp; \
} while(0) } while(0)
#define SWAP3(a,b) do { \ #define SWAP3(a,b) do { \
tmp = a[0]; a[0] = b[0]; b[0] = tmp; \ tmp = a[0]; a[0] = b[0]; b[0] = tmp; \
tmp = a[1]; a[1] = b[1]; b[1] = tmp; \ tmp = a[1]; a[1] = b[1]; b[1] = tmp; \
tmp = a[2]; a[2] = b[2]; b[2] = tmp; \ tmp = a[2]; a[2] = b[2]; b[2] = tmp; \
} while(0) } 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) {
if (ir == l+1 && arr[ir] < arr[l]) { if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir]); SWAP(arr[l],arr[ir]);
ISWAP(iarr[l],iarr[ir]); ISWAP(iarr[l],iarr[ir]);
SWAP3(arr3[l],arr3[ir]); SWAP3(arr3[l],arr3[ir]);
} }
return; return;
@ -342,17 +350,17 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (arr[l] > arr[ir]) { if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir]); SWAP(arr[l],arr[ir]);
ISWAP(iarr[l],iarr[ir]); ISWAP(iarr[l],iarr[ir]);
SWAP3(arr3[l],arr3[ir]); SWAP3(arr3[l],arr3[ir]);
} }
if (arr[l+1] > arr[ir]) { if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir]); SWAP(arr[l+1],arr[ir]);
ISWAP(iarr[l+1],iarr[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]) { if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1]); SWAP(arr[l],arr[l+1]);
ISWAP(iarr[l],iarr[l+1]); ISWAP(iarr[l],iarr[l+1]);
SWAP3(arr3[l],arr3[l+1]); SWAP3(arr3[l],arr3[l+1]);
} }
i = l+1; i = l+1;
j = ir; j = ir;
@ -367,7 +375,7 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (j < i) break; if (j < i) break;
SWAP(arr[i],arr[j]); SWAP(arr[i],arr[j]);
ISWAP(iarr[i],iarr[j]); ISWAP(iarr[i],iarr[j]);
SWAP3(arr3[i],arr3[j]); SWAP3(arr3[i],arr3[j]);
} }
arr[l+1] = arr[j]; arr[l+1] = arr[j];
arr[j] = a; 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 calculate the bond orientational order parameters
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void ComputeOrientOrderAtom::calc_boop(double **rlist, void ComputeOrientOrderAtom::calc_boop(double **rlist,
int ncount, double qn[], int ncount, double qn[],
int qlist[], int nqlist) { int qlist[], int nqlist) {
for (int iw = 0; iw < nqlist; iw++) { for (int iw = 0; iw < nqlist; iw++) {
int n = qlist[iw]; int n = qlist[iw];
@ -429,22 +437,22 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
double expphim_r = expphi_r; double expphim_r = expphi_r;
double expphim_i = expphi_i; double expphim_i = expphi_i;
for(int m = 1; m <= +n; m++) { for(int m = 1; m <= +n; m++) {
double prefactor = polar_prefactor(n, m, costheta); double prefactor = polar_prefactor(n, m, costheta);
double c_r = prefactor * expphim_r; double c_r = prefactor * expphim_r;
double c_i = prefactor * expphim_i; double c_i = prefactor * expphim_i;
qnm_r[iw][m+n] += c_r; qnm_r[iw][m+n] += c_r;
qnm_i[iw][m+n] += c_i; qnm_i[iw][m+n] += c_i;
if(m & 1) { if(m & 1) {
qnm_r[iw][-m+n] -= c_r; qnm_r[iw][-m+n] -= c_r;
qnm_i[iw][-m+n] += c_i; qnm_i[iw][-m+n] += c_i;
} else { } else {
qnm_r[iw][-m+n] += c_r; qnm_r[iw][-m+n] += c_r;
qnm_i[iw][-m+n] -= c_i; qnm_i[iw][-m+n] -= c_i;
} }
double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i; double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r; double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
expphim_r = tmp_r; expphim_r = tmp_r;
expphim_i = tmp_i; expphim_i = tmp_i;
} }
} }
@ -458,15 +466,15 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
for(int m = 0; m < 2*n+1; m++) { 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]; 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, // 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)); qn[iw] = fac * sqrt(qm_sum / (2*n+1));
if (qlcompflag && iqlcomp == iw) normfac = 1.0/sqrt(qm_sum); if (qlcompflag && iqlcomp == iw) normfac = 1.0/sqrt(qm_sum);
} }
// output of the complex vector // output of the complex vector
if (qlcompflag) { if (qlcompflag) {
int j = nqlist; int j = nqlist;
for(int m = 0; m < 2*qlcomp+1; m++) { 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) Y_l^m (theta, phi) = prefactor(l, m, cos(theta)) * exp(i*m*phi)
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

2
src/compute_orientorder_atom.h
View File

@ -49,7 +49,7 @@ class ComputeOrientOrderAtom : public Compute {
double **qnm_i; double **qnm_i;
void select3(int, int, double *, int *, double **); 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 qn[], int nlist[], int nnlist);
double dist(const double r[]); double dist(const double r[]);