From 12c8aaf29d5cac0ee92c4cc9b5bce76db5385554 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 11 Jan 2016 23:38:16 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14423 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_chunk_atom.cpp | 91 ++++++++++++++++++- src/compute_chunk_atom.h | 2 +- src/fix_store.cpp | 23 ++--- src/fix_store.h | 3 +- src/fix_store_state.cpp | 175 +++++++++++++++++++++++++++++++++++++ src/fix_store_state.h | 9 ++ 6 files changed, 288 insertions(+), 15 deletions(-) diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp index 9717604111..d4e15ba018 100644 --- a/src/compute_chunk_atom.cpp +++ b/src/compute_chunk_atom.cpp @@ -191,6 +191,7 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) : maxflag[1] = UPPER; maxflag[2] = UPPER; scaleflag = LATTICE; + pbcflag = 0; while (iarg < narg) { if (strcmp(arg[iarg],"region") == 0) { @@ -269,6 +270,12 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg+1],"reduced") == 0) scaleflag = REDUCED; else error->all(FLERR,"Illegal compute chunk/atom command"); iarg += 2; + } else if (strcmp(arg[iarg],"pbc") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command"); + if (strcmp(arg[iarg+1],"no") == 0) pbcflag = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) pbcflag = 1; + else error->all(FLERR,"Illegal compute chunk/atom command"); + iarg += 2; } else error->all(FLERR,"Illegal compute chunk/atom command"); } @@ -1265,7 +1272,7 @@ int ComputeChunkAtom::setup_xyz_bins() int ComputeChunkAtom::setup_sphere_bins() { // convert sorigin_user to sorigin - // sorigin is always in box units, for orthogonal or triclinic domains + // sorigin,srad are always in box units, for orthogonal or triclinic domains // lamda2x works for either orthogonal or triclinic if (scaleflag == REDUCED) { @@ -1280,6 +1287,23 @@ int ComputeChunkAtom::setup_sphere_bins() sradmax = sradmax_user; } + // if pbcflag set, sradmax must be < 1/2 box in any periodic dim + // treat orthongonal and triclinic the same + // check every time bins are created + + if (pbcflag) { + double *prd_half = domain->prd_half; + int *periodicity = domain->periodicity; + int flag = 0; + if (periodicity[0] && sradmax > prd_half[0]) flag = 1; + if (periodicity[1] && sradmax > prd_half[1]) flag = 1; + if (domain->dimension == 3 && + periodicity[2] && sradmax > prd_half[2]) flag = 1; + if (flag) + error->all(FLERR,"Compute chunk/atom bin/sphere radius " + "is too large for periodic box"); + } + sinvrad = nsbin / (sradmax-sradmin); // allocate and set bin coordinates @@ -1328,6 +1352,21 @@ int ComputeChunkAtom::setup_cylinder_bins() cradmax = cradmax_user; } + // if pbcflag set, sradmax must be < 1/2 box in any periodic non-axis dim + // treat orthongonal and triclinic the same + // check every time bins are created + + if (pbcflag) { + double *prd_half = domain->prd_half; + int *periodicity = domain->periodicity; + int flag = 0; + if (periodicity[cdim1] && sradmax > prd_half[cdim1]) flag = 1; + if (periodicity[cdim2] && sradmax > prd_half[cdim2]) flag = 1; + if (flag) + error->all(FLERR,"Compute chunk/atom bin/cylinder radius " + "is too large for periodic box"); + } + cinvrad = ncbin / (cradmax-cradmin); // allocate and set radial bin coordinates @@ -1743,6 +1782,7 @@ void ComputeChunkAtom::atom2binsphere() double *boxlo = domain->boxlo; double *boxhi = domain->boxhi; double *prd = domain->prd; + double *prd_half = domain->prd_half; int *periodicity = domain->periodicity; // remap each atom's relevant coords back into box via PBC if necessary @@ -1773,6 +1813,33 @@ void ComputeChunkAtom::atom2binsphere() dx = xremap - sorigin[0]; dy = yremap - sorigin[1]; dz = zremap - sorigin[2]; + + // if requested, apply PBC to distance from sphere center + // treat orthogonal and triclinic the same + // with dx,dy,dz = lengths independent of each other + // so do not use domain->minimum_image() which couples for triclinic + + if (pbcflag) { + if (periodicity[0]) { + if (fabs(dx) > prd_half[0]) { + if (dx < 0.0) dx += prd[0]; + else dx -= prd[0]; + } + } + if (periodicity[1]) { + if (fabs(dy) > prd_half[1]) { + if (dy < 0.0) dy += prd[1]; + else dy -= prd[1]; + } + } + if (periodicity[2]) { + if (fabs(dz) > prd_half[2]) { + if (dz < 0.0) dz += prd[2]; + else dz -= prd[2]; + } + } + } + r = sqrt(dx*dx + dy*dy + dz*dz); ibin = static_cast ((r - sradmin) * sinvrad); @@ -1792,7 +1859,6 @@ void ComputeChunkAtom::atom2binsphere() /* ---------------------------------------------------------------------- assign each atom to a cylindrical bin - ------------------------------------------------------------------------- */ void ComputeChunkAtom::atom2bincylinder() @@ -1812,6 +1878,7 @@ void ComputeChunkAtom::atom2bincylinder() double *boxlo = domain->boxlo; double *boxhi = domain->boxhi; double *prd = domain->prd; + double *prd_half = domain->prd_half; int *periodicity = domain->periodicity; // remap each atom's relevant coords back into box via PBC if necessary @@ -1837,6 +1904,26 @@ void ComputeChunkAtom::atom2bincylinder() d1 = remap1 - corigin[cdim1]; d2 = remap2 - corigin[cdim2]; + + // if requested, apply PBC to distance from cylinder axis + // treat orthogonal and triclinic the same + // with d1,d2 = lengths independent of each other + + if (pbcflag) { + if (periodicity[cdim1]) { + if (fabs(d1) > prd_half[cdim1]) { + if (d1 < 0.0) d1 += prd[cdim1]; + else d1 -= prd[cdim1]; + } + } + if (periodicity[cdim2]) { + if (fabs(d2) > prd_half[cdim2]) { + if (d2 < 0.0) d2 += prd[cdim2]; + else d2 -= prd[cdim2]; + } + } + } + r = sqrt(d1*d1 + d2*d2); rbin = static_cast ((r - cradmin) * cinvrad); diff --git a/src/compute_chunk_atom.h b/src/compute_chunk_atom.h index c1609ff6bc..6330a37ad8 100644 --- a/src/compute_chunk_atom.h +++ b/src/compute_chunk_atom.h @@ -51,7 +51,7 @@ class ComputeChunkAtom : public Compute { int which,binflag; int regionflag,nchunksetflag,nchunkflag,discard; int limit,limitstyle,limitfirst; - int scaleflag; + int scaleflag,pbcflag; double xscale,yscale,zscale; int argindex; char *cfvid; diff --git a/src/fix_store.cpp b/src/fix_store.cpp index 4e3d612cce..1375384f6c 100644 --- a/src/fix_store.cpp +++ b/src/fix_store.cpp @@ -15,9 +15,9 @@ #include #include "fix_store.h" #include "atom.h" +#include "force.h" #include "memory.h" #include "error.h" -#include "force.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -80,16 +80,6 @@ int FixStore::setmask() return mask; } -/* ---------------------------------------------------------------------- - memory usage of local atom-based array -------------------------------------------------------------------------- */ - -double FixStore::memory_usage() -{ - double bytes = atom->nmax*nvalues * sizeof(double); - return bytes; -} - /* ---------------------------------------------------------------------- allocate atom-based array ------------------------------------------------------------------------- */ @@ -189,3 +179,14 @@ int FixStore::size_restart(int nlocal) { return nvalues+1; } + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double FixStore::memory_usage() +{ + double bytes = atom->nmax*nvalues * sizeof(double); + return bytes; +} + diff --git a/src/fix_store.h b/src/fix_store.h index dab9899906..64b714c3da 100644 --- a/src/fix_store.h +++ b/src/fix_store.h @@ -33,7 +33,6 @@ class FixStore : public Fix { ~FixStore(); int setmask(); - double memory_usage(); void grow_arrays(int); void copy_arrays(int, int, int); int pack_exchange(int, double *); @@ -43,6 +42,8 @@ class FixStore : public Fix { int size_restart(int); int maxsize_restart(); + double memory_usage(); + private: int nvalues; // total # of values per atom int vecflag; // 1 if nvalues = 1 diff --git a/src/fix_store_state.cpp b/src/fix_store_state.cpp index efcc73f032..aec61e6980 100644 --- a/src/fix_store_state.cpp +++ b/src/fix_store_state.cpp @@ -105,6 +105,19 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) : if (domain->triclinic) pack_choice[nvalues++] = &FixStoreState::pack_zu_triclinic; else pack_choice[nvalues++] = &FixStoreState::pack_zu; + } else if (strcmp(arg[iarg],"xsu") == 0) { + if (domain->triclinic) + pack_choice[nvalues++] = &FixStoreState::pack_xsu_triclinic; + else pack_choice[nvalues++] = &FixStoreState::pack_xsu; + } else if (strcmp(arg[iarg],"ysu") == 0) { + if (domain->triclinic) + pack_choice[nvalues++] = &FixStoreState::pack_ysu_triclinic; + else pack_choice[nvalues++] = &FixStoreState::pack_ysu; + } else if (strcmp(arg[iarg],"zsu") == 0) { + if (domain->triclinic) + pack_choice[nvalues++] = &FixStoreState::pack_zsu_triclinic; + else pack_choice[nvalues++] = &FixStoreState::pack_zsu; + } else if (strcmp(arg[iarg],"ix") == 0) { pack_choice[nvalues++] = &FixStoreState::pack_ix; } else if (strcmp(arg[iarg],"iy") == 0) { @@ -145,12 +158,22 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Fix store/state for atom property that isn't allocated"); pack_choice[nvalues++] = &FixStoreState::pack_muz; + } else if (strcmp(arg[iarg],"mu") == 0) { + if (!atom->mu_flag) + error->all(FLERR, + "Fix store/state for atom property that isn't allocated"); + pack_choice[nvalues++] = &FixStoreState::pack_mu; } else if (strcmp(arg[iarg],"radius") == 0) { if (!atom->radius_flag) error->all(FLERR, "Fix store/state for atom property that isn't allocated"); pack_choice[nvalues++] = &FixStoreState::pack_radius; + } else if (strcmp(arg[iarg],"diameter") == 0) { + if (!atom->radius_flag) + error->all(FLERR, + "Fix store/state for atom property that isn't allocated"); + pack_choice[nvalues++] = &FixStoreState::pack_diameter; } else if (strcmp(arg[iarg],"omegax") == 0) { if (!atom->omega_flag) error->all(FLERR, @@ -989,6 +1012,128 @@ void FixStoreState::pack_zu_triclinic(int n) /* ---------------------------------------------------------------------- */ +void FixStoreState::pack_xsu(int n) +{ + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double boxxlo = domain->boxlo[0]; + double invxprd = 1.0/domain->xprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + vbuf[n] = (x[i][0]-boxxlo)*invxprd + ((image[i] & IMGMASK) - IMGMAX); + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixStoreState::pack_ysu(int n) +{ + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double boxylo = domain->boxlo[1]; + double invyprd = 1.0/domain->yprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + vbuf[n] = (x[i][1]-boxylo)*invyprd + (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixStoreState::pack_zsu(int n) +{ + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double boxzlo = domain->boxlo[2]; + double invzprd = 1.0/domain->zprd; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + vbuf[n] = (x[i][2]-boxzlo)*invzprd + (image[i] >> IMG2BITS) - IMGMAX; + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_xsu_triclinic(int n) +{ + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + vbuf[n] = h_inv[0]*(x[i][0]-boxlo[0]) + h_inv[5]*(x[i][1]-boxlo[1]) + + h_inv[4]*(x[i][2]-boxlo[2]) + (image[i] & IMGMASK) - IMGMAX; + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_ysu_triclinic(int n) +{ + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + vbuf[n] = h_inv[1]*(x[i][1]-boxlo[1]) + h_inv[3]*(x[i][2]-boxlo[2]) + + (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_zsu_triclinic(int n) +{ + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = h_inv[2]*(x[i][2]-boxlo[2]) + (image[i] >> IMG2BITS) - IMGMAX; + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void FixStoreState::pack_ix(int n) { imageint *image = atom->image; @@ -1184,6 +1329,21 @@ void FixStoreState::pack_muz(int n) /* ---------------------------------------------------------------------- */ +void FixStoreState::pack_mu(int n) +{ + double **mu = atom->mu; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) vbuf[n] = mu[i][3]; + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void FixStoreState::pack_radius(int n) { double *radius = atom->radius; @@ -1199,6 +1359,21 @@ void FixStoreState::pack_radius(int n) /* ---------------------------------------------------------------------- */ +void FixStoreState::pack_diameter(int n) +{ + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) vbuf[n] = 2.0*radius[i]; + else vbuf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void FixStoreState::pack_omegax(int n) { double **omega = atom->omega; diff --git a/src/fix_store_state.h b/src/fix_store_state.h index d7523d46bd..6397f16ef7 100644 --- a/src/fix_store_state.h +++ b/src/fix_store_state.h @@ -79,6 +79,13 @@ class FixStoreState : public Fix { void pack_xu_triclinic(int); void pack_yu_triclinic(int); void pack_zu_triclinic(int); + void pack_xsu(int); + void pack_ysu(int); + void pack_zsu(int); + void pack_xsu_triclinic(int); + void pack_ysu_triclinic(int); + void pack_zsu_triclinic(int); + void pack_ix(int); void pack_iy(int); void pack_iz(int); @@ -93,7 +100,9 @@ class FixStoreState : public Fix { void pack_mux(int); void pack_muy(int); void pack_muz(int); + void pack_mu(int); void pack_radius(int); + void pack_diameter(int); void pack_omegax(int); void pack_omegay(int); void pack_omegaz(int);