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) { + } 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; - } + 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[]);