From 9e1f8bcd6ad42714046afcdab318edf57bf06e93 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 8 Mar 2007 00:54:02 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@371 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/atom.cpp | 72 +++-- src/atom_vec.h | 8 +- src/atom_vec_atomic.cpp | 50 ++-- src/atom_vec_atomic.h | 8 +- src/atom_vec_charge.cpp | 50 ++-- src/atom_vec_charge.h | 8 +- src/atom_vec_hybrid.cpp | 43 ++- src/atom_vec_hybrid.h | 8 +- src/comm.cpp | 200 ++++++++++---- src/comm.h | 10 +- src/compute.h | 2 +- src/create_atoms.cpp | 158 ++++++----- src/create_atoms.h | 7 - src/create_box.cpp | 47 ++-- src/delete_atoms.cpp | 4 +- src/delete_bonds.cpp | 2 + src/displace_atoms.cpp | 2 + src/domain.cpp | 522 +++++++++++++++++++++++++++++-------- src/domain.h | 94 ++++--- src/dump.cpp | 19 ++ src/dump.h | 4 + src/dump_atom.cpp | 27 +- src/dump_custom.cpp | 32 ++- src/dump_dcd.cpp | 23 +- src/fix.cpp | 1 + src/fix.h | 3 +- src/fix_ave_spatial.cpp | 9 + src/fix_deposit.cpp | 73 +++--- src/fix_drag.cpp | 2 +- src/fix_nph.cpp | 10 +- src/fix_npt.cpp | 10 +- src/fix_orient_fcc.cpp | 3 +- src/fix_orient_fcc.h | 2 +- src/fix_recenter.cpp | 21 +- src/fix_shake.cpp | 19 +- src/fix_shake.h | 2 +- src/fix_uniaxial.cpp | 6 +- src/fix_volume_rescale.cpp | 6 + src/fix_wall_lj126.cpp | 8 + src/fix_wall_lj93.cpp | 8 + src/fix_wall_reflect.cpp | 7 + src/lattice.cpp | 76 ++++-- src/neigh_full.cpp | 15 +- src/neigh_gran.cpp | 127 ++++++--- src/neigh_half.cpp | 131 +++++++--- src/neigh_respa.cpp | 177 ++++++++++--- src/neighbor.cpp | 127 ++++++--- src/neighbor.h | 9 + src/pair.h | 2 +- src/read_data.cpp | 11 +- src/read_restart.cpp | 44 +++- src/region_block.cpp | 18 +- src/region_cylinder.cpp | 30 ++- src/region_prism.cpp | 43 +-- src/region_prism.h | 4 +- src/replicate.cpp | 95 +++++-- src/respa.cpp | 8 + src/respa.h | 1 + src/set.cpp | 2 + src/velocity.cpp | 28 +- src/velocity.h | 2 +- src/verlet.cpp | 8 + src/verlet.h | 1 + src/write_restart.cpp | 8 + 64 files changed, 1785 insertions(+), 772 deletions(-) diff --git a/src/atom.cpp b/src/atom.cpp index 0a96553168..7775f5cda3 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -35,18 +35,17 @@ using namespace LAMMPS_NS; -#define MIN(A,B) ((A) < (B)) ? (A) : (B) -#define MAX(A,B) ((A) > (B)) ? (A) : (B) - #define DELTA 1 #define DELTA_MEMSTR 1024 +#define EPSILON 1.0e-6 + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) /* ---------------------------------------------------------------------- */ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) { - // no atoms - natoms = nlocal = nghost = nmax = 0; ntypes = 0; nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; @@ -620,8 +619,9 @@ int Atom::count_words(char *line) void Atom::data_atoms(int n, char *buf, int ihybrid) { - int m,imagetmp,xptr,iptr; - double xtmp,ytmp,ztmp; + int m,imagedata,xptr,iptr; + double xdata[3],lamda[3],sublo[3],subhi[3]; + double *coord; char *next; next = strchr(buf,'\n'); @@ -635,12 +635,33 @@ void Atom::data_atoms(int n, char *buf, int ihybrid) char **values = new char*[nwords]; - double subxlo = domain->subxlo; - double subxhi = domain->subxhi; - double subylo = domain->subylo; - double subyhi = domain->subyhi; - double subzlo = domain->subzlo; - double subzhi = domain->subzhi; + // set bounds for my proc + // if periodic and I am lo/hi proc, adjust bounds by EPSILON + // insures all data atoms will be owned even with round-off + + int triclinic = domain->triclinic; + if (triclinic == 0) { + sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0]; + sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1]; + sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2]; + } else { + sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0]; + sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1]; + sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2]; + } + + if (domain->xperiodic) { + if (comm->myloc[0] == 0) sublo[0] -= EPSILON; + if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += EPSILON; + } + if (domain->yperiodic) { + if (comm->myloc[1] == 0) sublo[1] -= EPSILON; + if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += EPSILON; + } + if (domain->zperiodic) { + if (comm->myloc[2] == 0) sublo[2] -= EPSILON; + if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += EPSILON; + } // xptr = which word in line starts xyz coords // iptr = which word in line starts ix,iy,iz image flags @@ -662,20 +683,25 @@ void Atom::data_atoms(int n, char *buf, int ihybrid) for (m = 1; m < nwords; m++) values[m] = strtok(NULL," \t\n\r\f"); - xtmp = atof(values[xptr]); - ytmp = atof(values[xptr+1]); - ztmp = atof(values[xptr+2]); if (imageflag) - imagetmp = ((atoi(values[iptr+2]) + 512 & 1023) << 20) | + imagedata = ((atoi(values[iptr+2]) + 512 & 1023) << 20) | ((atoi(values[iptr+1]) + 512 & 1023) << 10) | (atoi(values[iptr]) + 512 & 1023); - else imagetmp = (512 << 20) | (512 << 10) | 512; + else imagedata = (512 << 20) | (512 << 10) | 512; - domain->remap(xtmp,ytmp,ztmp,imagetmp); - if (xtmp >= subxlo && xtmp < subxhi && - ytmp >= subylo && ytmp < subyhi && - ztmp >= subzlo && ztmp < subzhi) - avec->data_atom(xtmp,ytmp,ztmp,imagetmp,values,ihybrid); + xdata[0] = atof(values[xptr]); + xdata[1] = atof(values[xptr+1]); + xdata[2] = atof(values[xptr+2]); + domain->remap(xdata,imagedata); + if (triclinic) { + domain->x2lamda(xdata,lamda); + coord = lamda; + } else coord = xdata; + + if (coord[0] >= sublo[0] && coord[0] < subhi[0] && + coord[1] >= sublo[1] && coord[1] < subhi[1] && + coord[2] >= sublo[2] && coord[2] < subhi[2]) + avec->data_atom(xdata,imagedata,values,ihybrid); buf = next + 1; } diff --git a/src/atom_vec.h b/src/atom_vec.h index e23950c888..8278fdea5a 100644 --- a/src/atom_vec.h +++ b/src/atom_vec.h @@ -44,7 +44,7 @@ class AtomVec : protected Pointers { virtual void zero_ghost(int,int) {} virtual void copy(int, int) = 0; - virtual int pack_comm(int, int *, double *, int *) = 0; + virtual int pack_comm(int, int *, double *, int, double *) = 0; virtual int pack_comm_one(int, double *) {return 0;} virtual void unpack_comm(int, int, double *) = 0; virtual int unpack_comm_one(int, double *) {return 0;} @@ -54,7 +54,7 @@ class AtomVec : protected Pointers { virtual void unpack_reverse(int, int *, double *) = 0; virtual int unpack_reverse_one(int, double *) {return 0;} - virtual int pack_border(int, int *, double *, int *) = 0; + virtual int pack_border(int, int *, double *, int, double *) = 0; virtual int pack_border_one(int, double *) {return 0;} virtual void unpack_border(int, int, double *) = 0; virtual int unpack_border_one(int, double *) {return 0;} @@ -67,8 +67,8 @@ class AtomVec : protected Pointers { virtual int pack_restart(int, double *) = 0; virtual int unpack_restart(double *) = 0; - virtual void create_atom(int, double, double, double, int) = 0; - virtual void data_atom(double, double, double, int, char **, int) = 0; + virtual void create_atom(int, double *, int) = 0; + virtual void data_atom(double *, int, char **, int) = 0; virtual void data_vel(int, char *, int); virtual void data_params(int) {} diff --git a/src/atom_vec_atomic.cpp b/src/atom_vec_atomic.cpp index 4c8fa7333a..9921b30c70 100644 --- a/src/atom_vec_atomic.cpp +++ b/src/atom_vec_atomic.cpp @@ -14,7 +14,6 @@ #include "stdlib.h" #include "atom_vec_atomic.h" #include "atom.h" -#include "domain.h" #include "modify.h" #include "fix.h" #include "memory.h" @@ -102,12 +101,13 @@ void AtomVecAtomic::copy(int i, int j) /* ---------------------------------------------------------------------- */ -int AtomVecAtomic::pack_comm(int n, int *list, double *buf, int *pbc_flags) +int AtomVecAtomic::pack_comm(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -115,14 +115,11 @@ int AtomVecAtomic::pack_comm(int n, int *list, double *buf, int *pbc_flags) buf[m++] = x[j][2]; } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; } } return m; @@ -176,12 +173,13 @@ void AtomVecAtomic::unpack_reverse(int n, int *list, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecAtomic::pack_border(int n, int *list, double *buf, int *pbc_flags) +int AtomVecAtomic::pack_border(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -192,14 +190,11 @@ int AtomVecAtomic::pack_border(int n, int *list, double *buf, int *pbc_flags) buf[m++] = mask[j]; } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; @@ -378,21 +373,20 @@ int AtomVecAtomic::unpack_restart(double *buf) } /* ---------------------------------------------------------------------- - create one atom of type itype at x0,y0,z0 + create one atom of itype at coord set other values to defaults ------------------------------------------------------------------------- */ -void AtomVecAtomic::create_atom(int itype, double x0, double y0, double z0, - int ihybrid) +void AtomVecAtomic::create_atom(int itype, double *coord, int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); tag[nlocal] = 0; type[nlocal] = itype; - x[nlocal][0] = x0; - x[nlocal][1] = y0; - x[nlocal][2] = z0; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; mask[nlocal] = 1; image[nlocal] = (512 << 20) | (512 << 10) | 512; v[nlocal][0] = 0.0; @@ -407,8 +401,8 @@ void AtomVecAtomic::create_atom(int itype, double x0, double y0, double z0, initialize other atom quantities ------------------------------------------------------------------------- */ -void AtomVecAtomic::data_atom(double xtmp, double ytmp, double ztmp, - int imagetmp, char **values, int ihybrid) +void AtomVecAtomic::data_atom(double *coord, int imagetmp, char **values, + int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -421,9 +415,9 @@ void AtomVecAtomic::data_atom(double xtmp, double ytmp, double ztmp, if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) error->one("Invalid atom type in Atoms section of data file"); - x[nlocal][0] = xtmp; - x[nlocal][1] = ytmp; - x[nlocal][2] = ztmp; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; image[nlocal] = imagetmp; diff --git a/src/atom_vec_atomic.h b/src/atom_vec_atomic.h index ba98432b51..0cff9b0d0e 100644 --- a/src/atom_vec_atomic.h +++ b/src/atom_vec_atomic.h @@ -25,11 +25,11 @@ class AtomVecAtomic : public AtomVec { void grow(int); void reset_ptrs(); void copy(int, int); - virtual int pack_comm(int, int *, double *, int *); + virtual int pack_comm(int, int *, double *, int, double *); virtual void unpack_comm(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); - virtual int pack_border(int, int *, double *, int *); + virtual int pack_border(int, int *, double *, int, double *); virtual void unpack_border(int, int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); @@ -37,8 +37,8 @@ class AtomVecAtomic : public AtomVec { int size_restart_one(int); int pack_restart(int, double *); int unpack_restart(double *); - void create_atom(int, double, double, double, int); - void data_atom(double, double, double, int, char **, int); + void create_atom(int, double *, int); + void data_atom(double *, int, char **, int); int memory_usage(); protected: diff --git a/src/atom_vec_charge.cpp b/src/atom_vec_charge.cpp index 7d073afcff..20c36f2865 100644 --- a/src/atom_vec_charge.cpp +++ b/src/atom_vec_charge.cpp @@ -14,7 +14,6 @@ #include "stdlib.h" #include "atom_vec_charge.h" #include "atom.h" -#include "domain.h" #include "modify.h" #include "fix.h" #include "memory.h" @@ -134,12 +133,13 @@ void AtomVecCharge::copy(int i, int j) /* ---------------------------------------------------------------------- */ -int AtomVecCharge::pack_comm(int n, int *list, double *buf, int *pbc_flags) +int AtomVecCharge::pack_comm(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -147,14 +147,11 @@ int AtomVecCharge::pack_comm(int n, int *list, double *buf, int *pbc_flags) buf[m++] = x[j][2]; } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; } } return m; @@ -208,12 +205,13 @@ void AtomVecCharge::unpack_reverse(int n, int *list, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecCharge::pack_border(int n, int *list, double *buf, int *pbc_flags) +int AtomVecCharge::pack_border(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -225,14 +223,11 @@ int AtomVecCharge::pack_border(int n, int *list, double *buf, int *pbc_flags) buf[m++] = q[j]; } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; @@ -437,21 +432,20 @@ int AtomVecCharge::unpack_restart(double *buf) } /* ---------------------------------------------------------------------- - create one atom of type itype at x0,y0,z0 + create one atom of itype at coord set other values to defaults ------------------------------------------------------------------------- */ -void AtomVecCharge::create_atom(int itype, double x0, double y0, double z0, - int ihybrid) +void AtomVecCharge::create_atom(int itype, double *coord, int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); tag[nlocal] = 0; type[nlocal] = itype; - x[nlocal][0] = x0; - x[nlocal][1] = y0; - x[nlocal][2] = z0; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; mask[nlocal] = 1; image[nlocal] = (512 << 20) | (512 << 10) | 512; v[nlocal][0] = 0.0; @@ -468,8 +462,8 @@ void AtomVecCharge::create_atom(int itype, double x0, double y0, double z0, initialize other atom quantities ------------------------------------------------------------------------- */ -void AtomVecCharge::data_atom(double xtmp, double ytmp, double ztmp, - int imagetmp, char **values, int ihybrid) +void AtomVecCharge::data_atom(double *coord, int imagetmp, char **values, + int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -484,9 +478,9 @@ void AtomVecCharge::data_atom(double xtmp, double ytmp, double ztmp, q[nlocal] = atof(values[2]); - x[nlocal][0] = xtmp; - x[nlocal][1] = ytmp; - x[nlocal][2] = ztmp; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; image[nlocal] = imagetmp; diff --git a/src/atom_vec_charge.h b/src/atom_vec_charge.h index bca54faa67..955d74008e 100644 --- a/src/atom_vec_charge.h +++ b/src/atom_vec_charge.h @@ -26,11 +26,11 @@ class AtomVecCharge : public AtomVec { void zero_owned(int); void zero_ghost(int, int); void copy(int, int); - int pack_comm(int, int *, double *, int *); + int pack_comm(int, int *, double *, int, double *); void unpack_comm(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); - int pack_border(int, int *, double *, int *); + int pack_border(int, int *, double *, int, double *); int pack_border_one(int, double *); void unpack_border(int, int, double *); int unpack_border_one(int, double *); @@ -40,8 +40,8 @@ class AtomVecCharge : public AtomVec { int size_restart_one(int); int pack_restart(int, double *); int unpack_restart(double *); - void create_atom(int, double, double, double, int); - void data_atom(double, double, double, int, char **, int); + void create_atom(int, double *, int); + void data_atom(double *, int, char **, int); int memory_usage(); private: diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp index ff46cb824e..f5d01daf54 100644 --- a/src/atom_vec_hybrid.cpp +++ b/src/atom_vec_hybrid.cpp @@ -1,3 +1,4 @@ + /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories @@ -14,7 +15,6 @@ #include "string.h" #include "atom_vec_hybrid.h" #include "atom.h" -#include "domain.h" #include "modify.h" #include "fix.h" #include "memory.h" @@ -152,12 +152,13 @@ void AtomVecHybrid::copy(int i, int j) /* ---------------------------------------------------------------------- */ -int AtomVecHybrid::pack_comm(int n, int *list, double *buf, int *pbc_flags) +int AtomVecHybrid::pack_comm(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -166,14 +167,11 @@ int AtomVecHybrid::pack_comm(int n, int *list, double *buf, int *pbc_flags) m += styles[hybrid[j]]->pack_comm_one(j,&buf[m]); } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; m += styles[hybrid[j]]->pack_comm_one(j,&buf[m]); } } @@ -231,12 +229,13 @@ void AtomVecHybrid::unpack_reverse(int n, int *list, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecHybrid::pack_border(int n, int *list, double *buf, int *pbc_flags) +int AtomVecHybrid::pack_border(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -249,14 +248,11 @@ int AtomVecHybrid::pack_border(int n, int *list, double *buf, int *pbc_flags) m += styles[hybrid[j]]->pack_border_one(j,&buf[m]); } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; @@ -396,14 +392,13 @@ int AtomVecHybrid::unpack_restart(double *buf) } /* ---------------------------------------------------------------------- - create one atom of type itype at x0,y0,z0 for ihybrid style + create one atom of itype at coord for ihybrid style sub-style does create grow() must occur here so arrays for all sub-styles are grown zero auxiliary arrays for all other styles before create ------------------------------------------------------------------------- */ -void AtomVecHybrid::create_atom(int itype, double x0, double y0, double z0, - int ihybrid) +void AtomVecHybrid::create_atom(int itype, double *coord, int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -411,7 +406,7 @@ void AtomVecHybrid::create_atom(int itype, double x0, double y0, double z0, for (int m = 0; m < nstyles; m++) if (m != ihybrid) styles[m]->zero_owned(nlocal); hybrid[nlocal] = ihybrid; - styles[ihybrid]->create_atom(itype,x0,y0,z0,0); + styles[ihybrid]->create_atom(itype,coord,0); } /* ---------------------------------------------------------------------- @@ -420,8 +415,8 @@ void AtomVecHybrid::create_atom(int itype, double x0, double y0, double z0, sub-style will increment nlocal ------------------------------------------------------------------------- */ -void AtomVecHybrid::data_atom(double xtmp, double ytmp, double ztmp, - int imagetmp, char **values, int ihybrid) +void AtomVecHybrid::data_atom(double *coord, int imagetmp, char **values, + int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -429,7 +424,7 @@ void AtomVecHybrid::data_atom(double xtmp, double ytmp, double ztmp, for (int m = 0; m < nstyles; m++) if (m != ihybrid) styles[m]->zero_owned(nlocal); hybrid[nlocal] = ihybrid; - styles[ihybrid]->data_atom(xtmp,ytmp,ztmp,imagetmp,values,0); + styles[ihybrid]->data_atom(coord,imagetmp,values,0); } /* ---------------------------------------------------------------------- diff --git a/src/atom_vec_hybrid.h b/src/atom_vec_hybrid.h index 66b4025128..bb402d6e83 100644 --- a/src/atom_vec_hybrid.h +++ b/src/atom_vec_hybrid.h @@ -29,11 +29,11 @@ class AtomVecHybrid : public AtomVec { void grow(int); void reset_ptrs(); void copy(int, int); - int pack_comm(int, int *, double *, int *); + int pack_comm(int, int *, double *, int, double *); void unpack_comm(int, int, double *); int pack_reverse(int, int, double *); void unpack_reverse(int, int *, double *); - int pack_border(int, int *, double *, int *); + int pack_border(int, int *, double *, int, double *); void unpack_border(int, int, double *); int pack_exchange(int, double *); int unpack_exchange(double *); @@ -41,8 +41,8 @@ class AtomVecHybrid : public AtomVec { int size_restart_one(int) {return 0;} int pack_restart(int, double *); int unpack_restart(double *); - void create_atom(int, double, double, double, int); - void data_atom(double, double, double, int, char **, int); + void create_atom(int, double *, int); + void data_atom(double *, int, char **, int); void data_vel(int, char *, int); void data_params(int); int memory_usage(); diff --git a/src/comm.cpp b/src/comm.cpp index fab58a1dd7..be2116f190 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -11,10 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "mpi.h" +#include "math.h" #include "string.h" #include "stdio.h" #include "stdlib.h" +#include "mpi.h" #include "comm.h" #include "atom.h" #include "atom_vec.h" @@ -72,9 +73,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp) maxforward_pair = maxreverse_pair = 0; } -/* ---------------------------------------------------------------------- - destructor, free all memory -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ Comm::~Comm() { @@ -122,11 +121,13 @@ void Comm::init() } /* ---------------------------------------------------------------------- - setup 3d grid of procs based on box size + setup 3d grid of procs based on box size ------------------------------------------------------------------------- */ void Comm::set_procs() { + triclinic = domain->triclinic; + if (user_procgrid[0] == 0) procs2box(); else { procgrid[0] = user_procgrid[0]; @@ -151,6 +152,10 @@ void Comm::set_procs() MPI_Cart_shift(cartesian,2,1,&procneigh[2][0],&procneigh[2][1]); MPI_Comm_free(&cartesian); + // can't set lamda box params until procs are assigned + + if (triclinic) domain->set_lamda_box(); + if (me == 0) { if (screen) fprintf(screen," %d by %d by %d processor grid\n", procgrid[0],procgrid[1],procgrid[2]); @@ -166,25 +171,47 @@ void Comm::set_procs() void Comm::setup() { + // for triclinic, tilted planes are neigh cutoff apart + // but cutneigh is distance between those planes along xyz (non-tilted) axes + + double cutneigh[3]; + double *prd,*prd_border,*sublo,*subhi; + + if (triclinic == 0) { + prd = prd_border = domain->prd; + sublo = domain->sublo; + subhi = domain->subhi; + cutneigh[0] = cutneigh[1] = cutneigh[2] = neighbor->cutneigh; + } else { + prd = domain->prd; + prd_border = domain->prd_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + double *h_inv = domain->h_inv; + double length; + length = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]); + cutneigh[0] = neighbor->cutneigh * length; + length = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]); + cutneigh[1] = neighbor->cutneigh * length; + length = h_inv[2]; + cutneigh[2] = neighbor->cutneigh * length; + } + // need = # of procs I need atoms from in each dim - - need[0] = static_cast - (neighbor->cutneigh * procgrid[0] / domain->xprd) + 1; - need[1] = static_cast - (neighbor->cutneigh * procgrid[1] / domain->yprd) + 1; - need[2] = static_cast - (neighbor->cutneigh * procgrid[2] / domain->zprd) + 1; - // for 2d, don't communicate in z + need[0] = static_cast (cutneigh[0] * procgrid[0] / prd_border[0]) + 1; + need[1] = static_cast (cutneigh[1] * procgrid[1] / prd_border[1]) + 1; + need[2] = static_cast (cutneigh[2] * procgrid[2] / prd_border[2]) + 1; if (force->dimension == 2) need[2] = 0; // if non-periodic, do not communicate further than procgrid-1 away // this enables very large cutoffs in non-periodic systems - if (domain->xperiodic == 0) need[0] = MIN(need[0],procgrid[0]-1); - if (domain->yperiodic == 0) need[1] = MIN(need[1],procgrid[1]-1); - if (domain->zperiodic == 0) need[2] = MIN(need[2],procgrid[2]-1); + int *periodicity = domain->periodicity; + if (periodicity[0] == 0) need[0] = MIN(need[0],procgrid[0]-1); + if (periodicity[1] == 0) need[1] = MIN(need[1],procgrid[1]-1); + if (periodicity[2] == 0) need[2] = MIN(need[2],procgrid[2]-1); // allocate comm memory @@ -201,47 +228,62 @@ void Comm::setup() // set slablo > slabhi for swaps across non-periodic boundaries // this insures no atoms are swapped // only for procs owning sub-box at non-periodic end of global box - // pbc_flags = add-on factor for atoms sent across a periodic global boundary - // 0 = not across a boundary - // 1 = add box-length to coord when sending - // -1 = subtract box-length from coord when sending + // pbc_flag: 0 = not across a boundary, 1 = yes across a boundary + // pbc_dist/border: factors to add to atom coords across PBC for comm/borders + // for triclinic, slablo/hi and pbc_border will be used in lamda (0-1) coords // 1st part of if statement is sending to the west/south/down // 2nd part of if statement is sending to the east/north/up int iswap = 0; for (int dim = 0; dim < 3; dim++) { for (int ineed = 0; ineed < 2*need[dim]; ineed++) { - pbc_flags[iswap][0] = 0; - pbc_flags[iswap][1] = 0; - pbc_flags[iswap][2] = 0; - pbc_flags[iswap][3] = 0; + pbc_flag[iswap] = 0; + pbc_dist[iswap][0] = pbc_dist[iswap][1] = pbc_dist[iswap][2] = 0.0; + pbc_dist_border[iswap][0] = pbc_dist_border[iswap][1] = + pbc_dist_border[iswap][2] = 0.0; if (ineed % 2 == 0) { sendproc[iswap] = procneigh[dim][0]; recvproc[iswap] = procneigh[dim][1]; if (ineed < 2) slablo[iswap] = -BIG; - else slablo[iswap] = 0.5 * (domain->sublo[dim] + domain->subhi[dim]); - slabhi[iswap] = domain->sublo[dim] + neighbor->cutneigh; + else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]); + slabhi[iswap] = sublo[dim] + cutneigh[dim]; if (myloc[dim] == 0) { - if (domain->periodicity[dim] == 0) + if (periodicity[dim] == 0) slabhi[iswap] = slablo[iswap] - 1.0; else { - pbc_flags[iswap][0] = 1; - pbc_flags[iswap][dim+1] = 1; + pbc_flag[iswap] = 1; + pbc_dist[iswap][dim] = prd[dim]; + pbc_dist_border[iswap][dim] = prd_border[dim]; + if (triclinic) { + if (dim == 1) pbc_dist[iswap][0] += domain->xy; + else if (dim == 2) { + pbc_dist[iswap][0] += domain->xz; + pbc_dist[iswap][1] += domain->yz; + } + } } } } else { sendproc[iswap] = procneigh[dim][1]; recvproc[iswap] = procneigh[dim][0]; - slablo[iswap] = domain->subhi[dim] - neighbor->cutneigh; + slablo[iswap] = subhi[dim] - cutneigh[dim]; if (ineed < 2) slabhi[iswap] = BIG; - else slabhi[iswap] = 0.5 * (domain->sublo[dim] + domain->subhi[dim]); + else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]); if (myloc[dim] == procgrid[dim]-1) { - if (domain->periodicity[dim] == 0) + if (periodicity[dim] == 0) slabhi[iswap] = slablo[iswap] - 1.0; else { - pbc_flags[iswap][0] = 1; - pbc_flags[iswap][dim+1] = -1; + pbc_flag[iswap] = 1; + pbc_dist[iswap][dim] = -prd[dim]; + pbc_dist_border[iswap][dim] = -prd_border[dim]; + if (triclinic) { + if (dim == 1) pbc_dist[iswap][0] -= domain->xy; + else if (dim == 2) { + pbc_dist[iswap][0] -= domain->xz; + pbc_dist[iswap][1] -= domain->yz; + } + } } } } @@ -277,14 +319,14 @@ void Comm::communicate() MPI_Irecv(buf,size_comm_recv[iswap],MPI_DOUBLE, recvproc[iswap],0,world,&request); n = avec->pack_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flags[iswap]); + buf_send,pbc_flag[iswap],pbc_dist[iswap]); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); } else { MPI_Irecv(buf_recv,size_comm_recv[iswap],MPI_DOUBLE, recvproc[iswap],0,world,&request); n = avec->pack_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flags[iswap]); + buf_send,pbc_flag[iswap],pbc_dist[iswap]); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); MPI_Wait(&request,&status); avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_recv); @@ -294,10 +336,11 @@ void Comm::communicate() if (comm_x_only) { if (sendnum[iswap]) n = avec->pack_comm(sendnum[iswap],sendlist[iswap], - x[firstrecv[iswap]],pbc_flags[iswap]); + x[firstrecv[iswap]], + pbc_flag[iswap],pbc_dist[iswap]); } else { n = avec->pack_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flags[iswap]); + buf_send,pbc_flag[iswap],pbc_dist[iswap]); avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send); } } @@ -363,6 +406,7 @@ void Comm::reverse_communicate() can happen if atom moves outside of non-periodic bounary or if atom moves more than one proc away this routine called before every reneighboring + for triclinic, atoms must be in lamda coords (0-1) before exchange is called ------------------------------------------------------------------------- */ void Comm::exchange() @@ -370,7 +414,7 @@ void Comm::exchange() int i,m,nsend,nrecv,nrecv1,nrecv2,nlocal; double lo,hi,value; double **x; - double *buf; + double *sublo,*subhi,*buf; MPI_Request request; MPI_Status status; AtomVec *avec = atom->avec; @@ -379,6 +423,16 @@ void Comm::exchange() if (map_style) atom->map_clear(); + // subbox bounds for orthogonal or triclinic + + if (triclinic == 0) { + sublo = domain->sublo; + subhi = domain->subhi; + } else { + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + // loop over dimensions for (int dim = 0; dim < 3; dim++) { @@ -387,8 +441,8 @@ void Comm::exchange() // when atom is deleted, fill it in with last atom x = atom->x; - lo = domain->sublo[dim]; - hi = domain->subhi[dim]; + lo = sublo[dim]; + hi = subhi[dim]; nlocal = atom->nlocal; i = nsend = 0; @@ -457,6 +511,7 @@ void Comm::exchange() this does equivalent of a communicate (so don't need to explicitly call communicate routine on reneighboring timestep) this routine is called before every reneighboring + for triclinic, atoms must be in lamda coords (0-1) before borders is called ------------------------------------------------------------------------- */ void Comm::borders() @@ -509,7 +564,8 @@ void Comm::borders() if (nsend*size_border > maxsend) grow_send(nsend*size_border,0); - n = avec->pack_border(nsend,sendlist[iswap],buf_send,pbc_flags[iswap]); + n = avec->pack_border(nsend,sendlist[iswap],buf_send, + pbc_flag[iswap],pbc_dist_border[iswap]); // swap atoms with other proc // put incoming ghosts at end of my atom arrays @@ -577,7 +633,7 @@ void Comm::comm_pair(Pair *pair) // pack buffer n = pair->pack_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flags[iswap]); + buf_send,pbc_flag[iswap],pbc_dist[iswap]); // exchange with another proc // if self, set recv buffer to send buffer @@ -646,7 +702,7 @@ void Comm::comm_fix(Fix *fix) // pack buffer n = fix->pack_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flags[iswap]); + buf_send,pbc_flag[iswap],pbc_dist[iswap]); // exchange with another proc // if self, set recv buffer to send buffer @@ -715,7 +771,7 @@ void Comm::comm_compute(Compute *compute) // pack buffer n = compute->pack_comm(sendnum[iswap],sendlist[iswap], - buf_send,pbc_flags[iswap]); + buf_send,pbc_flag[iswap],pbc_dist[iswap]); // exchange with another proc // if self, set recv buffer to send buffer @@ -770,23 +826,38 @@ void Comm::reverse_comm_compute(Compute *compute) /* ---------------------------------------------------------------------- assign nprocs to 3d xprd,yprd,zprd box so as to minimize surface area + area = surface area of each of 3 faces of simulation box + for triclinic, area = cross product of 2 edge vectors stored in h matrix ------------------------------------------------------------------------- */ void Comm::procs2box() { - int ipx,ipy,ipz,nremain; - double boxx,boxy,boxz,surf; + double area[3]; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + if (triclinic == 0) { + area[0] = domain->xprd * domain->yprd; + area[1] = domain->xprd * domain->zprd; + area[2] = domain->yprd * domain->zprd; + } else { + double *h = domain->h; + double x,y,z; + cross(h[0],0.0,0.0,h[5],h[1],0.0,x,y,z); + area[0] = sqrt(x*x + y*y + z*z); + cross(h[0],0.0,0.0,h[4],h[3],h[2],x,y,z); + area[1] = sqrt(x*x + y*y + z*z); + cross(h[5],h[1],0.0,h[4],h[3],h[2],x,y,z); + area[2] = sqrt(x*x + y*y + z*z); + } - double bestsurf = 2.0 * (xprd*yprd + yprd*zprd + zprd*xprd); + double bestsurf = 2.0 * (area[0]+area[1]+area[2]); // loop thru all possible factorizations of nprocs // surf = surface area of a proc sub-domain // for 2d, insure ipz = 1 + int ipx,ipy,ipz,nremain; + double surf; + ipx = 1; while (ipx <= nprocs) { if (nprocs % ipx == 0) { @@ -796,10 +867,7 @@ void Comm::procs2box() if (nremain % ipy == 0) { ipz = nremain/ipy; if (force->dimension == 3 || ipz == 1) { - boxx = xprd/ipx; - boxy = yprd/ipy; - boxz = zprd/ipz; - surf = boxx*boxy + boxy*boxz + boxz*boxx; + surf = area[0]/ipx/ipy + area[1]/ipx/ipz + area[2]/ipy/ipz; if (surf < bestsurf) { bestsurf = surf; procgrid[0] = ipx; @@ -815,6 +883,19 @@ void Comm::procs2box() } } +/* ---------------------------------------------------------------------- + vector cross product: c = a x b +------------------------------------------------------------------------- */ + +void Comm::cross(double ax, double ay, double az, + double bx, double by, double bz, + double &cx, double &cy, double &cz) +{ + cx = ay*bz - az*by; + cy = az*bx - ax*bz; + cz = ax*by - ay*bx; +} + /* ---------------------------------------------------------------------- realloc the size of the send buffer as needed with BUFFACTOR & BUFEXTRA if flag = 1, realloc @@ -896,7 +977,10 @@ void Comm::allocate_swap(int n) slablo = (double *) memory->smalloc(n*sizeof(double),"comm:slablo"); slabhi = (double *) memory->smalloc(n*sizeof(double),"comm:slabhi"); firstrecv = (int *) memory->smalloc(n*sizeof(int),"comm:firstrecv"); - pbc_flags = (int **) memory->create_2d_int_array(n,4,"comm:pbc_flags"); + pbc_flag = (int *) memory->smalloc(n*sizeof(int),"comm:pbc_flag"); + pbc_dist = (double **) memory->create_2d_double_array(n,3,"comm:pbc_dist"); + pbc_dist_border = (double **) + memory->create_2d_double_array(n,3,"comm:pbc_dist_border"); } /* ---------------------------------------------------------------------- @@ -915,7 +999,9 @@ void Comm::free_swap() memory->sfree(slablo); memory->sfree(slabhi); memory->sfree(firstrecv); - memory->destroy_2d_int_array(pbc_flags); + memory->sfree(pbc_flag); + memory->destroy_2d_double_array(pbc_dist); + memory->destroy_2d_double_array(pbc_dist_border); } /* ---------------------------------------------------------------------- diff --git a/src/comm.h b/src/comm.h index 232a6aec31..a6eaf8ce04 100644 --- a/src/comm.h +++ b/src/comm.h @@ -51,6 +51,7 @@ class Comm : protected Pointers { int memory_usage(); // tally memory usage private: + int triclinic; // 0 if domain is orthog, 1 if triclinic int maxswap; // max # of swaps memory is allocated for int *sendnum,*recvnum; // # of atoms to send/recv in each swap int *sendproc,*recvproc; // proc to send/recv to/from at each swap @@ -58,9 +59,9 @@ class Comm : protected Pointers { int *size_reverse_send; // # to send in each reverse comm int *size_reverse_recv; // # to recv in each reverse comm double *slablo,*slabhi; // bounds of slab to send at each swap - int **pbc_flags; // flags for sending atoms thru PBC - // [0] = 1 if any dim is across PBC - // [123] = 1 if dim 012 is across PBC + int *pbc_flag; // flags for sending atoms thru PBC + double **pbc_dist; // distance to adjust atoms coords for PBC + double **pbc_dist_border; // PBC distance for borders() int comm_x_only,comm_f_only; // 1 if only exchange x,f in for/rev comm int map_style; // non-0 if global->local mapping is done @@ -74,6 +75,9 @@ class Comm : protected Pointers { int maxforward,maxreverse; // max # of datums in forward/reverse comm void procs2box(); // map procs to 3d box + void cross(double, double, double, + double, double, double, + double &, double &, double &); // cross product void grow_send(int,int); // reallocate send buffer void grow_recv(int); // free/allocate recv buffer void grow_list(int, int); // reallocate one sendlist diff --git a/src/compute.h b/src/compute.h index f0d9222e04..5e6c3cca86 100644 --- a/src/compute.h +++ b/src/compute.h @@ -55,7 +55,7 @@ class Compute : protected Pointers { virtual void compute_vector() {} virtual void compute_peratom() {} - virtual int pack_comm(int, int *, double *, int *) {return 0;} + virtual int pack_comm(int, int *, double *, int, double *) {return 0;} virtual void unpack_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index ad7c53226e..e422894525 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -54,8 +54,8 @@ void CreateAtoms::command(int narg, char **arg) error->all("Invalid atom type in create_atoms command"); for (int i = 0; i < nbasis; i++) basistype[i] = itype; - regionflag = -1; - nhybrid = 0; + int regionflag = -1; + int nhybrid = 0; int iarg = 1; while (iarg < narg) { @@ -93,30 +93,42 @@ void CreateAtoms::command(int narg, char **arg) // convert 8 corners of my sub-box from box coords to lattice coords // min to max = bounding box around the pts in lattice space - subxlo = domain->subxlo; - subxhi = domain->subxhi; - subylo = domain->subylo; - subyhi = domain->subyhi; - subzlo = domain->subzlo; - subzhi = domain->subzhi; + int triclinic = domain->triclinic; + double *bboxlo,*bboxhi; + + if (triclinic == 0) { + bboxlo = domain->sublo; + bboxhi = domain->subhi; + } else { + bboxlo = domain->sublo_bound; + bboxhi = domain->subhi_bound; + } double xmin,ymin,zmin,xmax,ymax,zmax; xmin = ymin = zmin = BIG; xmax = ymax = zmax = -BIG; - domain->lattice->bbox(1,subxlo,subylo,subzlo,xmin,ymin,zmin,xmax,ymax,zmax); - domain->lattice->bbox(1,subxhi,subylo,subzlo,xmin,ymin,zmin,xmax,ymax,zmax); - domain->lattice->bbox(1,subxlo,subyhi,subzlo,xmin,ymin,zmin,xmax,ymax,zmax); - domain->lattice->bbox(1,subxhi,subyhi,subzlo,xmin,ymin,zmin,xmax,ymax,zmax); - domain->lattice->bbox(1,subxlo,subylo,subzhi,xmin,ymin,zmin,xmax,ymax,zmax); - domain->lattice->bbox(1,subxhi,subylo,subzhi,xmin,ymin,zmin,xmax,ymax,zmax); - domain->lattice->bbox(1,subxlo,subyhi,subzhi,xmin,ymin,zmin,xmax,ymax,zmax); - domain->lattice->bbox(1,subxhi,subyhi,subzhi,xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxlo[2], + xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxlo[2], + xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxlo[2], + xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxlo[2], + xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxhi[2], + xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxhi[2], + xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxhi[2], + xmin,ymin,zmin,xmax,ymax,zmax); + domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxhi[2], + xmin,ymin,zmin,xmax,ymax,zmax); - // ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my sub-box - // overlap = any part of a unit cell (face,edge,pt) in common with my sub-box - // in lattice space, sub-box is a tilted box - // but bbox of sub-box is aligned with lattice axes + // ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my subbox + // overlap = any part of a unit cell (face,edge,pt) in common with my subbox + // in lattice space, subbox is a tilted box + // but bbox of subbox is aligned with lattice axes // so ilo:khi unit cells should completely tile bounding box // decrement lo values if min < 0, since static_cast(-1.5) = -1 @@ -132,25 +144,84 @@ void CreateAtoms::command(int narg, char **arg) if (ymin < 0.0) jlo--; if (zmin < 0.0) klo--; + // set bounds for my proc + // if periodic and I am lo/hi proc, adjust bounds by EPSILON + // on lower boundary, allows triclinic atoms just outside box to be added + // on upper boundary, prevents atoms with lower images from being added + + double sublo[3],subhi[3]; + + if (triclinic == 0) { + sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0]; + sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1]; + sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2]; + } else { + sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0]; + sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1]; + sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2]; + } + + if (domain->xperiodic) { + if (triclinic && comm->myloc[0] == 0) sublo[0] -= EPSILON; + if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] -= EPSILON; + } + if (domain->yperiodic) { + if (triclinic && comm->myloc[1] == 0) sublo[1] -= EPSILON; + if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] -= EPSILON; + } + if (domain->zperiodic) { + if (triclinic && comm->myloc[2] == 0) sublo[2] -= EPSILON; + if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] -= EPSILON; + } + // iterate on 3d periodic lattice using loop bounds // invoke add_atom for nbasis atoms in each unit cell - // add_atom converts lattice coords to box coords, checks if in my sub-box - - boxxhi = domain->boxxhi; - boxyhi = domain->boxyhi; - boxzhi = domain->boxzhi; + // converts lattice coords to box coords + // add atom if it meets all criteria double natoms_previous = atom->natoms; int nlocal_previous = atom->nlocal; double **basis = domain->lattice->basis; + double x[3],lamda[3]; + double *coord; int i,j,k,m; for (k = klo; k <= khi; k++) for (j = jlo; j <= jhi; j++) for (i = ilo; i <= ihi; i++) - for (m = 0; m < nbasis; m++) - add_atom(basistype[m],i+basis[m][0],j+basis[m][1],k+basis[m][2]); + for (m = 0; m < nbasis; m++) { + + x[0] = i + basis[m][0]; + x[1] = j + basis[m][1]; + x[2] = k + basis[m][2]; + + // convert from lattice coords to box coords + + domain->lattice->lattice2box(x[0],x[1],x[2]); + + // if a region was specified, test if atom is in it + + if (regionflag >= 0) + if (!domain->regions[regionflag]->match(x[0],x[1],x[2])) continue; + + // test if atom is in my subbox + + if (triclinic) { + domain->x2lamda(x,lamda); + coord = lamda; + } else coord = x; + + if (coord[0] < sublo[0] || coord[0] >= subhi[0] || + coord[1] < sublo[1] || coord[1] >= subhi[1] || + coord[2] < sublo[2] || coord[2] >= subhi[2]) continue; + + // add the atom to my list of atoms + + atom->avec->create_atom(basistype[m],x,nhybrid); + } + + // clean up delete [] basistype; @@ -189,36 +260,3 @@ void CreateAtoms::command(int narg, char **arg) } } } - -/* ---------------------------------------------------------------------- - add an atom of type at lattice coords x,y,z if it meets all criteria -------------------------------------------------------------------------- */ - -void CreateAtoms::add_atom(int ntype, double x, double y, double z) -{ - // convert from lattice coords to box coords - - domain->lattice->lattice2box(x,y,z); - - // if a region was specified, test if atom is in it - - if (regionflag >= 0) - if (!domain->regions[regionflag]->match(x,y,z)) return; - - // test if atom is in my subbox - - if (x < subxlo || x >= subxhi || - y < subylo || y >= subyhi || - z < subzlo || z >= subzhi) return; - - // don't put atoms within EPSILON of upper periodic boundary - // else may overlap image atom at lower boundary - - if (domain->xperiodic && fabs(x-boxxhi) < EPSILON) return; - if (domain->yperiodic && fabs(y-boxyhi) < EPSILON) return; - if (domain->zperiodic && fabs(z-boxzhi) < EPSILON) return; - - // add the atom to my list of atoms - - atom->avec->create_atom(ntype,x,y,z,nhybrid); -} diff --git a/src/create_atoms.h b/src/create_atoms.h index a5b3b51314..c91ab74e97 100644 --- a/src/create_atoms.h +++ b/src/create_atoms.h @@ -22,13 +22,6 @@ class CreateAtoms : protected Pointers { public: CreateAtoms(class LAMMPS *); void command(int, char **); - - private: - int regionflag,nhybrid; - double boxxhi,boxyhi,boxzhi; - double subxlo,subxhi,subylo,subyhi,subzlo,subzhi; - - void add_atom(int, double, double, double); }; } diff --git a/src/create_box.cpp b/src/create_box.cpp index 19c88aaf51..570b125ec8 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -18,6 +18,7 @@ #include "force.h" #include "domain.h" #include "region.h" +#include "region_prism.h" #include "comm.h" #include "update.h" #include "error.h" @@ -39,6 +40,8 @@ void CreateBox::command(int narg, char **arg) if (force->dimension == 2 && domain->zperiodic == 0) error->all("Cannot run 2d simulation with nonperiodic Z dimension"); + domain->box_exist = 1; + // find region ID int iregion; @@ -50,26 +53,33 @@ void CreateBox::command(int narg, char **arg) if (domain->regions[iregion]->interior == 0) error->all("Create_box region must be of type inside"); - // set global box from region extent + // if region not prism: + // setup orthogonal domain + // set simulation domain from region extent + // if region is prism: + // seutp triclinic domain + // set simulation domain params from prism params - domain->boxxlo = domain->regions[iregion]->extent_xlo; - domain->boxxhi = domain->regions[iregion]->extent_xhi; - domain->boxylo = domain->regions[iregion]->extent_ylo; - domain->boxyhi = domain->regions[iregion]->extent_yhi; - domain->boxzlo = domain->regions[iregion]->extent_zlo; - domain->boxzhi = domain->regions[iregion]->extent_zhi; - - domain->box_exist = 1; + if (strcmp(domain->regions[iregion]->style,"prism") != 0) { + domain->boxxlo = domain->regions[iregion]->extent_xlo; + domain->boxxhi = domain->regions[iregion]->extent_xhi; + domain->boxylo = domain->regions[iregion]->extent_ylo; + domain->boxyhi = domain->regions[iregion]->extent_yhi; + domain->boxzlo = domain->regions[iregion]->extent_zlo; + domain->boxzhi = domain->regions[iregion]->extent_zhi; - if (comm->me == 0) { - if (screen) - fprintf(screen,"Created box = (%g %g %g) to (%g %g %g)\n", - domain->boxxlo,domain->boxylo,domain->boxzlo, - domain->boxxhi,domain->boxyhi,domain->boxzhi); - if (logfile) - fprintf(logfile,"Created box = (%g %g %g) to (%g %g %g)\n", - domain->boxxlo,domain->boxylo,domain->boxzlo, - domain->boxxhi,domain->boxyhi,domain->boxzhi); + } else { + domain->triclinic = 1; + RegPrism *region = (RegPrism *) domain->regions[iregion]; + domain->boxxlo = region->xlo; + domain->boxxhi = region->xhi; + domain->boxylo = region->ylo; + domain->boxyhi = region->yhi; + domain->boxzlo = region->zlo; + domain->boxzhi = region->zhi; + domain->xy = region->xy; + domain->xz = region->xz; + domain->yz = region->yz; } // if molecular, zero out topology info @@ -100,6 +110,7 @@ void CreateBox::command(int narg, char **arg) atom->allocate_type_arrays(); + domain->print_box("Created "); domain->set_initial_box(); domain->set_global_box(); comm->set_procs(); diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index d21893c1c9..2eecabd960 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -196,13 +196,15 @@ void DeleteAtoms::delete_overlap(int narg, char **arg, int *list) // setup domain, communication and neighboring // acquire ghosts - + + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); comm->exchange(); comm->borders(); + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); // call to build() forces memory allocation for neighbor lists // build half list explicitly if build() doesn't do it diff --git a/src/delete_bonds.cpp b/src/delete_bonds.cpp index 79b964aebd..02f644e5a4 100644 --- a/src/delete_bonds.cpp +++ b/src/delete_bonds.cpp @@ -104,11 +104,13 @@ void DeleteBonds::command(int narg, char **arg) // border swap to insure type and mask is current for off-proc atoms // enforce PBC before in case atoms are outside box + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); comm->exchange(); comm->borders(); + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); // set topology interactions either off or on // criteria for an interaction to potentially be changed (set flag = 1) diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp index a88279ef84..dad3c57b0f 100644 --- a/src/displace_atoms.cpp +++ b/src/displace_atoms.cpp @@ -165,10 +165,12 @@ void DisplaceAtoms::command(int narg, char **arg) // move atoms to new processors // enforce PBC before in case atoms are outside box + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); comm->exchange(); + if (domain->triclinic) domain->lamda2x(atom->nlocal); // check if any atoms were lost diff --git a/src/domain.cpp b/src/domain.cpp index fa62258e1c..699f612c68 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -58,8 +58,14 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp) boundary[1][0] = boundary[1][1] = 0; boundary[2][0] = boundary[2][1] = 0; + triclinic = 0; boxxlo = boxylo = boxzlo = -0.5; boxxhi = boxyhi = boxzhi = 0.5; + xy = xz = yz = 0.0; + + prd_lamda[0] = prd_lamda[1] = prd_lamda[2] = 1.0; + boxlo_lamda[0] = boxlo_lamda[1] = boxlo_lamda[2] = 0.0; + boxhi_lamda[0] = boxhi_lamda[1] = boxhi_lamda[2] = 1.0; lattice = NULL; nregion = maxregion = 0; @@ -79,26 +85,47 @@ Domain::~Domain() void Domain::init() { - // set box_change if box dimensions ever change - // due to shrink-wrapping, fix nph, npt, volume/rescale, uniaxial + // set box_change if box dimensions/shape ever changes + // due to shrink-wrapping, fixes that change volume (npt, vol/rescale, etc) box_change = 0; if (nonperiodic == 2) box_change = 1; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"nph") == 0) box_change = 1; - if (strcmp(modify->fix[i]->style,"npt") == 0) box_change = 1; - if (strcmp(modify->fix[i]->style,"volume/rescale") == 0) box_change = 1; - if (strcmp(modify->fix[i]->style,"uniaxial") == 0) box_change = 1; - } + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->box_change) box_change = 1; } /* ---------------------------------------------------------------------- - set initial global box from boxlo/hi (already set by caller) - adjust for shrink-wrapped dims + set initial global box + assumes boxlo/hi and triclinic tilts are already set ------------------------------------------------------------------------- */ void Domain::set_initial_box() { + // error checks for orthogonal and triclinic domains + + if (boxxlo >= boxxhi || boxylo >= boxyhi || boxzlo >= boxzhi) + error->one("Box bounds are invalid"); + + if (triclinic) { + if (force->dimension == 2 && (xz != 0.0 || yz != 0.0)) + error->all("Cannot skew triclinic box in z for 2d simulation"); + if (xy != 0.0 && (!xperiodic || !yperiodic)) + error->all("Triclinic box must be periodic in skewed dimensions"); + if (xz != 0.0 && (!xperiodic || !zperiodic)) + error->all("Triclinic box must be periodic in skewed dimensions"); + if (yz != 0.0 && (!yperiodic || !zperiodic)) + error->all("Triclinic box must be periodic in skewed dimensions"); + + if (fabs(xy/(boxyhi-boxylo)) > 0.5) + error->all("Triclinic box skew is too large"); + if (fabs(xz/(boxzhi-boxzlo)) > 0.5) + error->all("Triclinic box skew is too large"); + if (fabs(yz/(boxzhi-boxzlo)) > 0.5) + error->all("Triclinic box skew is too large"); + } + + // adjust box lo/hi for shrink-wrapped dims + if (boundary[0][0] == 2) boxxlo -= SMALL; else if (boundary[0][0] == 3) minxlo = boxxlo; if (boundary[0][1] == 2) boxxhi += SMALL; @@ -116,54 +143,132 @@ void Domain::set_initial_box() } /* ---------------------------------------------------------------------- - set global prd, prd_half, prd[], boxlo/hi[] from boxlo/hi + set global box params + assumes boxlo/hi and triclinic tilts are already set ------------------------------------------------------------------------- */ void Domain::set_global_box() { - xprd = boxxhi - boxxlo; - yprd = boxyhi - boxylo; - zprd = boxzhi - boxzlo; + prd[0] = xprd = boxxhi - boxxlo; + prd[1] = yprd = boxyhi - boxylo; + prd[2] = zprd = boxzhi - boxzlo; xprd_half = 0.5*xprd; yprd_half = 0.5*yprd; zprd_half = 0.5*zprd; - prd[0] = xprd; prd[1] = yprd; prd[2] = zprd; boxlo[0] = boxxlo; boxlo[1] = boxylo; boxlo[2] = boxzlo; boxhi[0] = boxxhi; boxhi[1] = boxyhi; boxhi[2] = boxzhi; + + if (triclinic) { + h[0] = xprd; + h[1] = yprd; + h[2] = zprd; + h[3] = yz; + h[4] = xz; + h[5] = xy; + + h_inv[0] = 1.0/h[0]; + h_inv[1] = 1.0/h[1]; + h_inv[2] = 1.0/h[2]; + h_inv[3] = -h[3] / (h[1]*h[2]); + h_inv[4] = (h[3]*h[5] - h[1]*h[4]) / (h[0]*h[1]*h[2]); + h_inv[5] = -h[5] / (h[0]*h[1]); + + boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy); + boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz); + boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz); + boxlo_bound[2] = boxlo[2]; + + boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy); + boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz); + boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz); + boxhi_bound[2] = boxhi[2]; + } } /* ---------------------------------------------------------------------- - set local subxlo-subzhi, sublo/hi[] from global boxxlo-boxzhi and proc grid + set lamda box params, only need be done one time + assumes global box is defined and proc assignment has been made by comm + for uppermost proc, insure subhi = 1.0 (in case round-off occurs) +------------------------------------------------------------------------- */ + +void Domain::set_lamda_box() +{ + int *myloc = comm->myloc; + int *procgrid = comm->procgrid; + + sublo_lamda[0] = 1.0*myloc[0] / procgrid[0]; + sublo_lamda[1] = 1.0*myloc[1] / procgrid[1]; + sublo_lamda[2] = 1.0*myloc[2] / procgrid[2]; + + if (myloc[0] < procgrid[0]-1) + subhi_lamda[0] = 1.0*(myloc[0]+1) / procgrid[0]; + else subhi_lamda[0] = 1.0; + + if (myloc[1] < procgrid[1]-1) + subhi_lamda[1] = 1.0*(myloc[1]+1) / procgrid[1]; + else subhi_lamda[1] = 1.0; + + if (myloc[2] < procgrid[2]-1) + subhi_lamda[2] = 1.0*(myloc[2]+1) / procgrid[2]; + else subhi_lamda[2] = 1.0; +} + +/* ---------------------------------------------------------------------- + set local subbox params + assumes global box is defined and proc assignment has been made for uppermost proc, insure subhi = boxhi (in case round-off occurs) + for triclinic, lo/hi_bound is a bbox around all 8 tilted pts of sub-box ------------------------------------------------------------------------- */ void Domain::set_local_box() { - subxlo = boxxlo + comm->myloc[0] * xprd / comm->procgrid[0]; - if (comm->myloc[0] < comm->procgrid[0]-1) - subxhi = boxxlo + (comm->myloc[0]+1) * xprd / comm->procgrid[0]; - else subxhi = boxxhi; + int *myloc = comm->myloc; + int *procgrid = comm->procgrid; - subylo = boxylo + comm->myloc[1] * yprd / comm->procgrid[1]; - if (comm->myloc[1] < comm->procgrid[1]-1) - subyhi = boxylo + (comm->myloc[1]+1) * yprd / comm->procgrid[1]; - else subyhi = boxyhi; + if (domain->triclinic == 0) { + sublo[0] = boxlo[0] + myloc[0] * xprd / procgrid[0]; + if (myloc[0] < procgrid[0]-1) + subhi[0] = boxlo[0] + (myloc[0]+1) * xprd / procgrid[0]; + else subhi[0] = boxhi[0]; + + sublo[1] = boxlo[1] + myloc[1] * yprd / procgrid[1]; + if (myloc[1] < procgrid[1]-1) + subhi[1] = boxlo[1] + (myloc[1]+1) * yprd / procgrid[1]; + else subhi[1] = boxhi[1]; - subzlo = boxzlo + comm->myloc[2] * zprd / comm->procgrid[2]; - if (comm->myloc[2] < comm->procgrid[2]-1) - subzhi = boxzlo + (comm->myloc[2]+1) * zprd / comm->procgrid[2]; - else subzhi = boxzhi; + sublo[2] = boxlo[2] + myloc[2] * zprd / procgrid[2]; + if (myloc[2] < procgrid[2]-1) + subhi[2] = boxlo[2] + (myloc[2]+1) * zprd / procgrid[2]; + else subhi[2] = boxhi[2]; - sublo[0] = subxlo; sublo[1] = subylo; sublo[2] = subzlo; - subhi[0] = subxhi; subhi[1] = subyhi; subhi[2] = subzhi; + } else { + sublo[0] = h[0]*sublo_lamda[0] + h[5]*sublo_lamda[1] + + h[4]*sublo_lamda[2] + boxlo[0]; + sublo[1] = h[1]*sublo_lamda[1] + h[3]*sublo_lamda[2] + boxlo[1]; + sublo[2] = h[2]*sublo_lamda[2] + boxlo[2]; + + subhi[0] = sublo[0] + xprd/procgrid[0]; + subhi[1] = sublo[1] + yprd/procgrid[1]; + subhi[2] = sublo[2] + zprd/procgrid[2]; + + sublo_bound[0] = MIN(sublo[0],sublo[0] + xy/procgrid[1]); + sublo_bound[0] = MIN(sublo_bound[0],sublo_bound[0] + xz/procgrid[2]); + sublo_bound[1] = MIN(sublo[1],sublo[1] + yz/procgrid[2]); + sublo_bound[2] = sublo[2]; + + subhi_bound[0] = MAX(subhi[0],subhi[0] + xy/procgrid[1]); + subhi_bound[0] = MAX(subhi_bound[0],subhi_bound[0] + xz/procgrid[2]); + subhi_bound[1] = MAX(subhi[1],subhi[1] + yz/procgrid[2]); + subhi_bound[2] = subhi[2]; + } } /* ---------------------------------------------------------------------- reset global & local boxes due to global box boundary changes if shrink-wrapped, determine atom extent and reset boxlo/hi - set global & local boxes from new boxlo/hi values + shrink-wrapping only occurs in non-periodic, non-triclinic dims ------------------------------------------------------------------------- */ void Domain::reset_box() @@ -199,6 +304,7 @@ void Domain::reset_box() MPI_Allreduce(extent,all,6,MPI_DOUBLE,MPI_MAX,world); // in shrink-wrapped dims, set box by atom extent + // if set, observe min box size settings if (xperiodic == 0) { if (boundary[0][0] == 2) boxxlo = -all[0][0] - SMALL; @@ -229,9 +335,10 @@ void Domain::reset_box() /* ---------------------------------------------------------------------- enforce PBC and modify box image flags for each atom - called every reneighboring + called every reneighboring and by other commands that change atoms resulting coord must satisfy lo <= coord < hi MAX is important since coord - prd < lo can happen when coord = hi + for triclinic, atoms must be in lamda coords (0-1) before pbc is called image = 10 bits for each dimension increment/decrement in wrap-around fashion ------------------------------------------------------------------------- */ @@ -239,23 +346,34 @@ void Domain::reset_box() void Domain::pbc() { int i,idim,otherdims; + double *lo,*hi,*period; int nlocal = atom->nlocal; double **x = atom->x; int *image = atom->image; - if (xperiodic) { - for (i = 0; i < nlocal; i++) { - if (x[i][0] < boxxlo) { - x[i][0] += xprd; + if (triclinic == 0) { + lo = boxlo; + hi = boxhi; + period = prd; + } else { + lo = boxlo_lamda; + hi = boxhi_lamda; + period = prd_lamda; + } + + for (i = 0; i < nlocal; i++) { + if (xperiodic) { + if (x[i][0] < lo[0]) { + x[i][0] += period[0]; idim = image[i] & 1023; otherdims = image[i] ^ idim; idim--; idim &= 1023; image[i] = otherdims | idim; } - if (x[i][0] >= boxxhi) { - x[i][0] -= xprd; - x[i][0] = MAX(x[i][0],boxxlo); + if (x[i][0] >= hi[0]) { + x[i][0] -= period[0]; + x[i][0] = MAX(x[i][0],lo[0]); idim = image[i] & 1023; otherdims = image[i] ^ idim; idim++; @@ -263,21 +381,19 @@ void Domain::pbc() image[i] = otherdims | idim; } } - } - if (yperiodic) { - for (i = 0; i < nlocal; i++) { - if (x[i][1] < boxylo) { - x[i][1] += yprd; + if (yperiodic) { + if (x[i][1] < lo[1]) { + x[i][1] += period[1]; idim = (image[i] >> 10) & 1023; otherdims = image[i] ^ (idim << 10); idim--; idim &= 1023; image[i] = otherdims | (idim << 10); } - if (x[i][1] >= boxyhi) { - x[i][1] -= yprd; - x[i][1] = MAX(x[i][1],boxylo); + if (x[i][1] >= hi[1]) { + x[i][1] -= period[1]; + x[i][1] = MAX(x[i][1],lo[1]); idim = (image[i] >> 10) & 1023; otherdims = image[i] ^ (idim << 10); idim++; @@ -285,21 +401,19 @@ void Domain::pbc() image[i] = otherdims | (idim << 10); } } - } - if (zperiodic) { - for (i = 0; i < nlocal; i++) { - if (x[i][2] < boxzlo) { - x[i][2] += zprd; + if (zperiodic) { + if (x[i][2] < lo[2]) { + x[i][2] += period[2]; idim = image[i] >> 20; otherdims = image[i] ^ (idim << 20); idim--; idim &= 1023; image[i] = otherdims | (idim << 20); } - if (x[i][2] >= boxzhi) { - x[i][2] -= zprd; - x[i][2] = MAX(x[i][2],boxzlo); + if (x[i][2] >= hi[2]) { + x[i][2] -= period[2]; + x[i][2] = MAX(x[i][2],lo[2]); idim = image[i] >> 20; otherdims = image[i] ^ (idim << 20); idim++; @@ -313,53 +427,123 @@ void Domain::pbc() /* ---------------------------------------------------------------------- minimum image convention use 1/2 of box size as test + for triclinic, also add/subtract tilt factors in other dims as needed ------------------------------------------------------------------------- */ -void Domain::minimum_image(double *dx, double *dy, double *dz) +void Domain::minimum_image(double &dx, double &dy, double &dz) { - if (xperiodic) { - if (fabs(*dx) > xprd_half) { - if (*dx < 0.0) *dx += xprd; - else *dx -= xprd; + if (triclinic == 0) { + if (xperiodic) { + if (fabs(dx) > xprd_half) { + if (dx < 0.0) dx += xprd; + else dx -= xprd; + } } - } - if (yperiodic) { - if (fabs(*dy) > yprd_half) { - if (*dy < 0.0) *dy += yprd; - else *dy -= yprd; + if (yperiodic) { + if (fabs(dy) > yprd_half) { + if (dy < 0.0) dy += yprd; + else dy -= yprd; + } } - } - if (zperiodic) { - if (fabs(*dz) > zprd_half) { - if (*dz < 0.0) *dz += zprd; - else *dz -= zprd; + if (zperiodic) { + if (fabs(dz) > zprd_half) { + if (dz < 0.0) dz += zprd; + else dz -= zprd; + } + } + + } else { + if (zperiodic) { + if (fabs(dz) > zprd_half) { + if (dz < 0.0) { + dz += zprd; + dy += yz; + dx += xz; + } else { + dz -= zprd; + dy -= yz; + dx -= xz; + } + } + } + if (yperiodic) { + if (fabs(dy) > yprd_half) { + if (dy < 0.0) { + dy += yprd; + dx += xy; + } else { + dy -= yprd; + dx -= xy; + } + } + } + if (xperiodic) { + if (fabs(dx) > xprd_half) { + if (dx < 0.0) dx += xprd; + else dx -= xprd; + } } } } /* ---------------------------------------------------------------------- minimum image convention - use 1/2 of box size as test + use 1/2 of box size as test + for triclinic, also add/subtract tilt factors in other dims as needed ------------------------------------------------------------------------- */ -void Domain::minimum_image(double *x) +void Domain::minimum_image(double *delta) { - if (xperiodic) { - if (fabs(x[0]) > xprd_half) { - if (x[0] < 0.0) x[0] += xprd; - else x[0] -= xprd; + if (triclinic == 0) { + if (xperiodic) { + if (fabs(delta[0]) > xprd_half) { + if (delta[0] < 0.0) delta[0] += xprd; + else delta[0] -= xprd; + } } - } - if (yperiodic) { - if (fabs(x[1]) > yprd_half) { - if (x[1] < 0.0) x[1] += yprd; - else x[1] -= yprd; + if (yperiodic) { + if (fabs(delta[1]) > yprd_half) { + if (delta[1] < 0.0) delta[1] += yprd; + else delta[1] -= yprd; + } } - } - if (zperiodic) { - if (fabs(x[2]) > zprd_half) { - if (x[2] < 0.0) x[2] += zprd; - else x[2] -= zprd; + if (zperiodic) { + if (fabs(delta[2]) > zprd_half) { + if (delta[2] < 0.0) delta[2] += zprd; + else delta[2] -= zprd; + } + } + + } else { + if (zperiodic) { + if (fabs(delta[2]) > zprd_half) { + if (delta[2] < 0.0) { + delta[2] += zprd; + delta[1] += yz; + delta[0] += xz; + } else { + delta[2] -= zprd; + delta[1] -= yz; + delta[0] -= xz; + } + } + } + if (yperiodic) { + if (fabs(delta[1]) > yprd_half) { + if (delta[1] < 0.0) { + delta[1] += yprd; + delta[0] += xy; + } else { + delta[1] -= yprd; + delta[0] -= xy; + } + } + } + if (xperiodic) { + if (fabs(delta[0]) > xprd_half) { + if (delta[0] < 0.0) delta[0] += xprd; + else delta[0] -= xprd; + } } } } @@ -369,87 +553,113 @@ void Domain::minimum_image(double *x) adjust image accordingly resulting coord must satisfy lo <= coord < hi MAX is important since coord - prd < lo can happen when coord = hi + for triclinic, convert atom to lamda coords (0-1) before doing remap image = 10 bits for each dimension increment/decrement in wrap-around fashion ------------------------------------------------------------------------- */ -void Domain::remap(double &x, double &y, double &z, int &image) +void Domain::remap(double *x, int &image) { + double *lo,*hi,*period,*coord; + double lamda[3]; + + if (triclinic == 0) { + lo = boxlo; + hi = boxhi; + period = prd; + coord = x; + } else { + lo = boxlo_lamda; + hi = boxhi_lamda; + period = prd_lamda; + x2lamda(x,lamda); + coord = lamda; + } + if (xperiodic) { - while (x < boxxlo) { - x += xprd; + while (coord[0] < lo[0]) { + coord[0] += period[0]; int idim = image & 1023; int otherdims = image ^ idim; idim--; idim &= 1023; image = otherdims | idim; } - while (x >= boxxhi) { - x -= xprd; + while (coord[0] >= hi[0]) { + coord[0] -= period[0]; int idim = image & 1023; int otherdims = image ^ idim; idim++; idim &= 1023; image = otherdims | idim; } - x = MAX(x,boxxlo); + coord[0] = MAX(coord[0],lo[0]); } if (yperiodic) { - while (y < boxylo) { - y += yprd; + while (coord[1] < lo[1]) { + coord[1] += period[1]; int idim = (image >> 10) & 1023; int otherdims = image ^ (idim << 10); idim--; idim &= 1023; image = otherdims | (idim << 10); } - while (y >= boxyhi) { - y -= yprd; + while (coord[1] >= hi[1]) { + coord[1] -= period[1]; int idim = (image >> 10) & 1023; int otherdims = image ^ (idim << 10); idim++; idim &= 1023; image = otherdims | (idim << 10); } - y = MAX(y,boxylo); + coord[1] = MAX(coord[1],lo[1]); } if (zperiodic) { - while (z < boxzlo) { - z += zprd; + while (coord[2] < lo[2]) { + coord[2] += period[2]; int idim = image >> 20; int otherdims = image ^ (idim << 20); idim--; idim &= 1023; image = otherdims | (idim << 20); } - while (z >= boxzhi) { - z -= zprd; + while (coord[2] >= hi[2]) { + coord[2] -= period[2]; int idim = image >> 20; int otherdims = image ^ (idim << 20); idim++; idim &= 1023; image = otherdims | (idim << 20); } - z = MAX(z,boxzlo); + coord[2] = MAX(coord[2],lo[2]); } + + if (triclinic) lamda2x(coord,x); } /* ---------------------------------------------------------------------- unmap the point via image flags don't reset image flag + for triclinic, use h[] to add in tilt factors in other dims as needed ------------------------------------------------------------------------- */ -void Domain::unmap(double &x, double &y, double &z, int image) +void Domain::unmap(double *x, int image) { int xbox = (image & 1023) - 512; int ybox = (image >> 10 & 1023) - 512; int zbox = (image >> 20) - 512; - x = x + xbox*xprd; - y = y + ybox*yprd; - z = z + zbox*zprd; + if (triclinic == 0) { + x[0] += xbox*xprd; + x[1] += ybox*yprd; + x[2] += zbox*zprd; + } else { + x[0] += h[0]*xbox + h[5]*ybox + h[4]*zbox; + x[1] += h[1]*ybox + h[3]*zbox; + x[2] += h[2]*zbox; + } } /* ---------------------------------------------------------------------- @@ -551,3 +761,103 @@ void Domain::set_boundary(int narg, char **arg) boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2; } } + +/* ---------------------------------------------------------------------- + print box info, orthogonal or triclinic +------------------------------------------------------------------------- */ + +void Domain::print_box(char *str) +{ + if (comm->me == 0) { + if (screen) { + if (domain->triclinic == 0) + fprintf(screen,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n", + str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi); + else { + char *format = + "%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n"; + fprintf(screen,format, + str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi,xy,xz,yz); + } + } + if (logfile) { + if (triclinic == 0) + fprintf(logfile,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n", + str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi); + else { + char *format = + "%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n"; + fprintf(logfile,format, + str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi,xy,xz,yz); + } + } + } +} + +/* ---------------------------------------------------------------------- + convert triclinic 0-1 lamda coords to box coords for all N atoms + x = H lamda + x0; +------------------------------------------------------------------------- */ + +void Domain::lamda2x(int n) +{ + double **x = atom->x; + + for (int i = 0; i < n; i++) { + x[i][0] = h[0]*x[i][0] + h[5]*x[i][1] + h[4]*x[i][2] + boxlo[0]; + x[i][1] = h[1]*x[i][1] + h[3]*x[i][2] + boxlo[1]; + x[i][2] = h[2]*x[i][2] + boxlo[2]; + } +} + +/* ---------------------------------------------------------------------- + convert box coords to triclinic 0-1 lamda coords for all N atoms + lamda = H^-1 (x - x0) +------------------------------------------------------------------------- */ + +void Domain::x2lamda(int n) +{ + double delta[3]; + double **x = atom->x; + + for (int i = 0; i < n; i++) { + delta[0] = x[i][0] - boxlo[0]; + delta[1] = x[i][1] - boxlo[1]; + delta[2] = x[i][2] - boxlo[2]; + + x[i][0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; + x[i][1] = h_inv[1]*delta[1] + h_inv[3]*delta[2]; + x[i][2] = h_inv[2]*delta[2]; + } +} + +/* ---------------------------------------------------------------------- + convert triclinic 0-1 lamda coords to box coords for one atom + x = H lamda + x0; + lamda and x can point to same 3-vector +------------------------------------------------------------------------- */ + +void Domain::lamda2x(double *lamda, double *x) +{ + x[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2] + boxlo[0]; + x[1] = h[1]*lamda[1] + h[3]*lamda[2] + boxlo[1]; + x[2] = h[2]*lamda[2] + boxlo[2]; +} + +/* ---------------------------------------------------------------------- + convert box coords to triclinic 0-1 lamda coords for one atom + lamda = H^-1 (x - x0) + x and lamda can point to same 3-vector +------------------------------------------------------------------------- */ + +void Domain::x2lamda(double *x, double *lamda) +{ + double delta[3]; + delta[0] = x[0] - boxlo[0]; + delta[1] = x[1] - boxlo[1]; + delta[2] = x[2] - boxlo[2]; + + lamda[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2]; + lamda[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2]; + lamda[2] = h_inv[2]*delta[2]; +} diff --git a/src/domain.h b/src/domain.h index 11f88ec6b4..16a6949271 100644 --- a/src/domain.h +++ b/src/domain.h @@ -20,61 +20,91 @@ namespace LAMMPS_NS { class Domain : protected Pointers { public: - int box_exist; // 0 = not yet created, 1 = exists + int box_exist; // 0 = not yet created, 1 = exists - int nonperiodic; // 0 = periodic in all 3 dims - // 1 = periodic or fixed in all 6 - // 2 = shrink-wrap in any of 6 - int xperiodic,yperiodic,zperiodic; // 0 = not periodic, 1 = periodic - int boundary[3][2]; // settings for 6 boundaries - // 0 = periodic - // 1 = fixed non-periodic - // 2 = shrink-wrap non-periodic - // 3 = shrink-wrap non-per w/ min + int nonperiodic; // 0 = periodic in all 3 dims + // 1 = periodic or fixed in all 6 + // 2 = shrink-wrap in any of 6 + int xperiodic,yperiodic,zperiodic; // 0 = non-periodic, 1 = periodic + int periodicity[3]; // xyz periodicity as array - double minxlo,minxhi; // minimum size of global box - double minylo,minyhi; // when shrink-wrapping - double minzlo,minzhi; + int boundary[3][2]; // settings for 6 boundaries + // 0 = periodic + // 1 = fixed non-periodic + // 2 = shrink-wrap non-periodic + // 3 = shrink-wrap non-per w/ min - double boxxlo,boxxhi; // global box boundaries - double boxylo,boxyhi; + int triclinic; // 0 = orthog box, 1 = triclinic + + // orthogonal box + double xprd,yprd,zprd; // global box dimensions + double xprd_half,yprd_half,zprd_half; // half dimensions + double prd[3]; // array form of dimensions + + // triclinic box + // xprd,half,prd = same + // (as if untilted) + double prd_lamda[3]; // lamda box = (1,1,1) + + double boxxlo,boxxhi; // orthogonal box + double boxylo,boxyhi; // global box bounds double boxzlo,boxzhi; + double boxlo[3],boxhi[3]; // array form of box bounds - double xprd,yprd,zprd; // global box size - double xprd_half,yprd_half,zprd_half; + // triclinic box + // boxxlo/hi,boxlo/hi = same + // (as if untilted) + double boxlo_lamda[3],boxhi_lamda[3]; // lamda box = (0,1) + double boxlo_bound[3],boxhi_bound[3]; // bounding box of tilted domain - double subxlo,subxhi; // sub-box boudaries on this proc - double subylo,subyhi; - double subzlo,subzhi; + // orthogonal box + double minxlo,minxhi; // minimum size of global box + double minylo,minyhi; // when shrink-wrapping + double minzlo,minzhi; // no shrink-wrapping for triclinic - double boxlo[3],boxhi[3]; // global box bounds as arrays - double prd[3]; // global box size as array - double sublo[3],subhi[3]; // sub-box bounds as arrays - int periodicity[3]; // xyz periodic as array + // orthogonal box + double sublo[3],subhi[3]; // sub-box bounds on this proc - int box_change; // 1 if box bounds ever change, 0 if fixed + // triclinic box + // sublo = lower-left corner with tilt + // subhi = upper-right as if untilted + double sublo_lamda[3],subhi_lamda[3]; // bounds of subbox in lamda + double sublo_bound[3],subhi_bound[3]; // bounding box of tilted subdomain - class Lattice *lattice; // user-defined lattice + // triclinic box + double xy,xz,yz; // 3 tilt factors + double h[6],h_inv[6]; // shape matrix in Voigt notation - int nregion; // # of defined Regions - int maxregion; // max # list can hold - class Region **regions; // list of defined Regions + int box_change; // 1 if box bounds ever change, 0 if fixed + + class Lattice *lattice; // user-defined lattice + + int nregion; // # of defined Regions + int maxregion; // max # list can hold + class Region **regions; // list of defined Regions Domain(class LAMMPS *); ~Domain(); void init(); void set_initial_box(); void set_global_box(); + void set_lamda_box(); void set_local_box(); void reset_box(); void pbc(); - void remap(double &, double &, double &, int &); - void unmap(double &, double &, double &, int); - void minimum_image(double *, double *, double *); + void remap(double *, int &); + void unmap(double *, int); + void minimum_image(double &, double &, double &); void minimum_image(double *); void set_lattice(int, char **); void add_region(int, char **); void set_boundary(int, char **); + void print_box(char *); + + void lamda2x(int); + void x2lamda(int); + void lamda2x(double *, double *); + void x2lamda(double *, double *); }; } diff --git a/src/dump.cpp b/src/dump.cpp index 8b269f69c1..b8c7bdf51e 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -18,6 +18,7 @@ #include "dump.h" #include "atom.h" #include "update.h" +#include "domain.h" #include "group.h" #include "output.h" #include "memory.h" @@ -118,6 +119,24 @@ void Dump::write() if (multifile) openfile(); + // simulation box bounds + + if (domain->triclinic == 0) { + boxxlo = domain->boxxlo; + boxxhi = domain->boxxhi; + boxylo = domain->boxylo; + boxyhi = domain->boxyhi; + boxzlo = domain->boxzlo; + boxzhi = domain->boxzhi; + } else { + boxxlo = domain->boxlo_bound[0]; + boxxhi = domain->boxhi_bound[0]; + boxylo = domain->boxlo_bound[1]; + boxyhi = domain->boxhi_bound[1]; + boxzlo = domain->boxlo_bound[2]; + boxzhi = domain->boxhi_bound[2]; + } + // nme = # of dump lines this proc will contribute to dump // ntotal = total # of dump lines // nmax = max # of dump lines on any proc diff --git a/src/dump.h b/src/dump.h index 906fc69368..4fa6925c8a 100644 --- a/src/dump.h +++ b/src/dump.h @@ -51,6 +51,10 @@ class Dump : protected Pointers { virtual int memory_usage(); protected: + double boxxlo,boxxhi; // local copies of domain values + double boxylo,boxyhi; // adjusted to be bounding box for triclinic + double boxzlo,boxzhi; + virtual void openfile(); virtual int modify_param(int, char **) {return 0;} diff --git a/src/dump_atom.cpp b/src/dump_atom.cpp index ea7b25eed9..3640feec18 100644 --- a/src/dump_atom.cpp +++ b/src/dump_atom.cpp @@ -40,6 +40,9 @@ DumpAtom::DumpAtom(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) void DumpAtom::init() { + if (scale_flag && domain->triclinic) + error->all("Cannot dump scaled coords with triclinic box"); + if (image_flag == 0) size_one = 5; else size_one = 8; @@ -142,12 +145,12 @@ void DumpAtom::header_binary(int ndump) { fwrite(&update->ntimestep,sizeof(int),1,fp); fwrite(&ndump,sizeof(int),1,fp); - fwrite(&domain->boxxlo,sizeof(double),1,fp); - fwrite(&domain->boxxhi,sizeof(double),1,fp); - fwrite(&domain->boxylo,sizeof(double),1,fp); - fwrite(&domain->boxyhi,sizeof(double),1,fp); - fwrite(&domain->boxzlo,sizeof(double),1,fp); - fwrite(&domain->boxzhi,sizeof(double),1,fp); + fwrite(&boxxlo,sizeof(double),1,fp); + fwrite(&boxxhi,sizeof(double),1,fp); + fwrite(&boxylo,sizeof(double),1,fp); + fwrite(&boxyhi,sizeof(double),1,fp); + fwrite(&boxzlo,sizeof(double),1,fp); + fwrite(&boxzhi,sizeof(double),1,fp); fwrite(&size_one,sizeof(int),1,fp); if (multiproc) { int one = 1; @@ -164,9 +167,9 @@ void DumpAtom::header_item(int ndump) fprintf(fp,"ITEM: NUMBER OF ATOMS\n"); fprintf(fp,"%d\n",ndump); fprintf(fp,"ITEM: BOX BOUNDS\n"); - fprintf(fp,"%g %g\n",domain->boxxlo,domain->boxxhi); - fprintf(fp,"%g %g\n",domain->boxylo,domain->boxyhi); - fprintf(fp,"%g %g\n",domain->boxzlo,domain->boxzhi); + fprintf(fp,"%g %g\n",boxxlo,boxxhi); + fprintf(fp,"%g %g\n",boxylo,boxyhi); + fprintf(fp,"%g %g\n",boxzlo,boxzhi); fprintf(fp,"ITEM: ATOMS\n"); } @@ -181,9 +184,6 @@ int DumpAtom::pack_scale_image() double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; double invxprd = 1.0/domain->xprd; double invyprd = 1.0/domain->yprd; double invzprd = 1.0/domain->zprd; @@ -214,9 +214,6 @@ int DumpAtom::pack_scale_noimage() double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; double invxprd = 1.0/domain->xprd; double invyprd = 1.0/domain->yprd; double invzprd = 1.0/domain->zprd; diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 4ea7e77b86..e732af5ac4 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -195,12 +195,12 @@ void DumpCustom::header_binary(int ndump) { fwrite(&update->ntimestep,sizeof(int),1,fp); fwrite(&ndump,sizeof(int),1,fp); - fwrite(&domain->boxxlo,sizeof(double),1,fp); - fwrite(&domain->boxxhi,sizeof(double),1,fp); - fwrite(&domain->boxylo,sizeof(double),1,fp); - fwrite(&domain->boxyhi,sizeof(double),1,fp); - fwrite(&domain->boxzlo,sizeof(double),1,fp); - fwrite(&domain->boxzhi,sizeof(double),1,fp); + fwrite(&boxxlo,sizeof(double),1,fp); + fwrite(&boxxhi,sizeof(double),1,fp); + fwrite(&boxylo,sizeof(double),1,fp); + fwrite(&boxyhi,sizeof(double),1,fp); + fwrite(&boxzlo,sizeof(double),1,fp); + fwrite(&boxzhi,sizeof(double),1,fp); fwrite(&size_one,sizeof(int),1,fp); if (multiproc) { int one = 1; @@ -217,9 +217,9 @@ void DumpCustom::header_item(int ndump) fprintf(fp,"ITEM: NUMBER OF ATOMS\n"); fprintf(fp,"%d\n",ndump); fprintf(fp,"ITEM: BOX BOUNDS\n"); - fprintf(fp,"%g %g\n",domain->boxxlo,domain->boxxhi); - fprintf(fp,"%g %g\n",domain->boxylo,domain->boxyhi); - fprintf(fp,"%g %g\n",domain->boxzlo,domain->boxzhi); + fprintf(fp,"%g %g\n",boxxlo,boxxhi); + fprintf(fp,"%g %g\n",boxylo,boxyhi); + fprintf(fp,"%g %g\n",boxzlo,boxzhi); fprintf(fp,"ITEM: ATOMS\n"); } @@ -550,6 +550,7 @@ void DumpCustom::parse_fields(int narg, char **arg) int i; for (int iarg = 5; iarg < narg; iarg++) { i = iarg-5; + if (strcmp(arg[iarg],"tag") == 0) { pack_choice[i] = &DumpCustom::pack_tag; vtype[i] = INT; @@ -572,12 +573,18 @@ void DumpCustom::parse_fields(int narg, char **arg) pack_choice[i] = &DumpCustom::pack_z; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"xs") == 0) { + if (domain->triclinic) + error->all("Cannot dump scaled coords with triclinic box"); pack_choice[i] = &DumpCustom::pack_xs; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"ys") == 0) { + if (domain->triclinic) + error->all("Cannot dump scaled coords with triclinic box"); pack_choice[i] = &DumpCustom::pack_ys; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"zs") == 0) { + if (domain->triclinic) + error->all("Cannot dump scaled coords with triclinic box"); pack_choice[i] = &DumpCustom::pack_zs; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"xu") == 0) { @@ -844,6 +851,13 @@ int DumpCustom::modify_param(int narg, char **arg) memory->srealloc(thresh_value,(nthresh+1)*sizeof(double), "dump:thresh_value"); + // error check for triclinic box + + if (domain->triclinic && + (strcmp(arg[1],"xs") == 0 || strcmp(arg[1],"ys") == 0 || + strcmp(arg[1],"zs") == 0)) + error->all("Cannot dump scaled coords with triclinic box"); + // set keyword type of threshhold // customize by adding to if statement diff --git a/src/dump_dcd.cpp b/src/dump_dcd.cpp index ef9b73ff40..92c5a3382c 100644 --- a/src/dump_dcd.cpp +++ b/src/dump_dcd.cpp @@ -15,6 +15,7 @@ Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U) ------------------------------------------------------------------------- */ +#include "math.h" #include "inttypes.h" #include "stdio.h" #include "time.h" @@ -133,14 +134,26 @@ void DumpDCD::write_header(int n) nframes = 0; } - // dim[134] = angle cosines of periodic box + // dim[] = size and angle cosines of orthogonal or triclinic box + // dim[1] = alpha = cosine of angle between b and c + // dim[3] = beta = cosine of angle between a and c + // dim[4] = gamma = cosine of angle between a and b // 48 = 6 doubles double dim[6]; - dim[0] = domain->boxxhi - domain->boxxlo; - dim[2] = domain->boxyhi - domain->boxylo; - dim[5] = domain->boxzhi - domain->boxzlo; - dim[1] = dim[3] = dim[4] = 0.0; + dim[0] = domain->xprd; + dim[2] = domain->yprd; + dim[5] = domain->zprd; + if (domain->triclinic == 0) dim[1] = dim[3] = dim[4] = 0.0; + else { + double *h = domain->h; + double alen = h[0]; + double blen = sqrt(h[5]*h[5] + h[1]*h[1]); + double clen = sqrt(h[4]*h[4] + h[3]*h[3] + h[2]*h[2]); + dim[1] = (h[5]*h[4] + h[1]*h[3]) / blen/clen; + dim[3] = (h[0]*h[4]) / alen/clen; + dim[4] = (h[0]*h[5]) / alen/blen; + } uint32_t out_integer = 48; fwrite(&out_integer,sizeof(uint32_t),1,fp); diff --git a/src/fix.cpp b/src/fix.cpp index 9b0b2cc3c1..be6813c75f 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -36,6 +36,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) restart_global = 0; restart_peratom = 0; force_reneighbor = 0; + box_change = 0; thermo_energy = 0; pressure_every = 0; diff --git a/src/fix.h b/src/fix.h index 176cc11ec8..97b18bb578 100644 --- a/src/fix.h +++ b/src/fix.h @@ -26,6 +26,7 @@ class Fix : protected Pointers { int restart_global; // 1 if Fix saves global state, 0 if not int restart_peratom; // 1 if Fix saves peratom state, 0 if not int force_reneighbor; // 1 if Fix forces reneighboring, 0 if not + int box_change; // 1 if Fix changes box size, 0 if not int next_reneighbor; // next timestep to force a reneighboring int thermo_energy; // 1 if ThEng enabled via fix_modify, 0 if not int nevery; // how often to call an end_of_step fix @@ -81,7 +82,7 @@ class Fix : protected Pointers { virtual int min_dof() {return 0;} virtual void min_xinitial(double *) {} - virtual int pack_comm(int, int *, double *, int *) {return 0;} + virtual int pack_comm(int, int *, double *, int, double *) {return 0;} virtual void unpack_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index 69523e8789..f2c90eae07 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -147,6 +147,15 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) : if (delta <= 0.0) error->all("Illegal fix ave/spatial command"); invdelta = 1.0/delta; + if (domain->triclinic) { + if (dim == 0 && (domain->xy != 0.0 || domain->xz != 0.0)) + error->all("Cannot (yet) use fix ave/spatial with triclinic box"); + if (dim == 1 && (domain->xy != 0.0 || domain->yz != 0.0)) + error->all("Cannot (yet) use fix ave/spatial with triclinic box"); + if (dim == 2 && (domain->xz != 0.0 || domain->yz != 0.0)) + error->all("Cannot (yet) use fix ave/spatial with triclinic box"); + } + nvalues = 1; if (which == COMPUTE) { int icompute = modify->find_compute(id_compute); diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp index f24ee98da1..4420be115e 100644 --- a/src/fix_deposit.cpp +++ b/src/fix_deposit.cpp @@ -150,7 +150,8 @@ int FixDeposit::setmask() void FixDeposit::pre_exchange() { int flag,flagall; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double coord[3],lamda[3],delx,dely,delz,rsq; + double *newcoord; // just return if should not be called on this timestep @@ -161,6 +162,15 @@ void FixDeposit::pre_exchange() double offset = 0.0; if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate; + double *sublo,*subhi; + if (domain->triclinic == 0) { + sublo = domain->sublo; + subhi = domain->subhi; + } else { + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + // attempt an insertion until successful int success = 0; @@ -170,19 +180,19 @@ void FixDeposit::pre_exchange() // choose random position for new atom within region - xtmp = xlo + random->uniform() * (xhi-xlo); - ytmp = ylo + random->uniform() * (yhi-ylo); - ztmp = zlo + random->uniform() * (zhi-zlo); - while (domain->regions[iregion]->match(xtmp,ytmp,ztmp) == 0) { - xtmp = xlo + random->uniform() * (xhi-xlo); - ytmp = ylo + random->uniform() * (yhi-ylo); - ztmp = zlo + random->uniform() * (zhi-zlo); + coord[0] = xlo + random->uniform() * (xhi-xlo); + coord[1] = ylo + random->uniform() * (yhi-ylo); + coord[2] = zlo + random->uniform() * (zhi-zlo); + while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + coord[0] = xlo + random->uniform() * (xhi-xlo); + coord[1] = ylo + random->uniform() * (yhi-ylo); + coord[2] = zlo + random->uniform() * (zhi-zlo); } // adjust vertical coord by offset - if (force->dimension == 2) ytmp += offset; - else ztmp += offset; + if (force->dimension == 2) coord[1] += offset; + else coord[2] += offset; // if global, reset vertical coord to be lo-hi above highest atom // if local, reset vertical coord to be lo-hi above highest "nearby" atom @@ -204,10 +214,10 @@ void FixDeposit::pre_exchange() int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { if (localflag) { - delx = xtmp - x[i][0]; - dely = ytmp - x[i][1]; + delx = coord[0] - x[i][0]; + dely = coord[1] - x[i][1]; delz = 0.0; - domain->minimum_image(&delx,&dely,&delz); + domain->minimum_image(delx,dely,delz); if (force->dimension == 2) rsq = delx*delx; else rsq = delx*delx + dely*dely; if (rsq > deltasq) continue; @@ -217,12 +227,12 @@ void FixDeposit::pre_exchange() MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world); if (force->dimension == 2) - ytmp = maxall + lo + random->uniform()*(hi-lo); + coord[1] = maxall + lo + random->uniform()*(hi-lo); else - ztmp = maxall + lo + random->uniform()*(hi-lo); + coord[2] = maxall + lo + random->uniform()*(hi-lo); } - // now have final xtmp,ytmp,ztmp + // now have final coord // if distance to any atom is less than near, try again double **x = atom->x; @@ -230,10 +240,10 @@ void FixDeposit::pre_exchange() flag = 0; for (int i = 0; i < nlocal; i++) { - delx = xtmp - x[i][0]; - dely = ytmp - x[i][1]; - delz = ztmp - x[i][2]; - domain->minimum_image(&delx,&dely,&delz); + delx = coord[0] - x[i][0]; + dely = coord[1] - x[i][1]; + delz = coord[2] - x[i][2]; + domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; if (rsq < nearsq) flag = 1; } @@ -251,21 +261,26 @@ void FixDeposit::pre_exchange() // if so, add to my list via create_atom() // initialize info about the atoms // set group mask to "all" plus fix group + + if (domain->triclinic) { + domain->x2lamda(coord,lamda); + newcoord = lamda; + } else newcoord = coord; flag = 0; - if (xtmp >= domain->subxlo && xtmp < domain->subxhi && - ytmp >= domain->subylo && ytmp < domain->subyhi && - ztmp >= domain->subzlo && ztmp < domain->subzhi) flag = 1; - else if (force->dimension == 3 && ztmp >= domain->boxzhi && + if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] && + newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1; + else if (force->dimension == 3 && newcoord[2] >= domain->boxzhi && comm->myloc[2] == comm->procgrid[2]-1 && - xtmp >= domain->subxlo && xtmp < domain->subxhi && - ytmp >= domain->subylo && ytmp < domain->subyhi) flag = 1; - else if (force->dimension == 2 && ytmp >= domain->boxyhi && + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; + else if (force->dimension == 2 && newcoord[1] >= domain->boxyhi && comm->myloc[1] == comm->procgrid[1]-1 && - xtmp >= domain->subxlo && xtmp < domain->subxhi) flag = 1; + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; if (flag) { - atom->avec->create_atom(ntype,xtmp,ytmp,ztmp,0); + atom->avec->create_atom(ntype,coord,0); int m = atom->nlocal - 1; atom->type[m] = ntype; atom->mask[m] = 1 | groupbit; diff --git a/src/fix_drag.cpp b/src/fix_drag.cpp index d3220bae8d..0127f27db6 100644 --- a/src/fix_drag.cpp +++ b/src/fix_drag.cpp @@ -96,7 +96,7 @@ void FixDrag::post_force(int vflag) if (!xflag) dx = 0.0; if (!yflag) dy = 0.0; if (!zflag) dz = 0.0; - domain->minimum_image(&dx,&dy,&dz); + domain->minimum_image(dx,dy,dz); r = sqrt(dx*dx + dy*dy + dz*dz); if (r > delta) { prefactor = f_mag/r; diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp index a760575075..b2867c4c99 100644 --- a/src/fix_nph.cpp +++ b/src/fix_nph.cpp @@ -45,6 +45,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : restart_global = 1; pressure_every = 1; + box_change = 1; double p_period[3]; if (strcmp(arg[3],"xyz") == 0) { @@ -119,14 +120,15 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : } else error->all("Illegal fix nph command"); } - // check for periodicity in controlled dimensions + // error checks + if (domain->triclinic) error->all("Cannot use fix nph with triclinic box"); if (p_flag[0] && domain->xperiodic == 0) - error->all("Cannot fix nph on a non-periodic dimension"); + error->all("Cannot use fix nph on a non-periodic dimension"); if (p_flag[1] && domain->yperiodic == 0) - error->all("Cannot fix nph on a non-periodic dimension"); + error->all("Cannot use fix nph on a non-periodic dimension"); if (p_flag[2] && domain->zperiodic == 0) - error->all("Cannot fix nph on a non-periodic dimension"); + error->all("Cannot use fix nph on a non-periodic dimension"); // convert input periods to frequencies diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index cb73f8e296..3886234273 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -45,6 +45,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : restart_global = 1; pressure_every = 1; + box_change = 1; t_start = atof(arg[3]); t_stop = atof(arg[4]); @@ -126,14 +127,15 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : } else error->all("Illegal fix npt command"); } - // check for periodicity in controlled dimensions + // error checks + if (domain->triclinic) error->all("Cannot use fix npt with triclinic box"); if (p_flag[0] && domain->xperiodic == 0) - error->all("Cannot fix npt on a non-periodic dimension"); + error->all("Cannot use fix npt on a non-periodic dimension"); if (p_flag[1] && domain->yperiodic == 0) - error->all("Cannot fix npt on a non-periodic dimension"); + error->all("Cannot use fix npt on a non-periodic dimension"); if (p_flag[2] && domain->zperiodic == 0) - error->all("Cannot fix npt on a non-periodic dimension"); + error->all("Cannot use fix npt on a non-periodic dimension"); // convert input periods to frequencies diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp index 804e35dacb..d6ac434d2c 100644 --- a/src/fix_orient_fcc.cpp +++ b/src/fix_orient_fcc.cpp @@ -419,7 +419,8 @@ double FixOrientFCC::thermo(int n) /* ---------------------------------------------------------------------- */ -int FixOrientFCC::pack_comm(int n, int *list, double *buf, int *pbc_flags) +int FixOrientFCC::pack_comm(int n, int *list, double *buf, + int pbc_flag, double *pair_dist) { int i,j,k,id,num; int *tag = atom->tag; diff --git a/src/fix_orient_fcc.h b/src/fix_orient_fcc.h index 1ec09a6833..cdfc05673e 100644 --- a/src/fix_orient_fcc.h +++ b/src/fix_orient_fcc.h @@ -45,7 +45,7 @@ class FixOrientFCC : public Fix { void post_force(int); void post_force_respa(int, int, int); double thermo(int); - int pack_comm(int, int *, double *, int *); + int pack_comm(int, int *, double *, int, double *); void unpack_comm(int, int, double *); int memory_usage(); diff --git a/src/fix_recenter.cpp b/src/fix_recenter.cpp index 5f22435243..f870a00eab 100644 --- a/src/fix_recenter.cpp +++ b/src/fix_recenter.cpp @@ -135,22 +135,31 @@ void FixRecenter::init() void FixRecenter::initial_integrate() { // target COM + // bounding box around domain works for both orthog and triclinic double xtarget,ytarget,ztarget; + double *bboxlo,*bboxhi; + + if (scaleflag == 2) { + if (domain->triclinic == 0) { + bboxlo = domain->boxlo; + bboxhi = domain->boxhi; + } else { + bboxlo = domain->boxlo_bound; + bboxhi = domain->boxhi_bound; + } + } if (xinitflag) xtarget = xinit; - else if (scaleflag == 2) - xtarget = domain->boxxlo + xcom*(domain->boxxhi - domain->boxxlo); + else if (scaleflag == 2) xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]); else xtarget = xcom; if (yinitflag) ytarget = yinit; - else if (scaleflag == 2) - ytarget = domain->boxylo + ycom*(domain->boxyhi - domain->boxylo); + else if (scaleflag == 2) ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]); else ytarget = ycom; if (zinitflag) ztarget = zinit; - else if (scaleflag == 2) - ztarget = domain->boxzlo + zcom*(domain->boxzhi - domain->boxzlo); + else if (scaleflag == 2) ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]); else ztarget = zcom; // current COM diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp index e466e566ad..263d742d84 100644 --- a/src/fix_shake.cpp +++ b/src/fix_shake.cpp @@ -1911,7 +1911,7 @@ void FixShake::stats() delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; delz = x[iatom][2] - x[jatom][2]; - domain->minimum_image(&delx,&dely,&delz); + domain->minimum_image(delx,dely,delz); r = sqrt(delx*delx + dely*dely + delz*delz); m = shake_type[i][j-1]; @@ -1931,19 +1931,19 @@ void FixShake::stats() delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; delz = x[iatom][2] - x[jatom][2]; - domain->minimum_image(&delx,&dely,&delz); + domain->minimum_image(delx,dely,delz); r1 = sqrt(delx*delx + dely*dely + delz*delz); delx = x[iatom][0] - x[katom][0]; dely = x[iatom][1] - x[katom][1]; delz = x[iatom][2] - x[katom][2]; - domain->minimum_image(&delx,&dely,&delz); + domain->minimum_image(delx,dely,delz); r2 = sqrt(delx*delx + dely*dely + delz*delz); delx = x[jatom][0] - x[katom][0]; dely = x[jatom][1] - x[katom][1]; delz = x[jatom][2] - x[katom][2]; - domain->minimum_image(&delx,&dely,&delz); + domain->minimum_image(delx,dely,delz); r3 = sqrt(delx*delx + dely*dely + delz*delz); angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2)); @@ -2252,12 +2252,13 @@ void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop) /* ---------------------------------------------------------------------- */ -int FixShake::pack_comm(int n, int *list, double *buf, int *pbc_flags) +int FixShake::pack_comm(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = xshake[j][0]; @@ -2267,9 +2268,9 @@ int FixShake::pack_comm(int n, int *list, double *buf, int *pbc_flags) } else { for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = xshake[j][0] + pbc_flags[1]*domain->xprd; - buf[m++] = xshake[j][1] + pbc_flags[2]*domain->yprd; - buf[m++] = xshake[j][2] + pbc_flags[3]*domain->zprd; + buf[m++] = xshake[j][0] + pbc_dist[0]; + buf[m++] = xshake[j][1] + pbc_dist[1]; + buf[m++] = xshake[j][2] + pbc_dist[2]; } } return 3; diff --git a/src/fix_shake.h b/src/fix_shake.h index 47cdf7be29..09e78dffc2 100644 --- a/src/fix_shake.h +++ b/src/fix_shake.h @@ -34,7 +34,7 @@ class FixShake : public Fix { void copy_arrays(int, int); int pack_exchange(int, double *); int unpack_exchange(int, double *); - int pack_comm(int, int *, double *, int *); + int pack_comm(int, int *, double *, int, double *); void unpack_comm(int, int, double *); int dof(int); diff --git a/src/fix_uniaxial.cpp b/src/fix_uniaxial.cpp index a63a136b17..01496ad0e8 100644 --- a/src/fix_uniaxial.cpp +++ b/src/fix_uniaxial.cpp @@ -47,8 +47,11 @@ FixUniaxial::FixUniaxial(LAMMPS *lmp, int narg, char **arg) : lambda_final = atof(arg[5]); if (lambda_final <= 0) error->all("Illegal fix uniaxial command"); + if (domain->nonperiodic) - error->all("Cannot fix uniaxial on non-periodic system"); + error->all("Cannot use fix uniaxial on non-periodic system"); + if (domain->triclinic) + error->all("Cannot use fix uniaxial with triclinic box"); nrigid = 0; rfix = NULL; @@ -204,4 +207,3 @@ void FixUniaxial::end_of_step() if (kspace_flag) force->kspace->setup(); } - diff --git a/src/fix_volume_rescale.cpp b/src/fix_volume_rescale.cpp index 78ece5b58b..413172aa21 100644 --- a/src/fix_volume_rescale.cpp +++ b/src/fix_volume_rescale.cpp @@ -31,6 +31,9 @@ FixVolRescale::FixVolRescale(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 7) error->all("Illegal fix volume/rescale command"); + + box_change = 1; + nevery = atoi(arg[3]); if (nevery <= 0) error->all("Illegal fix volume/rescale command"); @@ -64,6 +67,9 @@ FixVolRescale::FixVolRescale(LAMMPS *lmp, int narg, char **arg) : } else error->all("Illegal fix volume/rescale command"); } + if (domain->triclinic) + error->all("Cannot use fix volume/rescale with triclinic box"); + nrigid = 0; rfix = NULL; } diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp index 6ef6e15191..a0ff286373 100644 --- a/src/fix_wall_lj126.cpp +++ b/src/fix_wall_lj126.cpp @@ -20,6 +20,7 @@ #include "string.h" #include "fix_wall_lj126.h" #include "atom.h" +#include "domain.h" #include "update.h" #include "output.h" #include "respa.h" @@ -67,6 +68,13 @@ FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) : double r2inv = 1.0/(cutoff*cutoff); double r6inv = r2inv*r2inv*r2inv; offset = r6inv*(coeff3*r6inv - coeff4); + + if (dim == 0 && domain->xperiodic) + error->all("Cannot use wall in periodic dimension"); + if (dim == 1 && domain->yperiodic) + error->all("Cannot use wall in periodic dimension"); + if (dim == 2 && domain->zperiodic) + error->all("Cannot use wall in periodic dimension"); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp index fa9099ca00..7d40af983d 100644 --- a/src/fix_wall_lj93.cpp +++ b/src/fix_wall_lj93.cpp @@ -16,6 +16,7 @@ #include "string.h" #include "fix_wall_lj93.h" #include "atom.h" +#include "domain.h" #include "update.h" #include "output.h" #include "respa.h" @@ -64,6 +65,13 @@ FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) : double r2inv = rinv*rinv; double r4inv = r2inv*r2inv; offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv; + + if (dim == 0 && domain->xperiodic) + error->all("Cannot use wall in periodic dimension"); + if (dim == 1 && domain->yperiodic) + error->all("Cannot use wall in periodic dimension"); + if (dim == 2 && domain->zperiodic) + error->all("Cannot use wall in periodic dimension"); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index d193400e99..da8ee9bb5c 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -38,6 +38,13 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg],"zhi") == 0) zhiflag = 1; else error->all("Illegal fix wall/reflect command"); } + + if ((xloflag || xhiflag) && domain->xperiodic) + error->all("Cannot use wall in periodic dimension"); + if ((yloflag || yhiflag) && domain->yperiodic) + error->all("Cannot use wall in periodic dimension"); + if ((zloflag || zhiflag) && domain->zperiodic) + error->all("Cannot use wall in periodic dimension"); } /* ---------------------------------------------------------------------- */ diff --git a/src/lattice.cpp b/src/lattice.cpp index 8890e82a84..70abc11280 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -27,7 +27,7 @@ using namespace LAMMPS_NS; #define MAX(A,B) ((A) > (B)) ? (A) : (B) #define BIG 1.0e30 -enum{NONE,SC,BCC,FCC,DIAMOND,SQ,SQ2,HEX,USER}; +enum{NONE,SC,BCC,FCC,DIAMOND,SQ,SQ2,HEX,CUSTOM}; /* ---------------------------------------------------------------------- */ @@ -45,7 +45,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) else if (strcmp(arg[0],"sq") == 0) style = SQ; else if (strcmp(arg[0],"sq2") == 0) style = SQ2; else if (strcmp(arg[0],"hex") == 0) style = HEX; - else if (strcmp(arg[0],"user") == 0) style = USER; + else if (strcmp(arg[0],"custom") == 0) style = CUSTOM; else error->all("Illegal lattice command"); if (style == NONE) { @@ -54,7 +54,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) } // check that lattice matches dimension - // style USER can be either 2d or 3d + // style CUSTOM can be either 2d or 3d int dimension = force->dimension; if (dimension == 2) { @@ -74,7 +74,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) // set basis atoms for each style // x,y,z = fractional coords within unit cell - // style USER will be defined by optional args + // style CUSTOM will be defined by optional args nbasis = 0; basis = NULL; @@ -116,6 +116,8 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) orienty[0] = 0; orienty[1] = 1; orienty[2] = 0; orientz[0] = 0; orientz[1] = 0; orientz[2] = 1; + int spaceflag = 0; + a1[0] = 1.0; a1[1] = 0.0; a1[2] = 0.0; a2[0] = 0.0; a2[1] = 1.0; a2[2] = 0.0; a3[0] = 0.0; a3[1] = 0.0; a3[2] = 1.0; @@ -153,26 +155,34 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) orient[2] = atoi(arg[iarg+4]); iarg += 5; + } else if (strcmp(arg[iarg],"spacings") == 0) { + if (iarg+4 > narg) error->all("Illegal lattice command"); + spaceflag = 1; + xlattice = atoi(arg[iarg+1]); + ylattice = atoi(arg[iarg+2]); + zlattice = atoi(arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg],"a1") == 0) { if (iarg+4 > narg) error->all("Illegal lattice command"); - if (style != USER) - error->all("Invalid option in lattice command for non-user style"); + if (style != CUSTOM) + error->all("Invalid option in lattice command for non-custom style"); a1[0] = atof(arg[iarg+1]); a1[1] = atof(arg[iarg+2]); a1[2] = atof(arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"a2") == 0) { if (iarg+4 > narg) error->all("Illegal lattice command"); - if (style != USER) - error->all("Invalid option in lattice command for non-user style"); + if (style != CUSTOM) + error->all("Invalid option in lattice command for non-custom style"); a2[0] = atof(arg[iarg+1]); a2[1] = atof(arg[iarg+2]); a2[2] = atof(arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"a3") == 0) { if (iarg+4 > narg) error->all("Illegal lattice command"); - if (style != USER) - error->all("Invalid option in lattice command for non-user style"); + if (style != CUSTOM) + error->all("Invalid option in lattice command for non-custom style"); a3[0] = atof(arg[iarg+1]); a3[1] = atof(arg[iarg+2]); a3[2] = atof(arg[iarg+3]); @@ -180,8 +190,8 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) } else if (strcmp(arg[iarg],"basis") == 0) { if (iarg+4 > narg) error->all("Illegal lattice command"); - if (style != USER) - error->all("Invalid option in lattice command for non-user style"); + if (style != CUSTOM) + error->all("Invalid option in lattice command for non-custom style"); double x = atof(arg[iarg+1]); double y = atof(arg[iarg+2]); double z = atof(arg[iarg+3]); @@ -212,6 +222,11 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) error->all("Lattice settings are not compatible with 2d simulation"); } + if (spaceflag) { + if (xlattice <= 0.0 || ylattice <= 0.0 || zlattice <= 0.0) + error->all("Lattice spacings are invalid"); + } + // reset scale for LJ units (input scale is rho*) // scale = (Nbasis/(Vprimitive * rho*)) ^ (1/dim) @@ -232,23 +247,30 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) // set xlattice,ylattice,zlattice to 0.0 initially // since bbox uses them to shift origin (irrelevant for this computation) - double xmin,ymin,zmin,xmax,ymax,zmax; - xmin = ymin = zmin = BIG; - xmax = ymax = zmax = -BIG; - xlattice = ylattice = zlattice = 0.0; + if (spaceflag == 0) { + double xmin,ymin,zmin,xmax,ymax,zmax; + xmin = ymin = zmin = BIG; + xmax = ymax = zmax = -BIG; + xlattice = ylattice = zlattice = 0.0; - bbox(0,0.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); - bbox(0,1.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); - bbox(0,0.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); - bbox(0,1.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); - bbox(0,0.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); - bbox(0,1.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); - bbox(0,0.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); - bbox(0,1.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,0.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,1.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,0.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,1.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,0.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,1.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,0.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); + bbox(0,1.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax); - xlattice = xmax - xmin; - ylattice = ymax - ymin; - zlattice = zmax - zmin; + xlattice = xmax - xmin; + ylattice = ymax - ymin; + zlattice = zmax - zmin; + + } else { + xlattice *= scale; + ylattice *= scale; + zlattice *= scale; + } // print lattice spacings diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp index 91a12246ca..a8ae6556b3 100644 --- a/src/neigh_full.cpp +++ b/src/neigh_full.cpp @@ -131,16 +131,9 @@ void Neighbor::full_bin() // only skip i = j for (k = 0; k < nstencil_full; k++) { - j = binhead[ibin+stencil_full[k]]; - while (j >= 0) { - if (i == j) { - j = bins[j]; - continue; - } - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + for (j = binhead[ibin+stencil_full[k]]; j >= 0; j = bins[j]) { + if (i == j) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; jtype = type[j]; delx = xtmp - x[j][0]; @@ -154,8 +147,6 @@ void Neighbor::full_bin() if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } - - j = bins[j]; } } diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp index 21d9b38085..5662c02a96 100644 --- a/src/neigh_gran.cpp +++ b/src/neigh_gran.cpp @@ -290,17 +290,9 @@ void Neighbor::granular_bin_no_newton() // stores own/ghost pairs on both procs for (k = 0; k < nstencil; k++) { - j = binhead[ibin+stencil[k]]; - while (j >= 0) { - if (j <= i) { - j = bins[j]; - continue; - } - - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -337,8 +329,6 @@ void Neighbor::granular_bin_no_newton() n++; } - - j = bins[j]; } } @@ -358,8 +348,8 @@ void Neighbor::granular_bin_no_newton() granular particles binned neighbor list construction with full Newton's 3rd law no shear history is allowed for this option - every pair stored exactly once by some processor each owned atom i checks its own bin and other bins in Newton stencil + every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::granular_bin_newton() @@ -405,19 +395,13 @@ void Neighbor::granular_bin_newton() // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i - j = bins[i]; - while (j >= 0) { - if (j >= nlocal) { - if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) || - (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) { - j = bins[j]; - continue; - } - } - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; } delx = xtmp - x[j][0]; @@ -428,20 +412,93 @@ void Neighbor::granular_bin_newton() cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) neighptr[n++] = j; - - j = bins[j]; } // loop over all atoms in other bins in stencil, store every pair ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { - j = binhead[ibin+stencil[k]]; - while (j >= 0) { - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (exclude && exclusion(i,j,type,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutsq) neighptr[n++] = j; + } + } + + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } +} + +/* ---------------------------------------------------------------------- + granular particles + binned neighbor list construction with Newton's 3rd law for triclinic + no shear history is allowed for this option + each owned atom i checks its own bin and other bins in triclinic stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::granular_bin_newton_tri() +{ + int i,j,k,n,ibin; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutsq; + int *neighptr; + + // bin local & ghost atoms + + bin_atoms(); + + // loop over each atom, storing neighbors + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + + int npage = 0; + int npnt = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == maxpage) add_pages(npage); + } + + n = 0; + neighptr = &pages[npage][npnt]; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in bins in stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and <= x) + // this excludes self-self interaction + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -451,8 +508,6 @@ void Neighbor::granular_bin_newton() cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) neighptr[n++] = j; - - j = bins[j]; } } diff --git a/src/neigh_half.cpp b/src/neigh_half.cpp index af0d3c0f11..37ca5ec665 100644 --- a/src/neigh_half.cpp +++ b/src/neigh_half.cpp @@ -219,17 +219,9 @@ void Neighbor::half_bin_no_newton() // stores own/ghost pairs on both procs for (k = 0; k < nstencil; k++) { - j = binhead[ibin+stencil[k]]; - while (j >= 0) { - if (j <= i) { - j = bins[j]; - continue; - } - - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; jtype = type[j]; delx = xtmp - x[j][0]; @@ -243,8 +235,6 @@ void Neighbor::half_bin_no_newton() if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } - - j = bins[j]; } } @@ -305,20 +295,14 @@ void Neighbor::half_bin_newton() // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i - j = bins[i]; - while (j >= 0) { + for (j = bins[i]; j >= 0; j = bins[j]) { if (j >= nlocal) { - if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) || - (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) { - j = bins[j]; - continue; - } + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + if (exclude && exclusion(i,j,type,mask,molecule)) continue; jtype = type[j]; delx = xtmp - x[j][0]; @@ -332,20 +316,14 @@ void Neighbor::half_bin_newton() if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } - - j = bins[j]; } // loop over all atoms in other bins in stencil, store every pair ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { - j = binhead[ibin+stencil[k]]; - while (j >= 0) { - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (exclude && exclusion(i,j,type,mask,molecule)) continue; jtype = type[j]; delx = xtmp - x[j][0]; @@ -359,8 +337,87 @@ void Neighbor::half_bin_newton() if (which == 0) neighptr[n++] = j; else if (which > 0) neighptr[n++] = which*nall + j; } + } + } - j = bins[j]; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::half_bin_newton_tri() +{ + int i,j,k,n,itype,jtype,ibin,which; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr; + + // bin local & ghost atoms + + bin_atoms(); + + // loop over each atom, storing neighbors + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int npage = 0; + int npnt = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == maxpage) add_pages(npage); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in bins in stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and <= x) + // this excludes self-self interaction + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; + + jtype = type[j]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } } } @@ -460,11 +517,9 @@ void Neighbor::half_full_newton() if (j < nlocal) { if (i > j) continue; } else { - if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) || - (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) { - j = bins[j]; - continue; - } + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } neighptr[n++] = j; } diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp index c27d8e6ba2..a37275f37b 100644 --- a/src/neigh_respa.cpp +++ b/src/neigh_respa.cpp @@ -342,17 +342,9 @@ void Neighbor::respa_bin_no_newton() // stores own/ghost pairs on both procs for (k = 0; k < nstencil; k++) { - j = binhead[ibin+stencil[k]]; - while (j >= 0) { - if (j <= i) { - j = bins[j]; - continue; - } - - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; jtype = type[j]; delx = xtmp - x[j][0]; @@ -376,8 +368,6 @@ void Neighbor::respa_bin_no_newton() else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } - - j = bins[j]; } } @@ -406,8 +396,8 @@ void Neighbor::respa_bin_no_newton() /* ---------------------------------------------------------------------- multiple respa lists binned neighbor list construction with full Newton's 3rd law - every pair stored exactly once by some processor each owned atom i checks its own bin and other bins in Newton stencil + every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void Neighbor::respa_bin_newton() @@ -477,19 +467,12 @@ void Neighbor::respa_bin_newton() // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i - j = bins[i]; - while (j >= 0) { + for (j = bins[i]; j >= 0; j = bins[j]) { if (j >= nlocal) { - if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) || - (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) { - j = bins[j]; - continue; - } - } - - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; } jtype = type[j]; @@ -514,20 +497,144 @@ void Neighbor::respa_bin_newton() else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } - - j = bins[j]; } // loop over all atoms in other bins in stencil, store every pair ibin = coord2bin(x[i]); for (k = 0; k < nstencil; k++) { - j = binhead[ibin+stencil[k]]; - while (j >= 0) { - if (exclude && exclusion(i,j,type,mask,molecule)) { - j = bins[j]; - continue; - } + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (exclude && exclusion(i,j,type,mask,molecule)) continue; + + jtype = type[j]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + + if (rsq < cut_inner_sq) { + if (which == 0) neighptr_inner[n_inner++] = j; + else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; + } + + if (respa == 2 && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { + if (which == 0) neighptr_middle[n_middle++] = j; + else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; + } + } + } + } + + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + + firstneigh_inner[i] = neighptr_inner; + numneigh_inner[i] = n_inner; + npnt_inner += n_inner; + if (npnt_inner >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + + if (respa == 2) { + firstneigh_middle[i] = neighptr_middle; + numneigh_middle[i] = n_middle; + npnt_middle += n_middle; + if (npnt_middle >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + } +} + +/* ---------------------------------------------------------------------- + multiple respa lists + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::respa_bin_newton_tri() +{ + int i,j,k,itype,jtype,ibin,which; + int n_inner,n_middle,n; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr_inner; + int *neighptr_middle; + int *neighptr; + + // bin local & ghost atoms + + bin_atoms(); + + // loop over each atom, storing neighbors + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int npage = 0; + int npnt = 0; + int npage_inner = 0; + int npnt_inner = 0; + int npage_middle = 0; + int npnt_middle = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == maxpage) add_pages(npage); + } + neighptr = &pages[npage][npnt]; + n = 0; + + if (pgsize - npnt_inner < oneatom) { + npnt_inner = 0; + npage_inner++; + if (npage_inner == maxpage_inner) add_pages_inner(npage_inner); + } + neighptr_inner = &pages_inner[npage_inner][npnt_inner]; + n_inner = 0; + + if (respa == 2) { + if (pgsize - npnt_middle < oneatom) { + npnt_middle = 0; + npage_middle++; + if (npage_middle == maxpage_middle) add_pages_middle(npage_middle); + } + neighptr_middle = &pages_middle[npage_middle][npnt_middle]; + n_middle = 0; + } + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in bins in stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and <= x) + // this excludes self-self interaction + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue; + if (exclude && exclusion(i,j,type,mask,molecule)) continue; jtype = type[j]; delx = xtmp - x[j][0]; @@ -551,8 +658,6 @@ void Neighbor::respa_bin_newton() else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; } } - - j = bins[j]; } } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index aadd4cac86..dea05f5244 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -206,6 +206,7 @@ void Neighbor::init() int i,j,m,n; ncalls = ndanger = 0; + triclinic = domain->triclinic; // error check @@ -544,7 +545,7 @@ void Neighbor::init() cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin); } - // set ptrs to correct half and full build functions + // set ptrs to correct half, full, triclinic build functions // cannot combine granular and rRESPA if (half) { @@ -556,6 +557,8 @@ void Neighbor::init() } else if (style == 1) { if (force->newton_pair == 0) half_build = &Neighbor::granular_bin_no_newton; + else if (triclinic) + half_build = &Neighbor::granular_bin_newton_tri; else half_build = &Neighbor::granular_bin_newton; } } else if (respa) { @@ -566,6 +569,8 @@ void Neighbor::init() } else if (style == 1) { if (force->newton_pair == 0) half_build = &Neighbor::respa_bin_no_newton; + else if (triclinic) + half_build = &Neighbor::respa_bin_newton_tri; else half_build = &Neighbor::respa_bin_newton; } } else { @@ -583,6 +588,8 @@ void Neighbor::init() else half_build = &Neighbor::half_bin_no_newton; } else { if (full_every) half_build = &Neighbor::half_full_newton; + else if (triclinic) + half_build = &Neighbor::half_bin_newton_tri; else half_build = &Neighbor::half_bin_newton; } } @@ -869,6 +876,7 @@ void Neighbor::build_full() binsize = 1/2 of cutoff is roughly optimal prd must be filled exactly by integer # of bins so procs on both sides of PBC see same bin boundary + triclinic is an exception, since stencil & neigh list built differently mbinlo = lowest global bin any of my ghost atoms could fall into mbinhi = highest global bin any of my ghost atoms could fall into mbin = number of bins I need in a dimension @@ -878,31 +886,53 @@ void Neighbor::build_full() void Neighbor::setup_bins() { - double cutneighinv = 1.0/cutneigh; + // bbox = bounding box of entire domain + // bboxlo,bboxhi = bounding box of proc's sub-domain + // for triclinic, bounding box surrounds all 8 corner pts + + double bbox[3]; + double *bboxlo,*bboxhi; + + if (triclinic == 0) { + boxlo = domain->boxlo; + boxhi = domain->boxhi; + bboxlo = domain->sublo; + bboxhi = domain->subhi; + } else { + boxlo = domain->boxlo_bound; + boxhi = domain->boxhi_bound; + bboxlo = domain->sublo_bound; + bboxhi = domain->subhi_bound; + } + + bbox[0] = boxhi[0] - boxlo[0]; + bbox[1] = boxhi[1] - boxlo[1]; + bbox[2] = boxhi[2] - boxlo[2]; // test for too many global bins in any dimension due to huge domain - if (2.0*domain->xprd*cutneighinv > INT_MAX || - 2.0*domain->yprd*cutneighinv > INT_MAX || - 2.0*domain->zprd*cutneighinv > INT_MAX) + double cutneighinv = 1.0/cutneigh; + + if (2.0*bbox[0]*cutneighinv > INT_MAX || 2.0*bbox[1]*cutneighinv > INT_MAX || + 2.0*bbox[2]*cutneighinv > INT_MAX) error->all("Domain too large for neighbor bins"); // divide box into bins // optimal size is roughly 1/2 the cutoff - nbinx = static_cast (2.0*domain->xprd*cutneighinv); - nbiny = static_cast (2.0*domain->yprd*cutneighinv); + nbinx = static_cast (2.0*bbox[0]*cutneighinv); + nbiny = static_cast (2.0*bbox[1]*cutneighinv); if (force->dimension == 3) - nbinz = static_cast (2.0*domain->zprd*cutneighinv); + nbinz = static_cast (2.0*bbox[2]*cutneighinv); else nbinz = 1; if (nbinx == 0) nbinx = 1; if (nbiny == 0) nbiny = 1; if (nbinz == 0) nbinz = 1; - binsizex = domain->xprd/nbinx; - binsizey = domain->yprd/nbiny; - binsizez = domain->zprd/nbinz; + binsizex = bbox[0]/nbinx; + binsizey = bbox[1]/nbiny; + binsizez = bbox[2]/nbinz; bininvx = 1.0 / binsizex; bininvy = 1.0 / binsizey; @@ -915,23 +945,23 @@ void Neighbor::setup_bins() double coord; int mbinxhi,mbinyhi,mbinzhi; - coord = domain->subxlo - cutneigh - SMALL*domain->xprd; - mbinxlo = static_cast ((coord-domain->boxxlo)*bininvx); + coord = bboxlo[0] - cutneigh - SMALL*bbox[0]; + mbinxlo = static_cast ((coord-boxlo[0])*bininvx); if (coord < 0.0) mbinxlo = mbinxlo - 1; - coord = domain->subxhi + cutneigh + SMALL*domain->xprd; - mbinxhi = static_cast ((coord-domain->boxxlo)*bininvx); + coord = bboxhi[0] + cutneigh + SMALL*bbox[0]; + mbinxhi = static_cast ((coord-boxlo[0])*bininvx); - coord = domain->subylo - cutneigh - SMALL*domain->yprd; - mbinylo = static_cast ((coord-domain->boxylo)*bininvy); + coord = bboxlo[1] - cutneigh - SMALL*bbox[1]; + mbinylo = static_cast ((coord-boxlo[1])*bininvy); if (coord < 0.0) mbinylo = mbinylo - 1; - coord = domain->subyhi + cutneigh + SMALL*domain->yprd; - mbinyhi = static_cast ((coord-domain->boxylo)*bininvy); + coord = bboxhi[1] + cutneigh + SMALL*bbox[1]; + mbinyhi = static_cast ((coord-boxlo[1])*bininvy); - coord = domain->subzlo - cutneigh - SMALL*domain->zprd; - mbinzlo = static_cast ((coord-domain->boxzlo)*bininvz); + coord = bboxlo[2] - cutneigh - SMALL*bbox[2]; + mbinzlo = static_cast ((coord-boxlo[2])*bininvz); if (coord < 0.0) mbinzlo = mbinzlo - 1; - coord = domain->subzhi + cutneigh + SMALL*domain->zprd; - mbinzhi = static_cast ((coord-domain->boxzlo)*bininvz); + coord = bboxhi[2] + cutneigh + SMALL*bbox[2]; + mbinzhi = static_cast ((coord-boxlo[2])*bininvz); // extend bins by 1 to insure stencil extent is included @@ -965,10 +995,15 @@ void Neighbor::setup_bins() // is within neighbor cutoff // next(xyz) = how far the stencil could possibly extend // for partial Newton (newton = 0) - // stencil is all surrounding bins including self + // stencil is all surrounding bins + // stencil includes self + // for triclinic + // stencil is all bins in z-plane of self and above, but not below + // in 2d is all bins in y-plane of self and above, but not below + // stencil includes self // for full Newton (newton = 1) // stencil is bins to the "upper right" of central bin - // stencil does NOT include self + // stencil does not include self // for full neighbor list (full = 1) // stencil is all surrounding bins including self, regardless of Newton // stored in stencil_full @@ -999,6 +1034,12 @@ void Neighbor::setup_bins() for (i = -nextx; i <= nextx; i++) if (bin_distance(i,j,k) < cutsq) stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i; + } else if (triclinic) { + for (k = 0; k <= nextz; k++) + for (j = -nexty; j <= nexty; j++) + for (i = -nextx; i <= nextx; i++) + if (bin_distance(i,j,k) < cutsq) + stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i; } else { for (k = 0; k <= nextz; k++) for (j = -nexty; j <= nexty; j++) @@ -1013,6 +1054,11 @@ void Neighbor::setup_bins() for (i = -nextx; i <= nextx; i++) if (bin_distance(i,j,0) < cutsq) stencil[nstencil++] = j*mbinx + i; + } else if (triclinic) { + for (j = 0; j <= nexty; j++) + for (i = -nextx; i <= nextx; i++) + if (bin_distance(i,j,0) < cutsq) + stencil[nstencil++] = j*mbinx + i; } else { for (j = 0; j <= nexty; j++) for (i = -nextx; i <= nextx; i++) @@ -1373,32 +1419,33 @@ void Neighbor::bin_atoms() take special care to insure ghosts are put in correct bins this is necessary so that both procs on either side of PBC treat a pair of atoms straddling the PBC in a consistent way + triclinic is an exception, since stencil & neigh list built differently ------------------------------------------------------------------------- */ int Neighbor::coord2bin(double *x) { int ix,iy,iz; - if (x[0] >= domain->boxxhi) - ix = static_cast ((x[0]-domain->boxxhi)*bininvx) + nbinx - mbinxlo; - else if (x[0] >= domain->boxxlo) - ix = static_cast ((x[0]-domain->boxxlo)*bininvx) - mbinxlo; + if (x[0] >= boxhi[0]) + ix = static_cast ((x[0]-boxhi[0])*bininvx) + nbinx - mbinxlo; + else if (x[0] >= boxlo[0]) + ix = static_cast ((x[0]-boxlo[0])*bininvx) - mbinxlo; else - ix = static_cast ((x[0]-domain->boxxlo)*bininvx) - mbinxlo - 1; + ix = static_cast ((x[0]-boxlo[0])*bininvx) - mbinxlo - 1; - if (x[1] >= domain->boxyhi) - iy = static_cast ((x[1]-domain->boxyhi)*bininvy) + nbiny - mbinylo; - else if (x[1] >= domain->boxylo) - iy = static_cast ((x[1]-domain->boxylo)*bininvy) - mbinylo; + if (x[1] >= boxhi[1]) + iy = static_cast ((x[1]-boxhi[1])*bininvy) + nbiny - mbinylo; + else if (x[1] >= boxlo[1]) + iy = static_cast ((x[1]-boxlo[1])*bininvy) - mbinylo; else - iy = static_cast ((x[1]-domain->boxylo)*bininvy) - mbinylo - 1; + iy = static_cast ((x[1]-boxlo[1])*bininvy) - mbinylo - 1; - if (x[2] >= domain->boxzhi) - iz = static_cast ((x[2]-domain->boxzhi)*bininvz) + nbinz - mbinzlo; - else if (x[2] >= domain->boxzlo) - iz = static_cast ((x[2]-domain->boxzlo)*bininvz) - mbinzlo; + if (x[2] >= boxhi[2]) + iz = static_cast ((x[2]-boxhi[2])*bininvz) + nbinz - mbinzlo; + else if (x[2] >= boxlo[2]) + iz = static_cast ((x[2]-boxlo[2])*bininvz) - mbinzlo; else - iz = static_cast ((x[2]-domain->boxzlo)*bininvz) - mbinzlo - 1; + iz = static_cast ((x[2]-boxlo[2])*bininvz) - mbinzlo - 1; return (iz*mbiny*mbinx + iy*mbinx + ix + 1); } diff --git a/src/neighbor.h b/src/neighbor.h index 794f96109d..a4354b49c0 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -117,6 +117,9 @@ class Neighbor : protected Pointers { double binsizex,binsizey,binsizez; // bin sizes and inverse sizes double bininvx,bininvy,bininvz; + int triclinic; // 0 if domain is orthog, 1 if triclinic + double *boxlo,*boxhi; // copies of domain bounds + int **pages; // half neighbor list pages int maxpage; // # of half pages currently allocated @@ -217,6 +220,12 @@ class Neighbor : protected Pointers { int coord2bin(double *); // mapping atom coord to a bin int find_special(int, int); // look for special neighbor int exclusion(int, int, int *, int *, int *); // test for pair exclusion + + // PJV + + void granular_bin_newton_tri(); + void respa_bin_newton_tri(); + void half_bin_newton_tri(); }; } diff --git a/src/pair.h b/src/pair.h index 5c002819f7..5479078f94 100644 --- a/src/pair.h +++ b/src/pair.h @@ -79,7 +79,7 @@ class Pair : protected Pointers { virtual void write_restart_settings(FILE *) {} virtual void read_restart_settings(FILE *) {} - virtual int pack_comm(int, int *, double *, int *) {return 0;} + virtual int pack_comm(int, int *, double *, int, double *) {return 0;} virtual void unpack_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} diff --git a/src/read_data.cpp b/src/read_data.cpp index a614a01211..539fa759eb 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "math.h" #include "mpi.h" #include "string.h" #include "stdlib.h" @@ -114,6 +115,7 @@ void ReadData::command(int narg, char **arg) atom->allocate_type_arrays(); atom->avec->grow(n); + domain->print_box(" "); domain->set_initial_box(); domain->set_global_box(); comm->set_procs(); @@ -356,7 +358,10 @@ void ReadData::header(int flag) sscanf(line,"%lg %lg",&domain->boxylo,&domain->boxyhi); else if (strstr(line,"zlo zhi")) sscanf(line,"%lg %lg",&domain->boxzlo,&domain->boxzhi); - else break; + else if (strstr(line,"xy xz yz")) { + domain->triclinic = 1; + sscanf(line,"%lg %lg %lg",&domain->xy,&domain->xz,&domain->yz); + } else break; } // error check on consistency of header values @@ -382,10 +387,6 @@ void ReadData::header(int flag) error->one("Dihedrals defined but no dihedral types"); if (atom->nimpropers > 0 && atom->nimpropertypes <= 0) error->one("Impropers defined but no improper types"); - - if (domain->boxxlo >= domain->boxxhi || domain->boxylo >= domain->boxyhi || - domain->boxzlo >= domain->boxzhi) - error->one("Box bounds are invalid"); } /* ---------------------------------------------------------------------- diff --git a/src/read_restart.cpp b/src/read_restart.cpp index 87b057a60e..81b67be97b 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -101,6 +101,7 @@ void ReadRestart::command(int narg, char **arg) atom->allocate_type_arrays(); atom->avec->grow(n); + domain->print_box(" "); domain->set_initial_box(); domain->set_global_box(); comm->set_procs(); @@ -124,21 +125,25 @@ void ReadRestart::command(int narg, char **arg) // proc 0 reads chunks from series of files (nprocs_file) // proc 0 bcasts each chunk to other procs // each proc unpacks the atoms, saving ones in it's sub-domain + // check for atom in sub-domain is different for orthogonal vs triclinic box AtomVec *avec = atom->avec; + int triclinic = domain->triclinic; int maxbuf = 0; double *buf = NULL; - double subxlo = domain->subxlo; - double subxhi = domain->subxhi; - double subylo = domain->subylo; - double subyhi = domain->subyhi; - double subzlo = domain->subzlo; - double subzhi = domain->subzhi; + double *x,lamda[3]; + double *coord,*sublo,*subhi; + if (triclinic == 0) { + sublo = domain->sublo; + subhi = domain->subhi; + } else { + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } int m; - double xtmp,ytmp,ztmp; char *perproc = new char[strlen(file) + 16]; char *ptr = strchr(file,'%'); @@ -171,12 +176,15 @@ void ReadRestart::command(int narg, char **arg) m = 0; while (m < n) { - xtmp = buf[m+1]; - ytmp = buf[m+2]; - ztmp = buf[m+3]; - if (xtmp >= subxlo && xtmp < subxhi && - ytmp >= subylo && ytmp < subyhi && - ztmp >= subzlo && ztmp < subzhi) + x = &buf[m+1]; + if (triclinic) { + domain->x2lamda(x,lamda); + coord = lamda; + } else coord = x; + + if (coord[0] >= sublo[0] && coord[0] < subhi[0] && + coord[1] >= sublo[1] && coord[1] < subhi[1] && + coord[2] >= sublo[2] && coord[2] < subhi[2]) m += avec->unpack_restart(&buf[m]); else m += static_cast (buf[m]); } @@ -520,6 +528,16 @@ void ReadRestart::header() } else if (flag == 45) { force->special_coul[3] = read_double(); + } else if (flag == 46) { + domain->triclinic = 1; + domain->xy = read_double(); + } else if (flag == 47) { + domain->triclinic = 1; + domain->xz = read_double(); + } else if (flag == 48) { + domain->triclinic = 1; + domain->yz = read_double(); + } else error->all("Invalid flag in header of restart file"); flag = read_int(); diff --git a/src/region_block.cpp b/src/region_block.cpp index fa4284981b..61ea4a727e 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -28,37 +28,43 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) if (strcmp(arg[2],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - xlo = domain->boxxlo; + if (domain->triclinic == 0) xlo = domain->boxlo[0]; + else xlo = domain->boxlo_bound[0]; } else xlo = xscale*atof(arg[2]); if (strcmp(arg[3],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - xhi = domain->boxxhi; + if (domain->triclinic == 0) xhi = domain->boxhi[0]; + else xhi = domain->boxhi_bound[0]; } else xhi = xscale*atof(arg[3]); if (strcmp(arg[4],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - ylo = domain->boxylo; + if (domain->triclinic == 0) ylo = domain->boxlo[1]; + else ylo = domain->boxlo_bound[1]; } else ylo = yscale*atof(arg[4]); if (strcmp(arg[5],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - yhi = domain->boxyhi; + if (domain->triclinic == 0) yhi = domain->boxhi[1]; + else yhi = domain->boxhi_bound[1]; } else yhi = yscale*atof(arg[5]); if (strcmp(arg[6],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - zlo = domain->boxzlo; + if (domain->triclinic == 0) zlo = domain->boxlo[2]; + else zlo = domain->boxlo_bound[2]; } else zlo = zscale*atof(arg[6]); if (strcmp(arg[7],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - zhi = domain->boxzhi; + if (domain->triclinic == 0) zhi = domain->boxhi[2]; + else zhi = domain->boxhi_bound[2]; } else zhi = zscale*atof(arg[7]); // error check diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp index 3bac413c89..012c2ecc3a 100644 --- a/src/region_cylinder.cpp +++ b/src/region_cylinder.cpp @@ -46,9 +46,18 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[6],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - if (axis == 'x') lo = domain->boxxlo; - if (axis == 'y') lo = domain->boxylo; - if (axis == 'z') lo = domain->boxzlo; + if (axis == 'x') { + if (domain->triclinic == 0) lo = domain->boxlo[0]; + else lo = domain->boxlo_bound[0]; + } + if (axis == 'y') { + if (domain->triclinic == 0) lo = domain->boxlo[1]; + else lo = domain->boxlo_bound[1]; + } + if (axis == 'z') { + if (domain->triclinic == 0) lo = domain->boxlo[2]; + else lo = domain->boxlo_bound[2]; + } } else { if (axis == 'x') lo = xscale*atof(arg[6]); if (axis == 'y') lo = yscale*atof(arg[6]); @@ -58,9 +67,18 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[7],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - if (axis == 'x') hi = domain->boxxhi; - if (axis == 'y') hi = domain->boxyhi; - if (axis == 'z') hi = domain->boxzhi; + if (axis == 'x') { + if (domain->triclinic == 0) hi = domain->boxhi[0]; + else hi = domain->boxhi_bound[0]; + } + if (axis == 'y') { + if (domain->triclinic == 0) hi = domain->boxhi[1]; + else hi = domain->boxhi_bound[1]; + } + if (axis == 'z') { + if (domain->triclinic == 0) hi = domain->boxhi[2]; + else hi = domain->boxhi_bound[2]; + } } else { if (axis == 'x') hi = xscale*atof(arg[7]); if (axis == 'y') hi = yscale*atof(arg[7]); diff --git a/src/region_prism.cpp b/src/region_prism.cpp index 33a56a4d48..b2206960eb 100644 --- a/src/region_prism.cpp +++ b/src/region_prism.cpp @@ -36,42 +36,42 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) if (strcmp(arg[2],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - xlo = domain->boxxlo; + xlo = domain->boxlo[0]; } else xlo = xscale*atof(arg[2]); if (strcmp(arg[3],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - xhi = domain->boxxhi; + xhi = domain->boxhi[0]; } else xhi = xscale*atof(arg[3]); if (strcmp(arg[4],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - ylo = domain->boxylo; + ylo = domain->boxlo[1]; } else ylo = yscale*atof(arg[4]); if (strcmp(arg[5],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - yhi = domain->boxyhi; + yhi = domain->boxhi[1]; } else yhi = yscale*atof(arg[5]); if (strcmp(arg[6],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - zlo = domain->boxzlo; + zlo = domain->boxlo[0]; } else zlo = zscale*atof(arg[6]); if (strcmp(arg[7],"INF") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF when box does not exist"); - zhi = domain->boxzhi; + zhi = domain->boxhi[1]; } else zhi = zscale*atof(arg[7]); - yxshift = xscale*atof(arg[8]); - zxshift = xscale*atof(arg[9]); - zyshift = yscale*atof(arg[10]); + xy = xscale*atof(arg[8]); + xz = xscale*atof(arg[9]); + yz = yscale*atof(arg[10]); // error check // prism cannot be 0 thickness in any dim, else inverse blows up @@ -81,13 +81,14 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) // extent of prism - extent_xlo = MIN(xlo,xlo+yxshift); - extent_xlo = MIN(extent_xlo,extent_xlo+zxshift); - extent_xhi = MAX(xhi,xhi+yxshift); - extent_xhi = MAX(extent_xhi,extent_xhi+zxshift); - extent_ylo = MIN(ylo,ylo+zyshift); - extent_yhi = MAX(yhi,yhi+zyshift); + extent_xlo = MIN(xlo,xlo+xy); + extent_xlo = MIN(extent_xlo,extent_xlo+xz); + extent_ylo = MIN(ylo,ylo+yz); extent_zlo = zlo; + + extent_xhi = MAX(xhi,xhi+xy); + extent_xhi = MAX(extent_xhi,extent_xhi+xz); + extent_yhi = MAX(yhi,yhi+yz); extent_zhi = zhi; // h = transformation matrix from tilt coords (0-1) to box coords (xyz) @@ -98,27 +99,27 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) // and bottom face of prism is in xy plane h[0][0] = xhi - xlo; - h[0][1] = yxshift; + h[0][1] = xy; + h[0][2] = xz; h[1][1] = yhi - ylo; - h[0][2] = zxshift; - h[1][2] = zyshift; + h[1][2] = yz; h[2][2] = zhi - zlo; hinv[0][0] = 1.0/h[0][0]; hinv[0][1] = -h[0][1] / (h[0][0]*h[1][1]); - hinv[1][1] = 1.0/h[1][1]; hinv[0][2] = (h[0][1]*h[1][2] - h[0][2]*h[1][1]) / (h[0][0]*h[1][1]*h[2][2]); + hinv[1][1] = 1.0/h[1][1]; hinv[1][2] = -h[1][2] / (h[1][1]*h[2][2]); hinv[2][2] = 1.0/h[2][2]; } /* ---------------------------------------------------------------------- check xyz against prism - abc = Hinv * (xyz - xyzlo) + abc = Hinv * (xyz - xyz/lo) abc = tilt coords (0-1) Hinv = transformation matrix from box coords to tilt coords xyz = box coords - xyzlo = lower-left corner of prism + xyz/lo = lower-left corner of prism ------------------------------------------------------------------------- */ int RegPrism::match(double x, double y, double z) diff --git a/src/region_prism.h b/src/region_prism.h index 98af0d96ba..2bdae8c9cd 100644 --- a/src/region_prism.h +++ b/src/region_prism.h @@ -19,13 +19,15 @@ namespace LAMMPS_NS { class RegPrism : public Region { + friend class CreateBox; + public: RegPrism(class LAMMPS *, int, char **); int match(double, double, double); private: double xlo,xhi,ylo,yhi,zlo,zhi; - double yxshift,zxshift,zyshift; + double xy,xz,yz; double h[3][3],hinv[3][3]; }; diff --git a/src/replicate.cpp b/src/replicate.cpp index 31138e7d2b..97c8235962 100644 --- a/src/replicate.cpp +++ b/src/replicate.cpp @@ -27,7 +27,8 @@ using namespace LAMMPS_NS; #define LB_FACTOR 1.1 -#define MAXATOMS 0x7FFFFFFF +#define MAXATOMS 0x7FFFFFFF +#define EPSILON 1.0e-6 #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) @@ -92,7 +93,7 @@ void Replicate::command(int narg, char **arg) // unmap existing atoms via image flags for (i = 0; i < atom->nlocal; i++) - domain->unmap(atom->x[i][0],atom->x[i][1],atom->x[i][2],atom->image[i]); + domain->unmap(atom->x[i],atom->image[i]); // communication buffer for all my atom's info // max_size = largest buffer needed by any proc @@ -157,15 +158,24 @@ void Replicate::command(int narg, char **arg) // store old simulation box + int triclinic = domain->triclinic; double old_xprd = domain->xprd; double old_yprd = domain->yprd; double old_zprd = domain->zprd; + double old_xy = domain->xy; + double old_xz = domain->xz; + double old_yz = domain->yz; // setup new simulation box domain->boxxhi = domain->boxxlo + nx*old_xprd; domain->boxyhi = domain->boxylo + ny*old_yprd; domain->boxzhi = domain->boxzlo + nz*old_zprd; + if (triclinic) { + domain->xy *= ny; + domain->xz *= nz; + domain->yz *= nz; + } // new problem setup using new box boundaries @@ -175,20 +185,12 @@ void Replicate::command(int narg, char **arg) atom->allocate_type_arrays(); atom->avec->grow(n); + domain->print_box(" "); domain->set_initial_box(); domain->set_global_box(); comm->set_procs(); domain->set_local_box(); - // sub xyz lo/hi = new processor sub-domain - - double subxlo = domain->subxlo; - double subxhi = domain->subxhi; - double subylo = domain->subylo; - double subyhi = domain->subyhi; - double subzlo = domain->subzlo; - double subzhi = domain->subzhi; - // copy type arrays to new atom class if (atom->mass) { @@ -206,12 +208,41 @@ void Replicate::command(int narg, char **arg) } } + // set bounds for my proc + // if periodic and I am lo/hi proc, adjust bounds by EPSILON + // insures all replicated atoms will be owned even with round-off + + double sublo[3],subhi[3]; + + if (triclinic == 0) { + sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0]; + sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1]; + sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2]; + } else { + sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0]; + sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1]; + sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2]; + } + + if (domain->xperiodic) { + if (comm->myloc[0] == 0) sublo[0] -= EPSILON; + if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += EPSILON; + } + if (domain->yperiodic) { + if (comm->myloc[1] == 0) sublo[1] -= EPSILON; + if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += EPSILON; + } + if (domain->zperiodic) { + if (comm->myloc[2] == 0) sublo[2] -= EPSILON; + if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += EPSILON; + } + // loop over all procs // if this iteration of loop is me: // pack my unmapped atom data into buf // bcast it to all other procs - // for each atom in buf, each proc performs 3d replicate loop: - // xnew,ynew,znew = new replicated position + // performs 3d replicate loop with while loop over atoms in buf + // x = new replicated position // unpack atom into new atom class from buf if I own it // adjust tag, mol #, coord, topology info as needed @@ -219,7 +250,8 @@ void Replicate::command(int narg, char **arg) AtomVec *avec = atom->avec; int ix,iy,iz,image,atom_offset,mol_offset; - double xnew,ynew,znew; + double x[3],lamda[3]; + double *coord; int tag_enable = atom->tag_enable; for (int iproc = 0; iproc < nprocs; iproc++) { @@ -234,17 +266,29 @@ void Replicate::command(int narg, char **arg) for (iy = 0; iy < ny; iy++) { for (iz = 0; iz < nz; iz++) { + // while loop over one proc's atom list + m = 0; while (m < n) { - xnew = buf[m+1] + ix*old_xprd; - ynew = buf[m+2] + iy*old_yprd; - znew = buf[m+3] + iz*old_zprd; image = (512 << 20) | (512 << 10) | 512; - domain->remap(xnew,ynew,znew,image); + if (triclinic == 0) { + x[0] = buf[m+1] + ix*old_xprd; + x[1] = buf[m+2] + iy*old_yprd; + x[2] = buf[m+3] + iz*old_zprd; + } else { + x[0] = buf[m+1] + ix*old_xprd + iy*old_xy + iz*old_xz; + x[1] = buf[m+2] + iy*old_yprd + iz*old_yz; + x[2] = buf[m+3] + iz*old_zprd; + } + domain->remap(x,image); + if (triclinic) { + domain->x2lamda(x,lamda); + coord = lamda; + } else coord = x; - if (xnew >= subxlo && xnew < subxhi && - ynew >= subylo && ynew < subyhi && - znew >= subzlo && znew < subzhi) { + if (coord[0] >= sublo[0] && coord[0] < subhi[0] && + coord[1] >= sublo[1] && coord[1] < subhi[1] && + coord[2] >= sublo[2] && coord[2] < subhi[2]) { m += avec->unpack_restart(&buf[m]); @@ -254,9 +298,9 @@ void Replicate::command(int narg, char **arg) else atom_offset = 0; mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol; - atom->x[i][0] = xnew; - atom->x[i][1] = ynew; - atom->x[i][2] = znew; + atom->x[i][0] = x[0]; + atom->x[i][1] = x[1]; + atom->x[i][2] = x[2]; atom->tag[i] += atom_offset; atom->image[i] = image; @@ -289,8 +333,7 @@ void Replicate::command(int narg, char **arg) } } } else m += static_cast (buf[m]); - } // end of while loop over one proc's atom list - + } } } } diff --git a/src/respa.cpp b/src/respa.cpp index becb5949ce..79a51be076 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -321,6 +321,10 @@ void Respa::init() newton[ilevel] = 1; } } + + // orthogonal vs triclinic simulation box + + triclinic = domain->triclinic; } /* ---------------------------------------------------------------------- @@ -335,12 +339,14 @@ void Respa::setup() // acquire ghosts // build neighbor lists + if (triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); comm->exchange(); comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); neighbor->build(); neighbor->ncalls = 0; @@ -464,6 +470,7 @@ void Respa::recurse(int ilevel) int nflag = neighbor->decide(); if (nflag) { if (modify->n_pre_exchange) modify->pre_exchange(); + if (triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); if (domain->box_change) { domain->reset_box(); @@ -473,6 +480,7 @@ void Respa::recurse(int ilevel) timer->stamp(); comm->exchange(); comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); timer->stamp(TIME_COMM); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(); diff --git a/src/respa.h b/src/respa.h index e592f77fd8..6cb0aa06e4 100644 --- a/src/respa.h +++ b/src/respa.h @@ -51,6 +51,7 @@ class Respa : public Integrate { int nfix_virial; // # of fixes that need virial occasionally int *fix_virial_every; // frequency they require it int *next_fix_virial; // next timestep they need it + int triclinic; // 0 if domain is orthog, 1 if triclinic int *newton; // newton flag at each level class FixRespa *fix_respa; // Fix to store the force level array diff --git a/src/set.cpp b/src/set.cpp index 8cb475ed62..de7ed0a4a0 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -103,11 +103,13 @@ void Set::command(int narg, char **arg) if (comm->me == 0 && screen) fprintf(screen,"System init for set ...\n"); lmp->init(); + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); comm->exchange(); comm->borders(); + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); } if (comm->me == 0 && screen) fprintf(screen,"Setting atom values ...\n"); diff --git a/src/velocity.cpp b/src/velocity.cpp index 9d03e31979..7f84f3e239 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -283,7 +283,7 @@ void Velocity::create(int narg, char **arg) for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - triple(x[i][0],x[i][1],x[i][2],&vx,&vy,&vz,seed,random); + triple(x[i],&vx,&vy,&vz,seed,random); if (mass) factor = 1.0/sqrt(mass[type[i]]); else factor = 1.0/sqrt(rmass[i]); v[i][0] = vx * factor; @@ -661,31 +661,37 @@ void Velocity::options(int narg, char **arg) #define IC3 54773 #define IM3 259200 -void Velocity::triple(double x, double y, double z, - double *vx, double *vy, double *vz, +void Velocity::triple(double *x, double *vx, double *vy, double *vz, int seed, RanPark *random) { + // lamda[3] = fractional position in box + // for triclinic box, just convert to lamda coords + + double lamda[3]; + + if (domain->triclinic == 0) { + lamda[0] = (x[0] - domain->boxlo[0]) / domain->prd[0]; + lamda[1] = (x[1] - domain->boxlo[1]) / domain->prd[1]; + lamda[2] = (x[2] - domain->boxlo[2]) / domain->prd[2]; + } else domain->x2lamda(x,lamda); + // seed 1,2,3 = combination of atom coord in each dim and user-input seed // map geometric extent into range of each of 3 RNGs // warm-up each RNG by calling it twice - double fraction; int seed1,seed2,seed3; - fraction = (x - domain->boxxlo) / domain->xprd; - seed1 = static_cast (fraction * IM1); + seed1 = static_cast (lamda[0] * IM1); seed1 = (seed1+seed) % IM1; seed1 = (seed1*IA1+IC1) % IM1; seed1 = (seed1*IA1+IC1) % IM1; - fraction = (y - domain->boxylo) / domain->yprd; - seed2 = static_cast (fraction * IM2); + seed2 = static_cast (lamda[1] * IM2); seed2 = (seed2+seed) % IM2; seed2 = (seed2*IA2+IC2) % IM2; seed2 = (seed2*IA2+IC2) % IM2; - fraction = (z - domain->boxzlo) / domain->zprd; - seed3 = static_cast (fraction * IM3); + seed3 = static_cast (lamda[2] * IM3); seed3 = (seed3+seed) % IM3; seed3 = (seed3*IA3+IC3) % IM3; seed3 = (seed3*IA3+IC3) % IM3; @@ -693,7 +699,7 @@ void Velocity::triple(double x, double y, double z, // fraction = 0-1 with giving each dim an equal weighting // use fraction to reset Park/Miller RNG seed - fraction = 1.0*seed1/(3*IM1) + 1.0*seed2/(3*IM2) + 1.0*seed3/(3*IM3); + double fraction = 1.0*seed1/(3*IM1) + 1.0*seed2/(3*IM2) + 1.0*seed3/(3*IM3); random->reset(fraction); // use RNG to set velocities after warming up twice diff --git a/src/velocity.h b/src/velocity.h index a1a0795321..bc11f863c2 100644 --- a/src/velocity.h +++ b/src/velocity.h @@ -41,7 +41,7 @@ class Velocity : protected Pointers { void zero_rotation(); void options(int, char **); - void triple(double, double, double, double *, double *, double *, + void triple(double *, double *, double *, double *, int, class RanPark *); }; diff --git a/src/verlet.cpp b/src/verlet.cpp index 8342407cfb..4154bfe4bd 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -103,6 +103,10 @@ void Verlet::init() pairflag = 1; if (strcmp(atom->atom_style,"granular") == 0) pairflag = 0; + // orthogonal vs triclinic simulation box + + triclinic = domain->triclinic; + // local copies of Update quantities maxpair = update->maxpair; @@ -121,12 +125,14 @@ void Verlet::setup() // acquire ghosts // build neighbor lists + if (triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); comm->exchange(); comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); neighbor->build(); neighbor->ncalls = 0; @@ -197,6 +203,7 @@ void Verlet::iterate(int n) timer->stamp(TIME_COMM); } else { if (modify->n_pre_exchange) modify->pre_exchange(); + if (triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); if (domain->box_change) { domain->reset_box(); @@ -206,6 +213,7 @@ void Verlet::iterate(int n) timer->stamp(); comm->exchange(); comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); timer->stamp(TIME_COMM); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(); diff --git a/src/verlet.h b/src/verlet.h index 1aae82c7f6..b05465bad2 100644 --- a/src/verlet.h +++ b/src/verlet.h @@ -33,6 +33,7 @@ class Verlet : public Integrate { int nfix_virial; // # of fixes that need virial occasionally int *fix_virial_every; // frequency they require it int *next_fix_virial; // next timestep they need it + int triclinic; // 0 if domain is orthog, 1 if triclinic int pairflag,torqueflag,granflag; // arrays to zero out every step int maxpair; // local copies of Update quantities diff --git a/src/write_restart.cpp b/src/write_restart.cpp index 5a87e07356..16cccc5973 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -77,10 +77,12 @@ void WriteRestart::command(int narg, char **arg) // move atoms to new processors before writing file // enforce PBC before in case atoms are outside box + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); comm->exchange(); + if (domain->triclinic) domain->lamda2x(atom->nlocal); write(file); delete [] file; @@ -284,6 +286,12 @@ void WriteRestart::header() write_double(44,force->special_coul[2]); write_double(45,force->special_coul[3]); + if (domain->triclinic) { + write_double(46,domain->xy); + write_double(47,domain->xz); + write_double(48,domain->yz); + } + // -1 flag signals end of header int flag = -1;