ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2013-01-31 16:31:07 +00:00
parent 0c6b84987c
commit b648f55b2d
9 changed files with 1261 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
src/BODY/Install.sh Normal file
View File

@ -0,0 +1,27 @@
# Install/unInstall package files in LAMMPS
if (test $1 = 1) then
cp body_nparticle.cpp ..
cp compute_body_local.cpp ..
cp fix_nve_body.cpp ..
cp pair_body.cpp ..
cp body_nparticle.h ..
cp compute_body_local.h ..
cp fix_nve_body.h ..
cp pair_body.h ..
elif (test $1 = 0) then
rm -f ../body_nparticle.cpp
rm -f ../compute_body_local.cpp
rm -f ../fix_nve_body.cpp
rm -f ../pair_body.cpp
rm -f ../body_nparticle.h
rm -f ../compute_body_local.h
rm -f ../fix_nve_body.h
rm -f ../pair_body.h
fi

201
src/BODY/body_nparticle.cpp Normal file
View File

@ -0,0 +1,201 @@
/* ----------------------------------------------------------------------
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 "stdlib.h"
#include "body_nparticle.h"
#include "math_extra.h"
#include "atom_vec_body.h"
#include "atom.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EPSILON 1.0e-7
/* ---------------------------------------------------------------------- */
BodyNparticle::BodyNparticle(LAMMPS *lmp, int narg, char **arg) :
Body(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Invalid body nparticle command");
int nmin = atoi(arg[1]);
int nmax = atoi(arg[2]);
if (nmin <= 0 || nmin > nmax)
error->all(FLERR,"Invalid body nparticle command");
size_forward = 0;
size_border = 1 + 3*nmax;
}
/* ---------------------------------------------------------------------- */
int BodyNparticle::nsub(AtomVecBody::Bonus *bonus)
{
return bonus->ivalue[0];
}
/* ---------------------------------------------------------------------- */
double *BodyNparticle::coords(AtomVecBody::Bonus *bonus)
{
return bonus->dvalue;
}
/* ---------------------------------------------------------------------- */
int BodyNparticle::pack_border_body(AtomVecBody::Bonus *bonus, double *buf)
{
int nsub = bonus->ivalue[0];
buf[0] = nsub;
memcpy(&buf[1],bonus->dvalue,3*nsub*sizeof(double));
return 1+3*nsub;
}
/* ---------------------------------------------------------------------- */
int BodyNparticle::unpack_border_body(AtomVecBody::Bonus *bonus, double *buf)
{
int nsub = static_cast<int> (buf[0]);
bonus->ivalue[0] = nsub;
memcpy(bonus->dvalue,&buf[1],3*nsub*sizeof(double));
return 1+3*nsub;
}
/* ----------------------------------------------------------------------
populate bonus data structure with data file values
------------------------------------------------------------------------- */
void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
char **ifile, char **dfile)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
// error in data file if any values are NULL
for (int i = 0; i < ninteger; i++)
if (ifile[0] == NULL) error->one(FLERR,"");
for (int i = 0; i < ndouble; i++)
if (dfile[0] == NULL) error->one(FLERR,"");
// set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
if (ninteger != 1) error->one(FLERR,"");
int nsub = atoi(ifile[0]);
if (nsub < 1) error->one(FLERR,"");
if (ndouble != 6 + 3*nsub) error->one(FLERR,"");
bonus->ninteger = 1;
memory->create(bonus->ivalue,bonus->ninteger,"body:ivalue");
bonus->ivalue[0] = nsub;
bonus->ndouble = 3*nsub;
memory->create(bonus->dvalue,3*nsub,"body:dvalue");
// diagonalize inertia tensor
double tensor[3][3];
tensor[0][0] = atof(dfile[0]);
tensor[1][1] = atof(dfile[1]);
tensor[2][2] = atof(dfile[2]);
tensor[0][1] = tensor[1][0] = atof(dfile[3]);
tensor[0][2] = tensor[2][0] = atof(dfile[4]);
tensor[1][2] = tensor[2][1] = atof(dfile[5]);
double *inertia = bonus->inertia;
double evectors[3][3];
int ierror = MathExtra::jacobi(tensor,inertia,evectors);
if (ierror) error->one(FLERR,
"Insufficient Jacobi rotations for body nparticle");
// if any principal moment < scaled EPSILON, set to 0.0
double max;
max = MAX(inertia[0],inertia[1]);
max = MAX(max,inertia[2]);
if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
// exyz_space = principal axes in space frame
double ex_space[3],ey_space[3],ez_space[3];
ex_space[0] = evectors[0][0];
ex_space[1] = evectors[1][0];
ex_space[2] = evectors[2][0];
ey_space[0] = evectors[0][1];
ey_space[1] = evectors[1][1];
ey_space[2] = evectors[2][1];
ez_space[0] = evectors[0][2];
ez_space[1] = evectors[1][2];
ez_space[2] = evectors[2][2];
// enforce 3 evectors as a right-handed coordinate system
// flip 3rd vector if needed
double cross[3];
MathExtra::cross3(ex_space,ey_space,cross);
if (MathExtra::dot3(cross,ez_space) < 0.0) MathExtra::negate3(ez_space);
// create initial quaternion
MathExtra::exyz_to_q(ex_space,ey_space,ez_space,bonus->quat);
// bonus->dvalue = sub-particle displacements in body frame
double delta[3],displace[3];
int j = 6;
int k = 0;
for (int i = 0; i < nsub; i++) {
delta[0] = atof(dfile[j]);
delta[1] = atof(dfile[j+1]);
delta[2] = atof(dfile[j+2]);
MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
delta,&bonus->dvalue[k]);
j += 3;
k += 3;
}
}
/* ---------------------------------------------------------------------- */
int BodyNparticle::noutcol()
{
return 3;
}
/* ---------------------------------------------------------------------- */
int BodyNparticle::noutrow(int ibonus)
{
return avec->bonus[ibonus].ivalue[0];
}
/* ---------------------------------------------------------------------- */
void BodyNparticle::output(int ibonus, int m, double *values)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
double p[3][3];
MathExtra::quat_to_mat(bonus->quat,p);
MathExtra::matvec(p,&bonus->dvalue[3*m],values);
double *x = atom->x[bonus->ilocal];
values[0] += x[0];
values[1] += x[1];
values[2] += x[2];
}

