From 1aefa8d62472b4a2fae18a57315aba640861fca1 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 5 Jan 2015 15:41:14 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12859 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/fix_srp.cpp | 58 ++++++++++++++++++++++-------------- src/USER-MISC/pair_srp.cpp | 61 ++++++++++++++++++++------------------ 2 files changed, 68 insertions(+), 51 deletions(-) diff --git a/src/USER-MISC/fix_srp.cpp b/src/USER-MISC/fix_srp.cpp index 1b22c8945f..c213edf2d3 100644 --- a/src/USER-MISC/fix_srp.cpp +++ b/src/USER-MISC/fix_srp.cpp @@ -26,6 +26,7 @@ #include "error.h" #include "neighbor.h" #include "atom_vec.h" +#include "modify.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -66,11 +67,11 @@ FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) array[i][m] = 0.0; // assume setup of fix is needed to insert particles - // yes if reading datafile - // maybe if reading restart - // a datafile cannot contain BPs - // a restart file written during the run has BPs as per atom data - // a restart file written after the run does not have BPs + // is true if reading from datafile + // since a datafile cannot contain bond particles + // might be true if reading from restart + // a restart file written during the run has bond particles as per atom data + // a restart file written after the run does not have bond particles setup = 1; } @@ -80,7 +81,7 @@ FixSRP::~FixSRP() { // unregister callbacks to this fix from Atom class atom->delete_callback(id,0); - atom->delete_callback(id,1); // for restart + atom->delete_callback(id,1); atom->delete_callback(id,2); memory->destroy(array); } @@ -104,18 +105,16 @@ void FixSRP::init() if (force->pair_match("hybrid",1) == NULL) error->all(FLERR,"Cannot use pair srp without pair_style hybrid"); - // reserve the highest numbered atom type for bond particles (BPs) - // set bptype if not read from restart + // the highest numbered atom type should be reserved for bond particles (bptype) + // set bptype, unless it will be read from restart if(!restart_reset) bptype = atom->ntypes; - // check for prescense of extra atom type for bond particle - // if reading a rst file written in the middle of run - // then bond particles are already present - // otherwise need to insert particles in setup + // check if bptype is already in use for(int i=0; i < atom->nlocal; i++) if(atom->type[i] == bptype){ - // error if missing extra atom type in either - // data file or rst file without this fix + // error if bptype is already in use + // unless starting from a rst file + // since rst can contain the bond particles as per atom data if(!restart_reset) error->all(FLERR,"Fix SRP requires an extra atom type"); else @@ -123,7 +122,7 @@ void FixSRP::init() } // setup neigh exclusions for diff atom types - // BPs do not interact with other types + // bond particles do not interact with other types // type bptype only interacts with itself char* arg1[4]; arg1[0] = (char *) "exclude"; @@ -155,8 +154,8 @@ void FixSRP::setup_pre_force(int zz) bigint nall = atom->nlocal + atom->nghost; double xold[nall][3]; - // make a copy of coordinates and tags - // create_atom overwrites ghost atoms + // make a copy of all coordinates and tags + // need this because create_atom overwrites ghost atoms for(int i = 0; i < nall; i++){ xold[i][0] = x[i][0]; xold[i][1] = x[i][1]; @@ -189,7 +188,8 @@ void FixSRP::setup_pre_force(int zz) xone[1] = (xold[i][1] + xold[j][1])*0.5; xone[2] = (xold[i][2] + xold[j][2])*0.5; - // record longest bond for ghostcut + // record longest bond + // this used to set ghost cutoff delx = xold[j][0] - xold[i][0]; dely = xold[j][1] - xold[i][1]; delz = xold[j][2] - xold[i][2]; @@ -262,11 +262,20 @@ void FixSRP::setup_pre_force(int zz) // move owned to new procs // get ghosts // build neigh lists again + + // if triclinic, lambda coords needed for pbc, exchange, borders + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); comm->setup(); + if (neighbor->style) neighbor->setup_bins(); comm->exchange(); if (atom->sortfreq > 0) atom->sort(); comm->borders(); + // back to box coords + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + domain->image_check(); + domain->box_too_small_check(); + modify->setup_pre_neighbor(); neighbor->build(); neighbor->ncalls = 0; @@ -295,7 +304,7 @@ void FixSRP::pre_exchange() comm->forward_comm(); // reassign bond particle coordinates to midpoint of bonds - // only needed before neigh rebuild + // only need to do this before neigh rebuild double **x=atom->x; int i,j; int nlocal = atom->nlocal; @@ -422,7 +431,7 @@ void FixSRP::post_run() { // all bond particles are removed after each run // useful for write_data and write_restart commands - // since it occurs between runs + // since those commands occur between runs bigint natoms_previous = atom->natoms; int nlocal = atom->nlocal; @@ -481,11 +490,16 @@ void FixSRP::post_run() // verlet calls box_too_small_check() in post_run // this check maps all bond partners // therefore need ghosts + + // need to convert to lambda coords before apply pbc + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); comm->setup(); comm->exchange(); if (atom->sortfreq > 0) atom->sort(); comm->borders(); + // change back to box coordinates + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); } /* ---------------------------------------------------------------------- @@ -572,8 +586,8 @@ void FixSRP::restart(char *buf) } /* ---------------------------------------------------------------------- - interface with other classes - pair srp sets the bond type for this fix + interface with pair class + pair srp sets the bond type in this fix ------------------------------------------------------------------------- */ int FixSRP::modify_param(int narg, char **arg) diff --git a/src/USER-MISC/pair_srp.cpp b/src/USER-MISC/pair_srp.cpp index 163949a247..f32dc96e01 100644 --- a/src/USER-MISC/pair_srp.cpp +++ b/src/USER-MISC/pair_srp.cpp @@ -45,7 +45,7 @@ Please contact Timothy Sirk for questions (tim.sirk@us.army.mil). using namespace LAMMPS_NS; -#define SMALL 1.0e-3 +#define SMALL 1.0e-10 #define BIG 1e10 #define ONETWOBIT 0x40000000 @@ -155,8 +155,9 @@ void PairSRP::compute(int eflag, int vflag) } // this pair style only used with hybrid - // so each atom i is type bptype - // each neigh j is type bptype due to exclusions + // due to exclusions + // each atom i is type bptype + // each neigh j is type bptype // using midpoint distance option if(midpoint){ @@ -179,7 +180,7 @@ void PairSRP::compute(int eflag, int vflag) continue; j &= NEIGHMASK; - //atoms inside bond particle + //retrieve atoms from bond particle i1 = segment[j][0]; j1 = segment[j][1]; @@ -190,13 +191,13 @@ void PairSRP::compute(int eflag, int vflag) dijsq = dx*dx + dy*dy + dz*dz; if (dijsq < cutsq[bptype][bptype]){ - dij = pow(dijsq, 0.5); + dij = sqrt(dijsq); - if (dij < SMALL) - dij = SMALL; // prevent explosions + if (dij < SMALL) + continue; // dij can be 0.0 with soft potentials wd = 1.0 - dij / cut[bptype][bptype]; - fpair = 0.5 * a0[bptype][bptype] * wd / dij; + fpair = 0.5 * a0[bptype][bptype] * wd / dij; // 0.5 factor for lever rule // force for bond 0, beads 0,1 //force between bonds @@ -266,13 +267,13 @@ void PairSRP::compute(int eflag, int vflag) if (dijsq < cutsq[bptype][bptype]){ - dij = pow(dijsq, 0.5); + dij = sqrt(dijsq); - if (dij < SMALL) - dij = SMALL; // prevent explosions + if (dij < SMALL) + continue; // dij can be 0.0 with soft potentials wd = 1.0 - dij / cut[bptype][bptype]; - fpair = a0[bptype][bptype] * wd / dij; + fpair = a0[bptype][bptype] * wd / dij; // force for bond 0, beads 0,1 lever0 = 0.5 + ti; // assign force according to lever rule @@ -281,6 +282,7 @@ void PairSRP::compute(int eflag, int vflag) fx = fpair * dx; fy = fpair * dy; fz = fpair * dz; + //decompose onto atoms fxlever0 = fx * lever0; fylever0 = fy * lever0; @@ -353,16 +355,18 @@ void PairSRP::settings(int narg, char **arg) while (iarg < narg) { if (strcmp(arg[iarg],"exclude") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp exclude command"); + if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command"); if (strcmp(arg[iarg+1],"yes") == 0) - exclude = 1; - if (strcmp(arg[iarg+1],"no") == 0) - exclude = 0; + exclude = 1; + if (strcmp(arg[iarg+1],"no") == 0){ + if (min) error->all(FLERR,"Illegal exclude option in pair srp command"); + exclude = 0; + } iarg += 2; } else error->all(FLERR,"Illegal pair srp command"); } - // highest atom type is saved for bond particles (BPs) + // highest atom type is saved for bond particles bptype = atom->ntypes; // reset cutoffs if explicitly set @@ -372,7 +376,6 @@ void PairSRP::settings(int narg, char **arg) for (j = i+1; j <= bptype; j++) if (setflag[i][j]) cut[i][j] = cut_global; } - } /* ---------------------------------------------------------------------- @@ -422,16 +425,16 @@ void PairSRP::init_style() // need fix srp // invoke here instead of constructor // to make restart possible - char **newarg = new char*[3]; - newarg[0] = (char *) "pair_srp"; - newarg[1] = (char *) "all"; - newarg[2] = (char *) "SRP"; - modify->add_fix(3,newarg); + char **fixarg = new char*[3]; + fixarg[0] = (char *) "mysrp"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "SRP"; + modify->add_fix(3,fixarg); f_srp = (FixSRP *) modify->fix[modify->nfix-1]; - delete [] newarg; + delete [] fixarg; // set bond type in fix srp - // bonds of this type will be represented by bond particles (BPs) + // bonds of this type will be represented by bond particles // btype = bond type // bptype = bond particle type char c0[20]; @@ -441,10 +444,10 @@ void PairSRP::init_style() arg0[1] = c0; f_srp->modify_params(2, arg0); - // BPs do not contribute to energy or virial - // BPs do not belong to group all - // but thermo normalization is by nall - // turn off normalization + // bond particles do not contribute to energy or virial + // bond particles do not belong to group all + // but thermo normalization is by nall + // therefore should turn off normalization int me; MPI_Comm_rank(world,&me); char *arg1[2];