jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12859 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-01-05 15:41:14 +00:00
parent 32f47e4a99
commit 1aefa8d624
2 changed files with 68 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
src/USER-MISC/fix_srp.cpp
View File

@ -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)

49
src/USER-MISC/pair_srp.cpp
View File

@ -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
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,10 +267,10 @@ 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
continue; // dij can be 0.0 with soft potentials
wd = 1.0 - dij / cut[bptype][bptype];
fpair = a0[bptype][bptype] * wd / dij;
@ -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)
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
// bond particles do not contribute to energy or virial
// bond particles do not belong to group all
// but thermo normalization is by nall
// turn off normalization
// therefore should turn off normalization
int me;
MPI_Comm_rank(world,&me);
char *arg1[2];