57
src/BODY/body_nparticle.h Normal file
View File

@ -0,0 +1,57 @@
/* ----------------------------------------------------------------------
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 BODY_CLASS
BodyStyle(nparticle,BodyNparticle)
#else
#ifndef LMP_BODY_NPARTICLE_H
#define LMP_BODY_NPARTICLE_H
#include "body.h"
#include "atom_vec_body.h"
namespace LAMMPS_NS {
class BodyNparticle : public Body {
public:
BodyNparticle(class LAMMPS *, int, char **);
~BodyNparticle() {}
int nsub(class AtomVecBody::Bonus *);
double *coords(class AtomVecBody::Bonus *);
int pack_border_body(class AtomVecBody::Bonus *, double *);
int unpack_border_body(class AtomVecBody::Bonus *, double *);
void data_body(int, int, int, char **, char **);
int noutrow(int);
int noutcol();
void output(int, int, double *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
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.
*/

206
src/BODY/compute_body_local.cpp Normal file
View File

@ -0,0 +1,206 @@
/* ----------------------------------------------------------------------
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 "math.h"
#include "string.h"
#include "compute_body_local.h"
#include "atom.h"
#include "atom_vec_body.h"
#include "body.h"
#include "update.h"
#include "domain.h"
#include "force.h"
#include "bond.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define DELTA 10000
enum{TYPE,INDEX};
/* ---------------------------------------------------------------------- */
ComputeBodyLocal::ComputeBodyLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal compute body/local command");
local_flag = 1;
nvalues = narg - 3;
if (nvalues == 1) size_local_cols = 0;
else size_local_cols = nvalues;
which = new int[nvalues];
index = new int[nvalues];
nvalues = 0;
for (int iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"type") == 0) which[nvalues++] = TYPE;
else {
which[nvalues] = INDEX;
index[nvalues] = force->inumeric(arg[iarg]) - 1;
nvalues++;
}
}
avec = (AtomVecBody *) atom->style_match("body");
if (!avec) error->all(FLERR,"Compute body/local requires atom style body");
bptr = avec->bptr;
int indexmax = bptr->noutcol();
for (int i = 0; i < nvalues; i++) {
if (which[i] == INDEX && (index[i] < 0 || index[i] >= indexmax))
error->all(FLERR,"Invalid index in compute body/local command");
}
nmax = 0;
vector = NULL;
array = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeBodyLocal::~ComputeBodyLocal()
{
delete [] which;
delete [] index;
memory->destroy(vector);
memory->destroy(array);
}
/* ---------------------------------------------------------------------- */
void ComputeBodyLocal::init()
{
// do initial memory allocation so that memory_usage() is correct
int ncount = compute_body(0);
if (ncount > nmax) reallocate(ncount);
size_local_rows = ncount;
}
/* ---------------------------------------------------------------------- */
void ComputeBodyLocal::compute_local()
{
invoked_local = update->ntimestep;
// count local entries and compute body info
int ncount = compute_body(0);
if (ncount > nmax) reallocate(ncount);
size_local_rows = ncount;
ncount = compute_body(1);
}
/* ----------------------------------------------------------------------
count body info and compute body info on this proc
flag = 0, only count
flag = 1, also compute
------------------------------------------------------------------------- */
int ComputeBodyLocal::compute_body(int flag)
{
// perform count
int *mask = atom->mask;
int *body = atom->body;
int nlocal = atom->nlocal;
int ncount = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (body[i] < 0) ncount++;
else ncount += bptr->noutrow(body[i]);
}
if (flag == 0) return ncount;
// perform computation and fill output vector/array
int m,n,ibonus;
double *values = new double[bptr->noutcol()];
double **x = atom->x;
int *type = atom->type;
ncount = 0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (body[i] < 0) {
if (nvalues == 1) {
if (which[0] == TYPE) vector[ncount] = type[i];
else vector[ncount] = x[i][index[0]];
} else {
for (m = 0; m < nvalues; m++) {
if (which[m] == TYPE) array[ncount][m] = type[i];
else array[ncount][m] = x[i][index[m]];
}
}
ncount++;
} else {
ibonus = body[i];
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
n = bptr->noutrow(ibonus);
for (int j = 0; j < n; j++) {
bptr->output(ibonus,j,values);
if (nvalues == 1) {
if (which[0] == TYPE) vector[ncount] = type[i];
else vector[ncount] = values[index[0]];
} else {
for (m = 0; m < nvalues; m++) {
if (which[m] == TYPE) array[ncount][m] = type[i];
else array[ncount][m] = values[index[m]];
}
}
ncount++;
}
}
}
}
delete [] values;
return ncount;
}
/* ---------------------------------------------------------------------- */
void ComputeBodyLocal::reallocate(int n)
{
// grow vector or array and indices array
while (nmax < n) nmax += DELTA;
if (nvalues == 1) {
memory->destroy(vector);
memory->create(vector,nmax,"body/local:vector");
vector_local = vector;
} else {
memory->destroy(array);
memory->create(array,nmax,nvalues,"body/local:array");
array_local = array;
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeBodyLocal::memory_usage()
{
double bytes = nmax*nvalues * sizeof(double);
return bytes;
}

53
src/BODY/compute_body_local.h Normal file
View File

@ -0,0 +1,53 @@
/* ----------------------------------------------------------------------
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 COMPUTE_CLASS
ComputeStyle(body/local,ComputeBodyLocal)
#else
#ifndef LMP_COMPUTE_BODY_LOCAL_H
#define LMP_COMPUTE_BODY_LOCAL_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeBodyLocal : public Compute {
public:
ComputeBodyLocal(class LAMMPS *, int, char **);
~ComputeBodyLocal();
void init();
void compute_local();
double memory_usage();
private:
int nvalues;
int *which,*index;
int nmax;
double *vector;
double **array;
class AtomVecBody *avec;
class Body *bptr;
int compute_body(int);
void reallocate(int);
};
}
#endif
#endif

132
src/BODY/fix_nve_body.cpp Normal file
View File

@ -0,0 +1,132 @@
/* ----------------------------------------------------------------------
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 "math.h"
#include "stdio.h"
#include "string.h"
#include "fix_nve_body.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_body.h"
#include "force.h"
#include "update.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixNVEBody::FixNVEBody(LAMMPS *lmp, int narg, char **arg) :
FixNVE(lmp, narg, arg) {}
/* ---------------------------------------------------------------------- */
void FixNVEBody::init()
{
avec = (AtomVecBody *) atom->style_match("body");
if (!avec) error->all(FLERR,"Fix nve/body requires atom style body");
// check that all particles are bodies
// no point particles allowed
int *body = atom->body;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (body[i] < 0) error->one(FLERR,"Fix nve/body requires bodies");
FixNVE::init();
}
/* ---------------------------------------------------------------------- */
void FixNVEBody::initial_integrate(int vflag)
{
double dtfm;
double omega[3];
double *quat,*inertia;
AtomVecBody::Bonus *bonus = avec->bonus;
int *body = atom->body;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **angmom = atom->angmom;
double **torque = atom->torque;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// set timestep here since dt may have changed or come via rRESPA
dtq = 0.5 * dtv;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
// update angular momentum by 1/2 step
angmom[i][0] += dtf * torque[i][0];
angmom[i][1] += dtf * torque[i][1];
angmom[i][2] += dtf * torque[i][2];
// compute omega at 1/2 step from angmom at 1/2 step and current q
// update quaternion a full step via Richardson iteration
// returns new normalized quaternion
inertia = bonus[body[i]].inertia;
quat = bonus[body[i]].quat;
MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
MathExtra::richardson(quat,angmom[i],omega,inertia,dtq);
}
}
/* ---------------------------------------------------------------------- */
void FixNVEBody::final_integrate()
{
double dtfm;
double **v = atom->v;
double **f = atom->f;
double **angmom = atom->angmom;
double **torque = atom->torque;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
angmom[i][0] += dtf * torque[i][0];
angmom[i][1] += dtf * torque[i][1];
angmom[i][2] += dtf * torque[i][2];
}
}

