From c81b8c5d042d4556096da5c76c6596c2698402eb Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 6 Feb 2012 16:58:46 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7702 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/change_box.cpp | 441 ++++++++++++++++++++++++++++++++++++++++--- src/change_box.h | 61 +++++- src/displace_box.cpp | 418 ---------------------------------------- src/displace_box.h | 105 ----------- src/domain.cpp | 31 ++- src/domain.h | 3 +- src/input.cpp | 2 +- src/pair_hybrid.cpp | 2 +- src/read_restart.cpp | 6 +- 9 files changed, 503 insertions(+), 566 deletions(-) delete mode 100644 src/displace_box.cpp delete mode 100644 src/displace_box.h diff --git a/src/change_box.cpp b/src/change_box.cpp index 34f2f67b0b..be83182742 100644 --- a/src/change_box.cpp +++ b/src/change_box.cpp @@ -11,17 +11,28 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "lmptype.h" +#include "mpi.h" +#include "math.h" +#include "stdlib.h" #include "string.h" #include "change_box.h" -#include "domain.h" +#include "atom.h" #include "modify.h" #include "fix.h" +#include "domain.h" +#include "lattice.h" +#include "comm.h" +#include "irregular.h" #include "output.h" +#include "group.h" #include "error.h" using namespace LAMMPS_NS; -enum{ORTHO,TRICLINIC}; +enum{XYZ,TILT,BOUNDARY,ORTHO,TRICLINIC,SET,REMAP}; +enum{FINAL,DELTA,SCALE}; +enum{X,Y,Z,YZ,XZ,XY}; /* ---------------------------------------------------------------------- */ @@ -31,34 +42,414 @@ ChangeBox::ChangeBox(LAMMPS *lmp) : Pointers(lmp) {} void ChangeBox::command(int narg, char **arg) { + int i; + if (domain->box_exist == 0) error->all(FLERR,"Change_box command before simulation box is defined"); - if (narg != 1) error->all(FLERR,"Illegal change_box command"); + if (narg < 2) error->all(FLERR,"Illegal change_box command"); + if (modify->nfix_restart_peratom) + error->all(FLERR,"Cannot change_box after " + "reading restart file with per-atom info"); - int style; - if (strcmp(arg[0],"ortho") == 0) style = ORTHO; - else if (strcmp(arg[0],"triclinic") == 0) style = TRICLINIC; - else error->all(FLERR,"Illegal change_box command"); + if (comm->me == 0 && screen) fprintf(screen,"Changing box ...\n"); - if (style == ORTHO && domain->triclinic == 0) - error->all(FLERR,"Change_box operation is invalid"); - if (style == TRICLINIC && domain->triclinic == 1) - error->all(FLERR,"Change_box operation is invalid"); - if (style == ORTHO && - (domain->xy != 0.0 || domain->yz != 0.0 || domain->xz != 0.0)) - error->all(FLERR,"Cannot change box to orthogonal when tilt is non-zero"); + // group - if (output->ndump) - error->all(FLERR,"Cannot change box with dumps defined"); - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->no_change_box) - error->all(FLERR,"Cannot change box with certain fixes defined"); + int igroup = group->find(arg[0]); + if (igroup == -1) error->all(FLERR,"Could not find change_box group ID"); + int groupbit = group->bitmask[igroup]; - if (style == ORTHO) domain->triclinic = 0; - else domain->triclinic = 1; - domain->xy = domain->yz = domain->xz = 0.0; + // parse operation arguments + // allocate ops to max possible length + // volume option does not increment nops - domain->set_global_box(); - if (style == TRICLINIC) domain->set_lamda_box(); - domain->set_local_box(); + int dimension = domain->dimension; + + ops = new Operation[narg-1]; + nops = 0; + + int index; + int iarg = 1; + while (iarg < narg) { + if (strcmp(arg[iarg],"x") == 0 || strcmp(arg[iarg],"y") == 0 || + strcmp(arg[iarg],"z") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].style = XYZ; + if (strcmp(arg[iarg],"x") == 0) ops[nops].dim = X; + else if (strcmp(arg[iarg],"y") == 0) ops[nops].dim = Y; + else if (strcmp(arg[iarg],"z") == 0) ops[nops].dim = Z; + + if (dimension == 2 && ops[nops].dim == Z) + error->all(FLERR,"Cannot change_box in z dimension for 2d simulation"); + + if (strcmp(arg[iarg+1],"final") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].flavor = FINAL; + ops[nops].flo = atof(arg[iarg+2]); + ops[nops].fhi = atof(arg[iarg+3]); + ops[nops].vdim1 = ops[nops].vdim2 = -1; + nops++; + iarg += 4; + } else if (strcmp(arg[iarg+1],"delta") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].flavor = DELTA; + ops[nops].dlo = atof(arg[iarg+2]); + ops[nops].dhi = atof(arg[iarg+3]); + ops[nops].vdim1 = ops[nops].vdim2 = -1; + nops++; + iarg += 4; + } else if (strcmp(arg[iarg+1],"scale") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].flavor = SCALE; + ops[nops].scale = atof(arg[iarg+2]); + ops[nops].vdim1 = ops[nops].vdim2 = -1; + nops++; + iarg += 3; + } else if (strcmp(arg[iarg+1],"volume") == 0) { + if (nops == 0 || ops[nops-1].style != XYZ || + ops[nops].dim == ops[nops-1].dim) + error->all(FLERR,"Change_box volume used incorrectly"); + if (ops[nops-1].vdim2 >= 0) + error->all(FLERR,"Change_box volume used incorrectly"); + else if (ops[nops-1].vdim1 >= 0) ops[nops-1].vdim2 = ops[nops].dim; + else ops[nops-1].vdim1 = ops[nops].dim; + iarg += 2; + + } else error->all(FLERR,"Illegal change_box command"); + + } else if (strcmp(arg[iarg],"xy") == 0 || strcmp(arg[iarg],"xz") == 0 || + strcmp(arg[iarg],"yz") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].style = TILT; + if (strcmp(arg[iarg],"xy") == 0) ops[nops].dim = XY; + else if (strcmp(arg[iarg],"xz") == 0) ops[nops].dim = XZ; + else if (strcmp(arg[iarg],"yz") == 0) ops[nops].dim = YZ; + + if (dimension == 2 && (ops[nops].dim == XZ || ops[nops].dim == YZ)) + error->all(FLERR,"Cannot change_box in xz or yz for 2d simulation"); + + if (strcmp(arg[iarg+1],"final") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].flavor = FINAL; + ops[nops].ftilt = atof(arg[iarg+2]); + nops++; + iarg += 3; + } else if (strcmp(arg[iarg+1],"delta") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].flavor = DELTA; + ops[nops].dtilt = atof(arg[iarg+2]); + nops++; + iarg += 3; + } else error->all(FLERR,"Illegal change_box command"); + + } else if (strcmp(arg[iarg],"boundary") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].style = BOUNDARY; + ops[nops].boundindex = iarg+1; + nops++; + iarg += 4; + + } else if (strcmp(arg[iarg],"ortho") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].style = ORTHO; + nops++; + iarg += 1; + + } else if (strcmp(arg[iarg],"triclinic") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].style = TRICLINIC; + nops++; + iarg += 1; + + } else if (strcmp(arg[iarg],"set") == 0) { + if (iarg+1 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].style = SET; + nops++; + iarg += 1; + + } else if (strcmp(arg[iarg],"remap") == 0) { + if (iarg+1 > narg) error->all(FLERR,"Illegal change_box command"); + ops[nops].style = REMAP; + nops++; + iarg += 1; + + } else break; + } + + if (nops == 0) error->all(FLERR,"Illegal change_box command"); + + // read options from end of input line + + options(narg-iarg,&arg[iarg]); + + // compute scale factors if FINAL,DELTA used since they have distance units + + int flag = 0; + for (int i = 0; i < nops; i++) + if (ops[i].style == FINAL || ops[i].style == DELTA) flag = 1; + + if (flag && scaleflag && domain->lattice == NULL) + error->all(FLERR,"Use of change_box with undefined lattice"); + + if (flag && scaleflag) { + scale[0] = domain->lattice->xlattice; + scale[1] = domain->lattice->ylattice; + scale[2] = domain->lattice->zlattice; + } + else scale[0] = scale[1] = scale[2] = 1.0; + + // perform sequence of operations + // first insure atoms are in current box & update box via shrink-wrap + // then save current box state so can remap atoms from it + + domain->pbc(); + domain->reset_box(); + save_box_state(); + + for (int i = 0; i < nops; i++) { + if (ops[i].style == XYZ) { + double volume; + if (domain->dimension == 2) volume = domain->xprd * domain->yprd; + else volume = domain->xprd * domain->yprd * domain->zprd; + + if (ops[i].flavor == FINAL) { + domain->boxlo[ops[i].dim] = scale[ops[i].dim]*ops[i].flo; + domain->boxhi[ops[i].dim] = scale[ops[i].dim]*ops[i].fhi; + if (ops[i].vdim1) + volume_preserve(ops[i].vdim1,ops[i].vdim2,volume); + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + domain->print_box(" "); + + } else if (ops[i].flavor == DELTA) { + domain->boxlo[ops[i].dim] += scale[ops[i].dim]*ops[i].dlo; + domain->boxhi[ops[i].dim] += scale[ops[i].dim]*ops[i].dhi; + if (ops[i].vdim1) + volume_preserve(ops[i].vdim1,ops[i].vdim2,volume); + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + domain->print_box(" "); + + } else if (ops[i].flavor == SCALE) { + double mid = 0.5 * + (domain->boxlo[ops[i].dim] + domain->boxhi[ops[i].dim]); + double delta = domain->boxlo[ops[i].dim] - mid; + domain->boxlo[ops[i].dim] = mid + ops[i].scale*delta; + delta = domain->boxhi[ops[i].dim] - mid; + domain->boxhi[ops[i].dim] = mid + ops[i].scale*delta; + if (ops[i].vdim1) + volume_preserve(ops[i].vdim1,ops[i].vdim2,volume); + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + domain->print_box(" "); + } + + } else if (ops[i].style == TILT) { + if (domain->triclinic == 0) + error->all(FLERR,"Cannot change box tilt factors for orthogonal box"); + + if (ops[i].flavor == FINAL) { + if (ops[i].dim == XY) domain->xy = scale[X]*ops[i].ftilt; + else if (ops[i].dim == XZ) domain->xz = scale[X]*ops[i].ftilt; + else if (ops[i].dim == YZ) domain->yz = scale[Y]*ops[i].ftilt; + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + domain->print_box(" "); + + } else if (ops[i].flavor == DELTA) { + if (ops[i].dim == XY) domain->xy += scale[X]*ops[i].dtilt; + else if (ops[i].dim == XZ) domain->xz += scale[X]*ops[i].dtilt; + else if (ops[i].dim == YZ) domain->yz += scale[Y]*ops[i].dtilt; + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + domain->print_box(" "); + } + + } else if (ops[i].style == BOUNDARY) { + domain->set_boundary(3,&arg[ops[i].boundindex],1); + if (domain->dimension == 2 && domain->zperiodic == 0) + error->all(FLERR, + "Cannot run 2d simulation with nonperiodic Z dimension"); + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + + } else if (ops[i].style == ORTHO) { + if (domain->xy != 0.0 || domain->yz != 0.0 || domain->xz != 0.0) + error->all(FLERR, + "Cannot change box to orthogonal when tilt is non-zero"); + if (output->ndump) + error->all(FLERR, + "Cannot change box ortho/triclinic with dumps defined"); + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->no_change_box) + error->all(FLERR, + "Cannot change box ortho/triclinic with " + "certain fixes defined"); + domain->triclinic = 0; + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + domain->print_box(" "); + + } else if (ops[i].style == TRICLINIC) { + if (output->ndump) + error->all(FLERR, + "Cannot change box ortho/triclinic with dumps defined"); + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->no_change_box) + error->all(FLERR, + "Cannot change box ortho/triclinic with " + "certain fixes defined"); + domain->triclinic = 1; + domain->set_lamda_box(); + domain->set_initial_box(); + domain->set_global_box(); + domain->set_local_box(); + domain->print_box(" "); + + } else if (ops[i].style == SET) { + save_box_state(); + + } else if (ops[i].style == REMAP) { + + // convert atoms to lamda coords, using last box state + // convert atoms back to box coords, using current box state + // save current box state + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + domain->x2lamda(x[i],x[i],boxlo,h_inv); + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + domain->lamda2x(x[i],x[i]); + + save_box_state(); + } + } + + // clean up + + delete [] ops; + + // apply shrink-wrap boundary conditions + + if (domain->nonperiodic == 2) domain->reset_box(); + + // move atoms back inside simulation box and to new processors + // use remap() instead of pbc() + // in case box moved a long distance relative to atoms + // use irregular() in case box moved a long distance relative to atoms + + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); + + if (domain->triclinic) domain->x2lamda(atom->nlocal); + domain->reset_box(); + Irregular *irregular = new Irregular(lmp); + irregular->migrate_atoms(); + delete irregular; + if (domain->triclinic) domain->lamda2x(atom->nlocal); + + // check if any atoms were lost + + domain->print_box(" "); + + bigint natoms; + bigint nblocal = atom->nlocal; + MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); + if (natoms != atom->natoms) { + char str[128]; + sprintf(str,"Lost atoms via change_box: original " BIGINT_FORMAT + " current " BIGINT_FORMAT,atom->natoms,natoms); + error->warning(FLERR,str); + } +} + +/* ---------------------------------------------------------------------- + parse optional parameters +------------------------------------------------------------------------- */ + +void ChangeBox::options(int narg, char **arg) +{ + if (narg < 0) error->all(FLERR,"Illegal change_box command"); + + scaleflag = 1; + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal change_box command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all(FLERR,"Illegal change_box command"); + iarg += 2; + } else error->all(FLERR,"Illegal change_box command"); + } +} + +/* ---------------------------------------------------------------------- + save current box state for converting atoms to lamda coords +------------------------------------------------------------------------- */ + +void ChangeBox::save_box_state() +{ + boxlo[0] = domain->boxlo[0]; + boxlo[1] = domain->boxlo[1]; + boxlo[2] = domain->boxlo[2]; + + for (int i = 0; i < 6; i++) + h_inv[i] = domain->h_inv[i]; +} + +/* ---------------------------------------------------------------------- + reset box lengths of dim1/2 to preserve old volume + which changed due to change in 3rd dim +------------------------------------------------------------------------- */ + +void ChangeBox::volume_preserve(int dim1, int dim2, double oldvol) +{ + double newvol; + if (domain->dimension == 2) newvol = domain->xprd * domain->yprd; + else newvol = domain->xprd * domain->yprd * domain->zprd; + + double scale = oldvol/newvol; + double mid,delta; + + // change dim1 only + + if (dim2 < 0) { + mid = 0.5 * (domain->boxlo[dim1] + domain->boxhi[dim1]); + delta = domain->boxlo[dim1] - mid; + domain->boxlo[dim1] = mid + scale*delta; + delta = domain->boxhi[dim1] - mid; + domain->boxhi[dim1] = mid + scale*delta; + + // change dim1 and dim2, keeping their relative aspect constant + // means both are scaled by sqrt(scale) + + } else { + mid = 0.5 * (domain->boxlo[dim1] + domain->boxhi[dim1]); + delta = domain->boxlo[dim1] - mid; + domain->boxlo[dim1] = mid + sqrt(scale)*delta; + delta = domain->boxhi[dim1] - mid; + domain->boxhi[dim1] = mid + sqrt(scale)*delta; + mid = 0.5 * (domain->boxlo[dim2] + domain->boxhi[dim2]); + delta = domain->boxlo[dim2] - mid; + domain->boxlo[dim2] = mid + sqrt(scale)*delta; + delta = domain->boxhi[dim2] - mid; + domain->boxhi[dim2] = mid + sqrt(scale)*delta; + } } diff --git a/src/change_box.h b/src/change_box.h index 9399dcad8e..9b2eac97ae 100644 --- a/src/change_box.h +++ b/src/change_box.h @@ -28,6 +28,28 @@ class ChangeBox : protected Pointers { public: ChangeBox(class LAMMPS *); void command(int, char **); + + private: + int scaleflag; + double scale[3]; + + struct Operation { + int style,flavor; + int dim,boundindex; + int vdim1,vdim2; + double flo,fhi,ftilt; + double dlo,dhi,dtilt; + double scale; + }; + + Operation *ops; + int nops; + + double boxlo[3],h_inv[6]; + + void options(int, char **); + void save_box_state(); + void volume_preserve(int, int, double); }; } @@ -37,7 +59,7 @@ class ChangeBox : protected Pointers { /* ERROR/WARNING messages: -E: Change_box command before simulation box is defined +E: Displace_box command before simulation box is defined Self-explanatory. @@ -47,22 +69,41 @@ Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. -E: Change_box operation is invalid +E: Cannot displace_box after reading restart file with per-atom info -Cannot change orthogonal box to orthogonal or a triclinic box to -triclinic. +This is because the restart file info cannot be migrated with the +atoms. You can get around this by performing a 0-timestep run which +will assign the restart file info to actual atoms. -E: Cannot change box to orthogonal when tilt is non-zero +E: Could not find displace_box group ID -Self-explanatory +Group ID used in the displace_box command does not exist. -E: Cannot change box with dumps defined +E: Displace_box tilt factors require triclinic box + +Cannot use tilt factors unless the simulation box is +non-orthogonal. + +E: Cannot displace_box on a non-periodic boundary Self-explanatory. -E: Cannot change box with certain fixes defined +E: Use of displace_box with undefined lattice -The change_box command cannot be used when fix ave/spatial or -fix/deform are defined . +Must use lattice command with displace_box command if units option is +set to lattice. + +E: Fix deform volume setting is invalid + +Cannot use volume style unless other dimensions are being controlled. + +E: Induced tilt by displace_box is too large + +The final tilt value must be between -1/2 and 1/2 of the perpendicular +box length. + +E: Lost atoms via displace_box: original %ld current %ld + +UNDOCUMENTED */ diff --git a/src/displace_box.cpp b/src/displace_box.cpp deleted file mode 100644 index 79f8865788..0000000000 --- a/src/displace_box.cpp +++ /dev/null @@ -1,418 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "lmptype.h" -#include "mpi.h" -#include "math.h" -#include "stdlib.h" -#include "string.h" -#include "displace_box.h" -#include "atom.h" -#include "modify.h" -#include "domain.h" -#include "lattice.h" -#include "comm.h" -#include "irregular.h" -#include "group.h" -#include "error.h" - -using namespace LAMMPS_NS; - -enum{NONE,FINAL,DELTA,SCALE,VOLUME}; -enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; -enum{NO_REMAP,X_REMAP}; - -/* ---------------------------------------------------------------------- */ - -DisplaceBox::DisplaceBox(LAMMPS *lmp) : Pointers(lmp) {} - -/* ---------------------------------------------------------------------- */ - -void DisplaceBox::command(int narg, char **arg) -{ - int i; - - if (domain->box_exist == 0) - error->all(FLERR,"Displace_box command before simulation box is defined"); - if (narg < 2) error->all(FLERR,"Illegal displace_box command"); - if (modify->nfix_restart_peratom) - error->all(FLERR,"Cannot displace_box after " - "reading restart file with per-atom info"); - - if (comm->me == 0 && screen) fprintf(screen,"Displacing box ...\n"); - - // group - - int igroup = group->find(arg[0]); - if (igroup == -1) error->all(FLERR,"Could not find displace_box group ID"); - int groupbit = group->bitmask[igroup]; - - // set defaults - - set = new Set[6]; - set[0].style = set[1].style = set[2].style = - set[3].style = set[4].style = set[5].style = NONE; - - // parse arguments - - int triclinic = domain->triclinic; - - int index; - int iarg = 1; - while (iarg < narg) { - if (strcmp(arg[iarg],"x") == 0 || strcmp(arg[iarg],"y") == 0 || - strcmp(arg[iarg],"z") == 0) { - if (strcmp(arg[iarg],"x") == 0) index = 0; - else if (strcmp(arg[iarg],"y") == 0) index = 1; - else if (strcmp(arg[iarg],"z") == 0) index = 2; - - if (iarg+2 > narg) error->all(FLERR,"Illegal displace_box command"); - if (strcmp(arg[iarg+1],"final") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal displace_box command"); - set[index].style = FINAL; - set[index].flo = atof(arg[iarg+2]); - set[index].fhi = atof(arg[iarg+3]); - iarg += 4; - } else if (strcmp(arg[iarg+1],"delta") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal displace_box command"); - set[index].style = DELTA; - set[index].dlo = atof(arg[iarg+2]); - set[index].dhi = atof(arg[iarg+3]); - iarg += 4; - } else if (strcmp(arg[iarg+1],"scale") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal displace_box command"); - set[index].style = SCALE; - set[index].scale = atof(arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg+1],"volume") == 0) { - set[index].style = VOLUME; - iarg += 2; - } else error->all(FLERR,"Illegal displace_box command"); - - } else if (strcmp(arg[iarg],"xy") == 0 || strcmp(arg[iarg],"xz") == 0 || - strcmp(arg[iarg],"yz") == 0) { - if (triclinic == 0) - error->all(FLERR,"Displace_box tilt factors require triclinic box"); - if (strcmp(arg[iarg],"xy") == 0) index = 5; - else if (strcmp(arg[iarg],"xz") == 0) index = 4; - else if (strcmp(arg[iarg],"yz") == 0) index = 3; - if (iarg+2 > narg) error->all(FLERR,"Illegal displace_box command"); - if (strcmp(arg[iarg+1],"final") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal displace_box command"); - set[index].style = FINAL; - set[index].ftilt = atof(arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg+1],"delta") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal displace_box command"); - set[index].style = DELTA; - set[index].dtilt = atof(arg[iarg+2]); - iarg += 3; - } else error->all(FLERR,"Illegal displace_box command"); - - } else break; - } - - // read options from end of input line - - options(narg-iarg,&arg[iarg]); - - // check periodicity - - if ((set[0].style && domain->xperiodic == 0) || - (set[1].style && domain->yperiodic == 0) || - (set[2].style && domain->zperiodic == 0)) - error->all(FLERR,"Cannot displace_box on a non-periodic boundary"); - - if (set[3].style && (domain->yperiodic == 0 || domain->zperiodic == 0)) - error->all(FLERR,"Cannot displace_box on a non-periodic boundary"); - if (set[4].style && (domain->xperiodic == 0 || domain->zperiodic == 0)) - error->all(FLERR,"Cannot displace_box on a non-periodic boundary"); - if (set[5].style && (domain->xperiodic == 0 || domain->yperiodic == 0)) - error->all(FLERR,"Cannot displace_box on a non-periodic boundary"); - - // apply scaling to FINAL,DELTA since they have distance units - - int flag = 0; - for (int i = 0; i < 6; i++) - if (set[i].style == FINAL || set[i].style == DELTA) flag = 1; - - if (flag && scaleflag && domain->lattice == NULL) - error->all(FLERR,"Use of displace_box with undefined lattice"); - - double xscale,yscale,zscale; - if (flag && scaleflag) { - xscale = domain->lattice->xlattice; - yscale = domain->lattice->ylattice; - zscale = domain->lattice->zlattice; - } - else xscale = yscale = zscale = 1.0; - - // for 3,4,5 scaling is in 1st dimension, e.g. x for xz - - double map[6]; - map[0] = xscale; map[1] = yscale; map[2] = zscale; - map[3] = yscale; map[4] = xscale; map[5] = xscale; - - for (int i = 0; i < 3; i++) { - if (set[i].style == FINAL) { - set[i].flo *= map[i]; - set[i].fhi *= map[i]; - } else if (set[i].style == DELTA) { - set[i].dlo *= map[i]; - set[i].dhi *= map[i]; - } - } - - for (int i = 3; i < 6; i++) { - if (set[i].style == FINAL) set[i].ftilt *= map[i]; - else if (set[i].style == DELTA) set[i].dtilt *= map[i]; - } - - // set initial/final values for box size and shape - // final = initial if no setting - - for (int i = 0; i < 3; i++) { - set[i].lo_stop = set[i].lo_start = domain->boxlo[i]; - set[i].hi_stop = set[i].hi_start = domain->boxhi[i]; - - if (set[i].style == FINAL) { - set[i].lo_stop = set[i].flo; - set[i].hi_stop = set[i].fhi; - } else if (set[i].style == DELTA) { - set[i].lo_stop = set[i].lo_start + set[i].dlo; - set[i].hi_stop = set[i].hi_start + set[i].dhi; - } else if (set[i].style == SCALE) { - set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); - set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); - } - } - - for (int i = 3; i < 6; i++) { - if (i == 5) set[i].tilt_start = domain->xy; - else if (i == 4) set[i].tilt_start = domain->xz; - else if (i == 3) set[i].tilt_start = domain->yz; - set[i].tilt_stop = set[i].tilt_start; - - if (set[i].style == FINAL) { - set[i].tilt_stop = set[i].ftilt; - } else if (set[i].style == DELTA) { - set[i].tilt_stop = set[i].tilt_start + set[i].dtilt; - } - } - - // for VOLUME, setup links to other dims - // fixed, dynamic1,2, vol_start - - for (int i = 0; i < 3; i++) { - set[i].vol_start = domain->xprd * domain->yprd * domain->zprd; - - if (set[i].style != VOLUME) continue; - int other1 = (i+1) % 3; - int other2 = (i+2) % 3; - - if (set[other1].style == NONE) { - if (set[other2].style == NONE || set[other2].style == VOLUME) - error->all(FLERR,"Fix deform volume setting is invalid"); - set[i].substyle = ONE_FROM_ONE; - set[i].fixed = other1; - set[i].dynamic1 = other2; - } else if (set[other2].style == NONE) { - if (set[other1].style == NONE || set[other1].style == VOLUME) - error->all(FLERR,"Fix deform volume setting is invalid"); - set[i].substyle = ONE_FROM_ONE; - set[i].fixed = other2; - set[i].dynamic1 = other1; - } else if (set[other1].style == VOLUME) { - if (set[other2].style == NONE || set[other2].style == VOLUME) - error->all(FLERR,"Fix deform volume setting is invalid"); - set[i].substyle = TWO_FROM_ONE; - set[i].fixed = other1; - set[i].dynamic1 = other2; - } else if (set[other2].style == VOLUME) { - if (set[other1].style == NONE || set[other1].style == VOLUME) - error->all(FLERR,"Fix deform volume setting is invalid"); - set[i].substyle = TWO_FROM_ONE; - set[i].fixed = other2; - set[i].dynamic1 = other1; - } else { - set[i].substyle = ONE_FROM_TWO; - set[i].dynamic2 = other1; - set[i].dynamic2 = other2; - } - } - - // set new box size for VOLUME dims that are linked to other dims - - for (int i = 0; i < 3; i++) { - if (set[i].style != VOLUME) continue; - - if (set[i].substyle == ONE_FROM_ONE) { - set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_stop - - set[set[i].dynamic1].lo_stop) / - (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); - set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_stop - - set[set[i].dynamic1].lo_stop) / - (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); - - } else if (set[i].substyle == ONE_FROM_TWO) { - set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_stop - - set[set[i].dynamic1].lo_stop) / - (set[set[i].dynamic2].hi_stop - - set[set[i].dynamic2].lo_stop)); - set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_stop - - set[set[i].dynamic1].lo_stop) / - (set[set[i].dynamic2].hi_stop - - set[set[i].dynamic2].lo_stop)); - - } else if (set[i].substyle == TWO_FROM_ONE) { - set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*sqrt(set[i].vol_start / - (set[set[i].dynamic1].hi_stop - - set[set[i].dynamic1].lo_stop) / - (set[set[i].fixed].hi_start - - set[set[i].fixed].lo_start) * - (set[i].hi_start - set[i].lo_start)); - set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*sqrt(set[i].vol_start / - (set[set[i].dynamic1].hi_stop - - set[set[i].dynamic1].lo_stop) / - (set[set[i].fixed].hi_start - - set[set[i].fixed].lo_start) * - (set[i].hi_start - set[i].lo_start)); - } - } - - // check that final tilt is not illegal value - - double xprd_stop = set[0].hi_stop - set[0].lo_stop; - double yprd_stop = set[0].hi_stop - set[0].lo_stop; - - if (set[3].tilt_stop < -0.5*yprd_stop || set[3].tilt_stop > 0.5*yprd_stop || - set[4].tilt_stop < -0.5*xprd_stop || set[4].tilt_stop > 0.5*xprd_stop || - set[5].tilt_stop < -0.5*xprd_stop || set[5].tilt_stop > 0.5*xprd_stop) - error->all(FLERR,"Induced tilt by displace_box is too large"); - - // convert atoms to lamda coords - - if (remapflag == X_REMAP) { - double **x = atom->x; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - domain->x2lamda(x[i],x[i]); - } - - // reset global and local box to new size/shape - - domain->boxlo[0] = set[0].lo_stop; - domain->boxlo[1] = set[1].lo_stop; - domain->boxlo[2] = set[2].lo_stop; - domain->boxhi[0] = set[0].hi_stop; - domain->boxhi[1] = set[1].hi_stop; - domain->boxhi[2] = set[2].hi_stop; - - if (triclinic) { - domain->yz = set[3].tilt_stop; - domain->xz = set[4].tilt_stop; - domain->xy = set[5].tilt_stop; - } - - domain->set_global_box(); - domain->set_local_box(); - - // convert atoms back to box coords - - if (remapflag == X_REMAP) { - double **x = atom->x; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - domain->lamda2x(x[i],x[i]); - } - - // move atoms back inside simulation box and to new processors - // use remap() instead of pbc() - // in case box moved a long distance relative to atoms - // use irregular() in case box moved a long distance relative to atoms - - double **x = atom->x; - int *image = atom->image; - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); - - if (domain->triclinic) domain->x2lamda(atom->nlocal); - domain->reset_box(); - Irregular *irregular = new Irregular(lmp); - irregular->migrate_atoms(); - delete irregular; - if (domain->triclinic) domain->lamda2x(atom->nlocal); - - // clean up - - delete [] set; - - // check if any atoms were lost - - bigint natoms; - bigint nblocal = atom->nlocal; - MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); - if (natoms != atom->natoms) { - char str[128]; - sprintf(str,"Lost atoms via displace_box: original " BIGINT_FORMAT - " current " BIGINT_FORMAT,atom->natoms,natoms); - error->all(FLERR,str); - } -} - -/* ---------------------------------------------------------------------- - parse optional parameters at end of displace_box input line -------------------------------------------------------------------------- */ - -void DisplaceBox::options(int narg, char **arg) -{ - if (narg < 0) error->all(FLERR,"Illegal displace_box command"); - - remapflag = X_REMAP; - scaleflag = 1; - - int iarg = 0; - while (iarg < narg) { - if (strcmp(arg[iarg],"remap") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal displace_box command"); - if (strcmp(arg[iarg+1],"x") == 0) remapflag = X_REMAP; - else if (strcmp(arg[iarg+1],"none") == 0) remapflag = NO_REMAP; - else error->all(FLERR,"Illegal displace_box command"); - iarg += 2; - } else if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal displace_box command"); - if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; - else error->all(FLERR,"Illegal displace_box command"); - iarg += 2; - } else error->all(FLERR,"Illegal displace_box command"); - } -} diff --git a/src/displace_box.h b/src/displace_box.h deleted file mode 100644 index 44df201c36..0000000000 --- a/src/displace_box.h +++ /dev/null @@ -1,105 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef COMMAND_CLASS - -CommandStyle(displace_box,DisplaceBox) - -#else - -#ifndef LMP_DISPLACE_BOX_H -#define LMP_DISPLACE_BOX_H - -#include "pointers.h" - -namespace LAMMPS_NS { - -class DisplaceBox : protected Pointers { - public: - DisplaceBox(class LAMMPS *); - void command(int, char **); - - private: - int remapflag,scaleflag; - - struct Set { - int style,substyle; - double flo,fhi,ftilt; - double dlo,dhi,dtilt; - double scale; - double lo_start,hi_start; - double lo_stop,hi_stop; - double tilt_start,tilt_stop; - double vol_start; - int fixed,dynamic1,dynamic2; - }; - Set *set; - - void options(int, char **); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Displace_box command before simulation box is defined - -Self-explanatory. - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Cannot displace_box after reading restart file with per-atom info - -This is because the restart file info cannot be migrated with the -atoms. You can get around this by performing a 0-timestep run which -will assign the restart file info to actual atoms. - -E: Could not find displace_box group ID - -Group ID used in the displace_box command does not exist. - -E: Displace_box tilt factors require triclinic box - -Cannot use tilt factors unless the simulation box is -non-orthogonal. - -E: Cannot displace_box on a non-periodic boundary - -Self-explanatory. - -E: Use of displace_box with undefined lattice - -Must use lattice command with displace_box command if units option is -set to lattice. - -E: Fix deform volume setting is invalid - -Cannot use volume style unless other dimensions are being controlled. - -E: Induced tilt by displace_box is too large - -The final tilt value must be between -1/2 and 1/2 of the perpendicular -box length. - -E: Lost atoms via displace_box: original %ld current %ld - -UNDOCUMENTED - -*/ diff --git a/src/domain.cpp b/src/domain.cpp index c729ac16bb..a83f958f6a 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1109,10 +1109,12 @@ int Domain::find_region(char *name) } /* ---------------------------------------------------------------------- - boundary settings from the input script + (re)set boundary settings + flag = 0, called from the input script + flag = 1, called from change box command ------------------------------------------------------------------------- */ -void Domain::set_boundary(int narg, char **arg) +void Domain::set_boundary(int narg, char **arg, int flag) { if (narg != 3) error->all(FLERR,"Illegal boundary command"); @@ -1127,7 +1129,10 @@ void Domain::set_boundary(int narg, char **arg) else if (c == 'f') boundary[idim][iside] = 1; else if (c == 's') boundary[idim][iside] = 2; else if (c == 'm') boundary[idim][iside] = 3; - else error->all(FLERR,"Illegal boundary command"); + else { + if (flag == 0) error->all(FLERR,"Illegal boundary command"); + if (flag == 1) error->all(FLERR,"Illegal change_box command"); + } } for (int idim = 0; idim < 3; idim++) @@ -1257,6 +1262,26 @@ void Domain::x2lamda(double *x, double *lamda) lamda[2] = h_inv[2]*delta[2]; } +/* ---------------------------------------------------------------------- + convert box coords to triclinic 0-1 lamda coords for one atom + use my_boxlo & my_h_inv stored by caller for previous state of box + lamda = H^-1 (x - x0) + x and lamda can point to same 3-vector +------------------------------------------------------------------------- */ + +void Domain::x2lamda(double *x, double *lamda, + double *my_boxlo, double *my_h_inv) +{ + double delta[3]; + delta[0] = x[0] - my_boxlo[0]; + delta[1] = x[1] - my_boxlo[1]; + delta[2] = x[2] - my_boxlo[2]; + + lamda[0] = my_h_inv[0]*delta[0] + my_h_inv[5]*delta[1] + my_h_inv[4]*delta[2]; + lamda[1] = my_h_inv[1]*delta[1] + my_h_inv[3]*delta[2]; + lamda[2] = my_h_inv[2]*delta[2]; +} + /* ---------------------------------------------------------------------- convert 8 lamda corner pts of lo/hi box to box coords return bboxlo/hi = bounding box around 8 corner pts in box coords diff --git a/src/domain.h b/src/domain.h index 0496783c1f..d21119cffb 100644 --- a/src/domain.h +++ b/src/domain.h @@ -107,13 +107,14 @@ class Domain : protected Pointers { void add_region(int, char **); void delete_region(int, char **); int find_region(char *); - void set_boundary(int, char **); + void set_boundary(int, char **, int); void print_box(const char *); virtual void lamda2x(int); virtual void x2lamda(int); void lamda2x(double *, double *); void x2lamda(double *, double *); + void x2lamda(double *, double *, double *, double *); void bbox(double *, double *, double *, double *); void box_corners(); diff --git a/src/input.cpp b/src/input.cpp index 865975c68c..9c1905168e 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -930,7 +930,7 @@ void Input::boundary() { if (domain->box_exist) error->all(FLERR,"Boundary command after simulation box is defined"); - domain->set_boundary(narg,arg); + domain->set_boundary(narg,arg,0); } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 2ffd9c096e..73dd0709d5 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -16,6 +16,7 @@ #include "string.h" #include "ctype.h" #include "pair_hybrid.h" +#include "style_pair.h" #include "atom.h" #include "force.h" #include "pair.h" @@ -763,4 +764,3 @@ double PairHybrid::memory_usage() for (int m = 0; m < nstyles; m++) bytes += styles[m]->memory_usage(); return bytes; } - diff --git a/src/read_restart.cpp b/src/read_restart.cpp index 02e313a8d0..c015449ff4 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -464,7 +464,8 @@ void ReadRestart::header() if (flag == VERSION) { char *version = read_char(); if (strcmp(version,universe->version) != 0 && me == 0) { - error->warning(FLERR,"Restart file version does not match LAMMPS version"); + error->warning(FLERR, + "Restart file version does not match LAMMPS version"); if (screen) fprintf(screen," restart file = %s, LAMMPS = %s\n", version,universe->version); } @@ -581,7 +582,8 @@ void ReadRestart::header() boundary[2][0] != domain->boundary[2][0] || boundary[2][1] != domain->boundary[2][1]) { if (me == 0) - error->warning(FLERR,"Restart file used different boundary settings, " + error->warning(FLERR, + "Restart file used different boundary settings, " "using restart file values"); } }