Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78637100
pair_srp.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Aug 22, 02:46
Size
21 KB
Mime Type
text/x-c
Expires
Sat, Aug 24, 02:46 (2 d)
Engine
blob
Format
Raw Data
Handle
18459561
Attached To
rLAMMPS lammps
pair_srp.cpp
View Options
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Timothy Sirk (ARL), Pieter in't Veld (BASF)
This pair style srp command calculates a segmental repulsive force
between bonds. This is useful for preventing the crossing of bonds if
soft non-bonded potentials are used, such as DPD polymer chains.
See the doc page for pair_style srp command for usage instructions.
There is an example script for this package in examples/USER/srp.
Please contact Timothy Sirk for questions (tim.sirk@us.army.mil).
------------------------------------------------------------------------- */
#include <stdlib.h>
#include "pair_srp.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "fix_srp.h"
#include "thermo.h"
#include "output.h"
#include <string.h>
#include "citeme.h"
using namespace LAMMPS_NS;
#define SMALL 1.0e-10
#define BIG 1e10
#define ONETWOBIT 0x40000000
static const char cite_srp[] =
"@Article{Sirk2012\n"
" author = {T. Sirk and Y. Sliozberg and J. Brennan and M. Lisal and J. Andzelm},\n"
" title = {An enhanced entangled polymer model for dissipative particle dynamics},\n"
" journal = {J.~Chem.~Phys.},\n"
" year = 2012,\n"
" volume = 136,\n"
" pages = {134903}\n"
"}\n\n";
static int srp_instance = 0;
/* ----------------------------------------------------------------------
set size of pair comms in constructor
---------------------------------------------------------------------- */
PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp)
{
writedata = 1;
if (lmp->citeme) lmp->citeme->add(cite_srp);
nextra = 1;
segment = NULL;
// generate unique fix-id for this pair style instance
fix_id = strdup("XX_FIX_SRP");
fix_id[0] = '0' + srp_instance / 10;
fix_id[1] = '0' + srp_instance % 10;
++srp_instance;
// create fix SRP instance here, as it has to
// be executed before all other fixes
char **fixarg = new char*[3];
fixarg[0] = fix_id;
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "SRP";
modify->add_fix(3,fixarg);
f_srp = (FixSRP *) modify->fix[modify->nfix-1];
delete [] fixarg;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSRP::allocate()
{
allocated = 1;
// particles of bptype inserted by fix srp
// bptype is the highest numbered atom type
int n = bptype;
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(cut, n + 1, n + 1, "pair:cut");
memory->create(a0, n + 1, n + 1, "pair:a0");
// setflag for atom types
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
maxcount = 0;
}
/* ----------------------------------------------------------------------
free
------------------------------------------------------------------------- */
PairSRP::~PairSRP()
{
if (allocated)
{
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(a0);
memory->destroy(segment);
}
// check nfix in case all fixes have already been deleted
if (modify->nfix) modify->delete_fix(fix_id);
free(fix_id);
}
/* ----------------------------------------------------------------------
compute bond-bond repulsions
------------------------------------------------------------------------- */
void PairSRP::compute(int eflag, int vflag)
{
// setup energy and virial
if (eflag || vflag)
ev_setup(eflag, vflag);
else
evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int i0, i1, j0, j1;
int i,j,ii,jj,inum,jnum;
double dijsq, dij;
int *ilist,*jlist,*numneigh,**firstneigh;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double dx,dy,dz,ti,tj;
double wd, lever0, lever1, evdwl, fpair;
double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1;
double fx, fy, fz;
evdwl = 0.0;
// mapping global to local for atoms inside bond particles
// exclude 1-2 neighs if requested
if (neighbor->ago == 0){
remapBonds(nall);
if(exclude) onetwoexclude(ilist, inum, jlist, numneigh, firstneigh);
}
// this pair style only used with hybrid
// due to exclusions
// each atom i is type bptype
// each neigh j is type bptype
// using midpoint distance option
if(midpoint){
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jnum = numneigh[i];
// two atoms inside bond particle
i0 = segment[i][0];
j0 = segment[i][1];
for (jj = 0; jj < jnum; jj++) {
jlist = firstneigh[i];
j = jlist[jj];
// enforce 1-2 exclusions
if( (sbmask(j) & exclude) )
continue;
j &= NEIGHMASK;
//retrieve atoms from bond particle
i1 = segment[j][0];
j1 = segment[j][1];
// midpt dist bond 0 and 1
dx = 0.5*(x[i0][0] - x[i1][0] + x[j0][0] - x[j1][0]);
dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]);
dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]);
dijsq = dx*dx + dy*dy + dz*dz;
if (dijsq < cutsq[bptype][bptype]){
dij = sqrt(dijsq);
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; // 0.5 factor for lever rule
// force for bond 0, beads 0,1
//force between bonds
fx = fpair * dx;
fy = fpair * dy;
fz = fpair * dz;
f[i0][0] += fx; //keep force sign for bond 0
f[i0][1] += fy;
f[i0][2] += fz;
f[j0][0] += fx;
f[j0][1] += fy;
f[j0][2] += fz;
f[i1][0] -= fx; //flip force sign for bond 1
f[i1][1] -= fy;
f[i1][2] -= fz;
f[j1][0] -= fx;
f[j1][1] -= fy;
f[j1][2] -= fz;
// ************************************************* //
if (eflag){
evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd;
}
if (evflag){
ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz);
ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz);
}
if (vflag_fdotr) virial_fdotr_compute();
}
}
}
}
else{
// using min distance option
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jnum = numneigh[i];
i0 = segment[i][0];
j0 = segment[i][1];
for (jj = 0; jj < jnum; jj++) {
jlist = firstneigh[i];
j = jlist[jj];
// enforce 1-2 exclusions
if( (sbmask(j) & exclude) )
continue;
j &= NEIGHMASK;
i1 = segment[j][0];
j1 = segment[j][1];
getMinDist(x, dx, dy, dz, ti, tj, i0, j0, i1, j1);
dijsq = dx*dx + dy*dy + dz*dz;
if (dijsq < cutsq[bptype][bptype]){
dij = sqrt(dijsq);
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;
// force for bond 0, beads 0,1
lever0 = 0.5 + ti; // assign force according to lever rule
lever1 = 0.5 + tj; // assign force according to lever rule
//force between bonds
fx = fpair * dx;
fy = fpair * dy;
fz = fpair * dz;
//decompose onto atoms
fxlever0 = fx * lever0;
fylever0 = fy * lever0;
fzlever0 = fz * lever0;
fxlever1 = fx * lever1;
fylever1 = fy * lever1;
fzlever1 = fz * lever1;
f[i0][0] += fxlever0; //keep force sign for bond 0
f[i0][1] += fylever0;
f[i0][2] += fzlever0;
f[j0][0] += (fx - fxlever0);
f[j0][1] += (fy - fylever0);
f[j0][2] += (fz - fzlever0);
f[i1][0] -= fxlever1; //flip force sign for bond 1
f[i1][1] -= fylever1;
f[i1][2] -= fzlever1;
f[j1][0] -= (fx - fxlever1);
f[j1][1] -= (fy - fylever1);
f[j1][2] -= (fz - fzlever1);
// ************************************************* //
if (eflag){
evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd;
}
if (evflag){
ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz);
ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz);
}
if (vflag_fdotr) virial_fdotr_compute();
}
}
}
}
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSRP::settings(int narg, char **arg)
{
if (narg < 3 || narg > 7)
error->all(FLERR,"Illegal pair_style command");
if (atom->tag_enable == 0)
error->all(FLERR,"Pair_style srp requires atom IDs");
cut_global = force->numeric(FLERR,arg[0]);
// wildcard
if (strcmp(arg[1],"*") == 0)
btype = 0;
else {
btype = force->inumeric(FLERR,arg[1]);
if ((btype > atom->nbondtypes) || (btype <= 0))
error->all(FLERR,"Illegal pair_style command");
}
// settings
midpoint = 0;
min = 0;
if (strcmp(arg[2],"min") == 0) min = 1;
else if (strcmp(arg[2],"mid") == 0) midpoint = 1;
else
error->all(FLERR,"Illegal pair_style command");
int iarg = 3;
// default exclude 1-2
// scaling for 1-2, etc not supported
exclude = 1;
// use last atom type by default for bond particles
bptype = atom->ntypes;
while (iarg < narg) {
if (strcmp(arg[iarg],"exclude") == 0) {
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 (min) error->all(FLERR,"Illegal exclude option in pair srp command");
exclude = 0;
}
iarg += 2;
} else if (strcmp(arg[iarg],"bptype") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command");
bptype = force->inumeric(FLERR,arg[iarg+1]);
if ((bptype < 1) || (bptype > atom->ntypes))
error->all(FLERR,"Illegal bond particle type for srp");
iarg += 2;
} else error->all(FLERR,"Illegal pair srp command");
}
// reset cutoffs if explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= bptype; i++)
for (j = i+1; j <= bptype; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ----------------------------------------------------------------------
set coeffs
------------------------------------------------------------------------- */
void PairSRP::coeff(int narg, char **arg)
{
if (narg < 3 || narg > 4)
error->all(FLERR,"PairSRP: Incorrect args for pair coeff");
if (!allocated) allocate();
// set ij bond-bond cutoffs
int ilo, ihi, jlo, jhi;
force->bounds(arg[0], bptype, ilo, ihi);
force->bounds(arg[1], bptype, jlo, jhi);
double a0_one = force->numeric(FLERR,arg[2]);
double cut_one = cut_global;
if (narg == 4) cut_one = force->numeric(FLERR,arg[3]);
int count = 0;
for (int i = ilo; i <= ihi; i++)
{
for (int j = MAX(jlo,i); j <= jhi; j++)
{
a0[i][j] = a0_one;
cut[i][j] = cut_one;
cutsq[i][j] = cut_one * cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->warning(FLERR,"PairSRP: No pair coefficients were set");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairSRP::init_style()
{
if (!force->newton_pair)
error->all(FLERR,"PairSRP: Pair srp requires newton pair on");
// verify that fix SRP is still defined and has not been changed.
int ifix = modify->find_fix(fix_id);
if (f_srp != (FixSRP *)modify->fix[ifix])
error->all(FLERR,"Fix SRP has been changed unexpectedly");
if (comm->me == 0) {
if (screen) fprintf(screen,"Using type %d for bond particles\n",bptype);
if (logfile) fprintf(logfile,"Using type %d for bond particles\n",bptype);
}
// set bond and bond particle types in fix srp
// bonds of this type will be represented by bond particles
// if bond type is 0, then all bonds have bond particles
// btype = bond type
char c0[20];
char* arg0[2];
sprintf(c0, "%d", btype);
arg0[0] = (char *) "btype";
arg0[1] = c0;
f_srp->modify_params(2, arg0);
// bptype = bond particle type
sprintf(c0, "%d", bptype);
arg0[0] = (char *) "bptype";
arg0[1] = c0;
f_srp->modify_params(2, arg0);
// 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];
arg1[0] = (char *) "norm";
arg1[1] = (char *) "no";
output->thermo->modify_params(2, arg1);
if (me == 0)
error->message(FLERR,"Thermo normalization turned off by pair srp");
neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSRP::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"PairSRP: All pair coeffs are not set");
cut[j][i] = cut[i][j];
a0[j][i] = a0[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
find min distance for bonds i0/j0 and i1/j1
------------------------------------------------------------------------- */
inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, double &ti, double &tj, int &i0, int &j0, int &i1, int &j1)
{
// move these outside the loop
double diffx0, diffy0, diffz0, diffx1, diffy1, diffz1, dPx, dPy, dPz, RiRi, RiRj, RjRj;
double denom, termx0, termy0, termz0, num0, termx1, termy1, termz1, num1;
// compute midpt dist from 1st atom, 1st bond
diffx0 = x[j0][0] - x[i0][0]; // x,y,z from bond 0
diffy0 = x[j0][1] - x[i0][1];
diffz0 = x[j0][2] - x[i0][2];
// compute midpt dist from 1st atom, 2nd bond
diffx1 = x[j1][0] - x[i1][0];
diffy1 = x[j1][1] - x[i1][1];
diffz1 = x[j1][2] - x[i1][2];
// midpoint distance
dPx = 0.5*(diffx0-diffx1) + x[i0][0]-x[i1][0];
dPy = 0.5*(diffy0-diffy1) + x[i0][1]-x[i1][1];
dPz = 0.5*(diffz0-diffz1) + x[i0][2]-x[i1][2];
// Ri^2 Rj^2
RiRi = diffx0*diffx0 + diffy0*diffy0 + diffz0*diffz0;
RiRj = diffx0*diffx1 + diffy0*diffy1 + diffz0*diffz1;
RjRj = diffx1*diffx1 + diffy1*diffy1 + diffz1*diffz1;
denom = RiRj*RiRj - RiRi*RjRj;
// handle case of parallel lines
// reduce to midpt distance
if (fabs(denom) < SMALL){
if(denom < 0) denom = -BIG;
else denom = BIG;
}
// calc ti
termx0 = RiRj*diffx1 - RjRj*diffx0;
termy0 = RiRj*diffy1 - RjRj*diffy0;
termz0 = RiRj*diffz1 - RjRj*diffz0;
num0 = dPx*termx0 + dPy*termy0 + dPz*termz0;
ti = num0 / denom;
if (ti > 0.5) ti = 0.5;
if (ti < -0.5) ti = -0.5;
// calc tj
termx1 = RiRj*diffx0 - RiRi*diffx1;
termy1 = RiRj*diffy0 - RiRi*diffy1;
termz1 = RiRj*diffz0 - RiRi*diffz1;
num1 = dPx*termx1 + dPy*termy1 + dPz*termz1;
tj = -num1/ denom;
if (tj > 0.5) tj = 0.5;
if (tj < -0.5) tj = -0.5;
// min dist
dx = dPx - ti*diffx0 + tj*diffx1;
dy = dPy - ti*diffy0 + tj*diffy1;
dz = dPz - ti*diffz0 + tj*diffz1;
}
/* --------------------------------------------------------
map global id of atoms in stored by each bond particle
------------------------------------------------------- */
inline void PairSRP::remapBonds(int &nall)
{
if(nall > maxcount){
memory->grow(segment, nall, 2, "pair:segment");
maxcount = nall;
}
// loop over all bond particles
// each bond paricle holds two bond atoms
// map global ids of bond atoms to local ids
// might not be able to map both bond atoms of j, if j is outside neighcut
// these are not on neighlist, so are not used
int tmp;
srp = f_srp->array_atom;
for (int i = 0; i < nall; i++) {
if(atom->type[i] == bptype){
// tmp is local id
// tmp == -1 is ok
tmp = atom->map((int)srp[i][0]);
segment[i][0] = domain->closest_image(i,tmp);
// repeat with other id
tmp = atom->map((int)srp[i][1]);
segment[i][1] = domain->closest_image(i,tmp);
}
}
}
/* --------------------------------------------------------
add exclusions for 1-2 neighs, if requested
more complex exclusions or scaling probably not needed
------------------------------------------------------- */
inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &numneigh, int** &firstneigh)
{
int i0, i1, j0, j1;
int i,j,ii,jj,jnum;
// encode neighs with exclusions
// only need 1-2 info for normal uses of srp
// add 1-3, etc later if ever needed
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jnum = numneigh[i];
// two atoms inside bond particle
i0 = segment[i][0];
j0 = segment[i][1];
for (jj = 0; jj < jnum; jj++) {
jlist = firstneigh[i];
j = jlist[jj];
j &= NEIGHMASK;
//two atoms inside bond particle
i1 = segment[j][0];
j1 = segment[j][1];
// check for a 1-2 neigh
if(i0 == i1 || i0 == j1 || i1 == j0 || j0 == j1){
j |= ONETWOBIT;
jlist[jj] = j;
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairSRP::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d %g\n",i,a0[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairSRP::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g\n",i,j,a0[i][j],cut[i][j]);
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSRP::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&a0[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSRP::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
printf(" i %d j %d \n",i,j);
fread(&a0[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSRP::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&bptype,sizeof(int),1,fp);
fwrite(&btype,sizeof(int),1,fp);
fwrite(&min,sizeof(int),1,fp);
fwrite(&midpoint,sizeof(int),1,fp);
fwrite(&exclude,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSRP::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_global,sizeof(double),1,fp);
fread(&bptype,sizeof(int),1,fp);
fread(&btype,sizeof(int),1,fp);
fread(&min,sizeof(int),1,fp);
fread(&midpoint,sizeof(int),1,fp);
fread(&exclude,sizeof(int),1,fp);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
}
Event Timeline
Log In to Comment