41
src/BODY/fix_nve_body.h Normal file
View File

@ -0,0 +1,41 @@
/* ----------------------------------------------------------------------
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 FIX_CLASS
FixStyle(nve/body,FixNVEBody)
#else
#ifndef LMP_FIX_NVE_BODY_H
#define LMP_FIX_NVE_BODY_H
#include "fix_nve.h"
namespace LAMMPS_NS {
class FixNVEBody : public FixNVE {
public:
FixNVEBody(class LAMMPS *, int, char **);
void init();
void initial_integrate(int);
void final_integrate();
private:
double dtq;
class AtomVecBody *avec;
};
}
#endif
#endif

484
src/BODY/pair_body.cpp Normal file
View File

@ -0,0 +1,484 @@
/* ----------------------------------------------------------------------
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 "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_body.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_body.h"
#include "body_nparticle.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define DELTA 10000
/* ---------------------------------------------------------------------- */
PairBody::PairBody(LAMMPS *lmp) : Pair(lmp)
{
dmax = nmax = 0;
discrete = NULL;
dnum = dfirst = NULL;
single_enable = 0;
restartinfo = 0;
}
/* ---------------------------------------------------------------------- */
PairBody::~PairBody()
{
memory->destroy(discrete);
memory->destroy(dnum);
memory->destroy(dfirst);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
}
}
/* ---------------------------------------------------------------------- */
void PairBody::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
int ni,nj,npi,npj,ifirst,jfirst;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,r6inv,forcelj;
double xi[3],xj[3],fi[3],fj[3],ti[3],tj[3];
double *dxi,*dxj;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
double **torque = atom->torque;
int *body = atom->body;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// grow discrete list if necessary and initialize
if (nall > nmax) {
nmax = nall;
memory->destroy(dnum);
memory->destroy(dfirst);
memory->create(dnum,nall,"pair:dnum");
memory->create(dfirst,nall,"pair:dfirst");
}
for (i = 0; i < nall; i++) dnum[i] = 0;
ndiscrete = 0;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq >= cutsq[itype][jtype]) continue;
// body/body interactions = NxM sub-particles
evdwl = 0.0;
if (body[i] >= 0 && body[j] >= 0) {
if (dnum[i] == 0) body2space(i);
npi = dnum[i];
ifirst = dfirst[i];
if (dnum[j] == 0) body2space(j);
npj = dnum[j];
jfirst = dfirst[j];
for (ni = 0; ni < npi; ni++) {
dxi = discrete[ifirst+ni];
for (nj = 0; nj < npj; nj++) {
dxj = discrete[jfirst+nj];
xi[0] = x[i][0] + dxi[0];
xi[1] = x[i][1] + dxi[1];
xi[2] = x[i][2] + dxi[2];
xj[0] = x[j][0] + dxj[0];
xj[1] = x[j][1] + dxj[1];
xj[2] = x[j][2] + dxj[2];
delx = xi[0] - xj[0];
dely = xi[1] - xj[1];
delz = xi[2] - xj[2];
rsq = delx*delx + dely*dely + delz*delz;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
fpair = forcelj*r2inv;
if (eflag)
evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
fi[0] = delx*fpair;
fi[1] = dely*fpair;
fi[2] = delz*fpair;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
ti[0] = dxi[1]*fi[2] - dxi[2]*fi[1];
ti[1] = dxi[2]*fi[0] - dxi[0]*fi[2];
ti[2] = dxi[0]*fi[1] - dxi[1]*fi[0];
torque[i][0] += ti[0];
torque[i][1] += ti[1];
torque[i][2] += ti[2];
if (newton_pair || j < nlocal) {
fj[0] = -delx*fpair;
fj[1] = -dely*fpair;
fj[2] = -delz*fpair;
f[j][0] += fj[0];
f[j][1] += fj[1];
f[j][2] += fj[2];
tj[0] = dxj[1]*fj[2] - dxj[2]*fj[1];
tj[1] = dxj[2]*fj[0] - dxj[0]*fj[2];
tj[2] = dxj[0]*fj[1] - dxj[1]*fj[0];
torque[j][0] += tj[0];
torque[j][1] += tj[1];
torque[j][2] += tj[2];
}
}
}
// body/particle interaction = Nx1 sub-particles
} else if (body[i] >= 0) {
if (dnum[i] == 0) body2space(i);
npi = dnum[i];
ifirst = dfirst[i];
for (ni = 0; ni < npi; ni++) {
dxi = discrete[ifirst+ni];
xi[0] = x[i][0] + dxi[0];
xi[1] = x[i][1] + dxi[1];
xi[2] = x[i][2] + dxi[2];
xj[0] = x[j][0];
xj[1] = x[j][1];
xj[2] = x[j][2];
delx = xi[0] - xj[0];
dely = xi[1] - xj[1];
delz = xi[2] - xj[2];
rsq = delx*delx + dely*dely + delz*delz;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
fpair = forcelj*r2inv;
if (eflag)
evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
fi[0] = delx*fpair;
fi[1] = dely*fpair;
fi[2] = delz*fpair;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
ti[0] = dxi[1]*fi[2] - dxi[2]*fi[1];
ti[1] = dxi[2]*fi[0] - dxi[0]*fi[2];
ti[2] = dxi[0]*fi[1] - dxi[1]*fi[0];
torque[i][0] += ti[0];
torque[i][1] += ti[1];
torque[i][2] += ti[2];
if (newton_pair || j < nlocal) {
fj[0] = -delx*fpair;
fj[1] = -dely*fpair;
fj[2] = -delz*fpair;
f[j][0] += fj[0];
f[j][1] += fj[1];
f[j][2] += fj[2];
}
}
// particle/body interaction = Nx1 sub-particles
} else if (body[j] >= 0) {
if (dnum[j] == 0) body2space(j);
npj = dnum[j];
jfirst = dfirst[j];
for (nj = 0; nj < npj; nj++) {
dxj = discrete[jfirst+nj];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
xj[0] = x[j][0] + dxj[0];
xj[1] = x[j][1] + dxj[1];
xj[2] = x[j][2] + dxj[2];
delx = xi[0] - xj[0];
dely = xi[1] - xj[1];
delz = xi[2] - xj[2];
rsq = delx*delx + dely*dely + delz*delz;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
fpair = forcelj*r2inv;
if (eflag)
evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
fi[0] = delx*fpair;
fi[1] = dely*fpair;
fi[2] = delz*fpair;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
fj[0] = -delx*fpair;
fj[1] = -dely*fpair;
fj[2] = -delz*fpair;
f[j][0] += fj[0];
f[j][1] += fj[1];
f[j][2] += fj[2];
tj[0] = dxj[1]*fj[2] - dxj[2]*fj[1];
tj[1] = dxj[2]*fj[0] - dxj[0]*fj[2];
tj[2] = dxj[0]*fj[1] - dxj[1]*fj[0];
torque[j][0] += tj[0];
torque[j][1] += tj[1];
torque[j][2] += tj[2];
}
}
// particle/particle interaction = 1x1 particles
} else {
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
fpair = forcelj*r2inv;
if (eflag)
evdwl += r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairBody::allocate()
{
allocated = 1;
int n = atom->ntypes;
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;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairBody::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
cut_global = force->numeric(arg[0]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairBody::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double epsilon_one = force->numeric(arg[2]);
double sigma_one = force->numeric(arg[3]);
double cut_one = cut_global;
if (narg == 5) cut_one = force->numeric(arg[4]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairBody::init_style()
{
avec = (AtomVecBody *) atom->style_match("body");
if (!avec) error->all(FLERR,"Pair body requires atom style body");
if (strcmp(avec->bptr->style,"nparticle") != 0)
error->all(FLERR,"Pair body requires body style nparticle");
bptr = (BodyNparticle *) avec->bptr;
neighbor->request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairBody::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
}
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
epsilon[j][i] = epsilon[i][j];
sigma[j][i] = sigma[i][j];
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
convert N sub-particles in body I to space frame using current quaternion
store sub-particle space-frame displacements from COM in discrete list
------------------------------------------------------------------------- */
void PairBody::body2space(int i)
{
int ibonus = atom->body[i];
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
int nsub = bptr->nsub(bonus);
double *coords = bptr->coords(bonus);
dnum[i] = nsub;
dfirst[i] = ndiscrete;
if (ndiscrete + nsub > dmax) {
dmax += DELTA;
memory->grow(discrete,dmax,3,"pair:discrete");
}
double p[3][3];
MathExtra::quat_to_mat(bonus->quat,p);
for (int m = 0; m < nsub; m++) {
MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
ndiscrete++;
}
}

60
src/BODY/pair_body.h Normal file
View File

@ -0,0 +1,60 @@
/* ----------------------------------------------------------------------
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 PAIR_CLASS
PairStyle(body,PairBody)
#else
#ifndef LMP_PAIR_BODY_H
#define LMP_PAIR_BODY_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBody : public Pair {
public:
PairBody(class LAMMPS *);
~PairBody();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
protected:
double cut_global;
double **cut;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4;
class AtomVecBody *avec;
class BodyNparticle *bptr;
double **discrete; // list of all sub-particles for all bodies
int ndiscrete; // number of discretes in list
int dmax; // allocated size of discrete list
int *dnum; // number of discretes per line, 0 if uninit
int *dfirst; // index of first discrete per each line
int nmax; // allocated size of dnum,dfirst vectors
void allocate();
void body2space(int);
};
}
#endif
#endif