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@7252 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2011-12-01 14:45:19 +00:00
parent b7f57f079d
commit badb8c439d
7 changed files with 2005 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

494
src/USER-CG-CMM/angle_sdk.cpp Normal file
View File

@ -0,0 +1,494 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Variant of the harmonic angle potential for use with the
lj/sdk potential for coarse grained MD simulations.
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "angle_sdk.h"
#include "atom.h"
#include "neighbor.h"
#include "pair.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "lj_sdk_common.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace LJSDKParms;
#define SMALL 0.001
/* ---------------------------------------------------------------------- */
AngleSDK::AngleSDK(LAMMPS *lmp) : Angle(lmp) { repflag = 0;}
/* ---------------------------------------------------------------------- */
AngleSDK::~AngleSDK()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(k);
memory->destroy(theta0);
memory->destroy(repscale);
allocated = 0;
}
}
/* ---------------------------------------------------------------------- */
void AngleSDK::compute(int eflag, int vflag)
{
int i1,i2,i3,n,type;
double delx1,dely1,delz1,delx2,dely2,delz2,delx3,dely3,delz3;
double eangle,f1[3],f3[3],e13,f13;
double dtheta,tk;
double rsq1,rsq2,rsq3,r1,r2,c,s,a,a11,a12,a22;
eangle = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = 0;
double **x = atom->x;
double **f = atom->f;
int **anglelist = neighbor->anglelist;
int nanglelist = neighbor->nanglelist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
for (n = 0; n < nanglelist; n++) {
i1 = anglelist[n][0];
i2 = anglelist[n][1];
i3 = anglelist[n][2];
type = anglelist[n][3];
// 1st bond
delx1 = x[i1][0] - x[i2][0];
dely1 = x[i1][1] - x[i2][1];
delz1 = x[i1][2] - x[i2][2];
domain->minimum_image(delx1,dely1,delz1);
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[i3][0] - x[i2][0];
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
domain->minimum_image(delx2,dely2,delz2);
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL;
s = 1.0/s;
// 1-3 LJ interaction.
// we only want to use the repulsive part,
// and it can be scaled (or off).
// so this has to be done here and not in the
// general non-bonded code.
if (repflag) {
delx3 = x[i1][0] - x[i3][0];
dely3 = x[i1][1] - x[i3][1];
delz3 = x[i1][2] - x[i3][2];
domain->minimum_image(delx3,dely3,delz3);
rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3;
const int type1 = atom->type[i1];
const int type3 = atom->type[i3];
f13=0.0;
e13=0.0;
if (rsq3 < rminsq[type1][type3]) {
const int ljt = lj_type[type1][type3];
const double r2inv = 1.0/rsq3;
if (ljt == LJ12_4) {
const double r4inv=r2inv*r2inv;
f13 = r4inv*(lj1[type1][type3]*r4inv*r4inv - lj2[type1][type3]);
if (eflag) e13 = r4inv*(lj3[type1][type3]*r4inv*r4inv - lj4[type1][type3]);
} else if (ljt == LJ9_6) {
const double r3inv = r2inv*sqrt(r2inv);
const double r6inv = r3inv*r3inv;
f13 = r6inv*(lj1[type1][type3]*r3inv - lj2[type1][type3]);
if (eflag) e13 = r6inv*(lj3[type1][type3]*r3inv - lj4[type1][type3]);
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
f13 = r6inv*(lj1[type1][type3]*r6inv - lj2[type1][type3]);
if (eflag) e13 = r6inv*(lj3[type1][type3]*r6inv - lj4[type1][type3]);
}
// make sure energy is 0.0 at the cutoff.
if (eflag) e13 -= emin[type1][type3];
f13 *= r2inv;
}
}
// force & energy
dtheta = acos(c) - theta0[type];
tk = k[type] * dtheta;
if (eflag) eangle = tk*dtheta;
a = -2.0 * tk * s;
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
a22 = a*c / rsq2;
f1[0] = a11*delx1 + a12*delx2;
f1[1] = a11*dely1 + a12*dely2;
f1[2] = a11*delz1 + a12*delz2;
f3[0] = a22*delx2 + a12*delx1;
f3[1] = a22*dely2 + a12*dely1;
f3[2] = a22*delz2 + a12*delz1;
// apply force to each of the 3 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0] + f13*delx3;
f[i1][1] += f1[1] + f13*dely3;
f[i1][2] += f1[2] + f13*delz3;
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= f1[0] + f3[0];
f[i2][1] -= f1[1] + f3[1];
f[i2][2] -= f1[2] + f3[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0] - f13*delx3;
f[i3][1] += f3[1] - f13*dely3;
f[i3][2] += f3[2] - f13*delz3;
}
if (evflag) {
ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2);
ev_tally13(i1,i3,nlocal,newton_bond,e13,f13,delx3,dely3,delz3);
}
}
}
/* ---------------------------------------------------------------------- */
void AngleSDK::allocate()
{
allocated = 1;
int n = atom->nangletypes;
memory->create(k,n+1,"angle:k");
memory->create(theta0,n+1,"angle:theta0");
memory->create(repscale,n+1,"angle:repscale");
memory->create(setflag,n+1,"angle:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
void AngleSDK::coeff(int narg, char **arg)
{
if ((narg < 3) || (narg > 6))
error->all(FLERR,"Incorrect args for angle coefficients");
if (!allocated) allocate();
int ilo,ihi;
force->bounds(arg[0],atom->nangletypes,ilo,ihi);
double k_one = force->numeric(arg[1]);
double theta0_one = force->numeric(arg[2]);
double repscale_one;
// backward compatibility with old cg/cmm style input:
// this had <lj_type> <epsilon> <sigma>
// if epsilon is set to 0.0 we accept it as repscale 0.0
// otherwise assume repscale 1.0, since we were using
// epsilon to turn repulsion on or off.
if (narg == 6) {
repscale_one = force->numeric(arg[4]);
if (repscale_one > 0.0) repscale_one = 1.0;
} else if (narg == 4) repscale_one = force->numeric(arg[3]);
else if (narg == 3) repscale_one = 1.0;
else error->all(FLERR,"Incorrect args for angle coefficients");
// convert theta0 from degrees to radians and store coefficients
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
theta0[i] = theta0_one/180.0 * MY_PI;
repscale[i] = repscale_one;
setflag[i] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
}
/* ----------------------------------------------------------------------
error check and initialize all values needed for force computation
------------------------------------------------------------------------- */
void AngleSDK::init_style()
{
// make sure we use an SDK pair_style and that we need the 1-3 repulsion
repflag = 0;
for (int i = 1; i <= atom->nangletypes; i++)
if (repscale[i] > 0.0) repflag = 1;
// set up pointers to access SDK LJ parameters for 1-3 interactions
if (repflag) {
int itmp;
if (force->pair == NULL)
error->all(FLERR,"Angle style SDK requires use of a compatible with Pair style");
lj1 = (double **) force->pair->extract("lj1",itmp);
lj2 = (double **) force->pair->extract("lj2",itmp);
lj3 = (double **) force->pair->extract("lj3",itmp);
lj4 = (double **) force->pair->extract("lj4",itmp);
lj_type = (int **) force->pair->extract("lj_type",itmp);
rminsq = (double **) force->pair->extract("rminsq",itmp);
emin = (double **) force->pair->extract("emin",itmp);
if (!lj1 || !lj2 || !lj3 || !lj4 || !lj_type || !rminsq || !emin)
error->all(FLERR,"Angle style SDK is incompatible with Pair style");
}
}
/* ---------------------------------------------------------------------- */
double AngleSDK::equilibrium_angle(int i)
{
return theta0[i];
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void AngleSDK::write_restart(FILE *fp)
{
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
fwrite(&repscale[1],sizeof(double),atom->nangletypes,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void AngleSDK::read_restart(FILE *fp)
{
allocate();
if (comm->me == 0) {
fread(&k[1],sizeof(double),atom->nangletypes,fp);
fread(&theta0[1],sizeof(double),atom->nangletypes,fp);
fread(&repscale[1],sizeof(double),atom->nangletypes,fp);
}
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&repscale[1],atom->nangletypes,MPI_DOUBLE,0,world);
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
}
/* ---------------------------------------------------------------------- */
void AngleSDK::ev_tally13(int i, int j, int nlocal, int newton_bond,
double evdwl, double fpair,
double delx, double dely, double delz)
{
double v[6];
if (eflag_either) {
if (eflag_global) {
if (newton_bond) {
energy += evdwl;
} else {
if (i < nlocal)
energy += 0.5*evdwl;
if (j < nlocal)
energy += 0.5*evdwl;
}
}
if (eflag_atom) {
if (newton_bond || i < nlocal) eatom[i] += 0.5*evdwl;
if (newton_bond || j < nlocal) eatom[i] += 0.5*evdwl;
}
}
if (vflag_either) {
v[0] = delx*delx*fpair;
v[1] = dely*dely*fpair;
v[2] = delz*delz*fpair;
v[3] = delx*dely*fpair;
v[4] = delx*delz*fpair;
v[5] = dely*delz*fpair;
if (vflag_global) {
if (newton_bond) {
virial[0] += v[0];
virial[1] += v[1];
virial[2] += v[2];
virial[3] += v[3];
virial[4] += v[4];
virial[5] += v[5];
} else {
if (i < nlocal) {
virial[0] += 0.5*v[0];
virial[1] += 0.5*v[1];
virial[2] += 0.5*v[2];
virial[3] += 0.5*v[3];
virial[4] += 0.5*v[4];
virial[5] += 0.5*v[5];
}
if (j < nlocal) {
virial[0] += 0.5*v[0];
virial[1] += 0.5*v[1];
virial[2] += 0.5*v[2];
virial[3] += 0.5*v[3];
virial[4] += 0.5*v[4];
virial[5] += 0.5*v[5];
}
}
}
if (vflag_atom) {
if (newton_bond || i < nlocal) {
vatom[i][0] += 0.5*v[0];
vatom[i][1] += 0.5*v[1];
vatom[i][2] += 0.5*v[2];
vatom[i][3] += 0.5*v[3];
vatom[i][4] += 0.5*v[4];
vatom[i][5] += 0.5*v[5];
}
if (newton_bond || j < nlocal) {
vatom[j][0] += 0.5*v[0];
vatom[j][1] += 0.5*v[1];
vatom[j][2] += 0.5*v[2];
vatom[j][3] += 0.5*v[3];
vatom[j][4] += 0.5*v[4];
vatom[j][5] += 0.5*v[5];
}
}
}
}
/* ---------------------------------------------------------------------- */
double AngleSDK::single(int type, int i1, int i2, int i3)
{
double **x = atom->x;
double delx1 = x[i1][0] - x[i2][0];
double dely1 = x[i1][1] - x[i2][1];
double delz1 = x[i1][2] - x[i2][2];
domain->minimum_image(delx1,dely1,delz1);
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
double delx2 = x[i3][0] - x[i2][0];
double dely2 = x[i3][1] - x[i2][1];
double delz2 = x[i3][2] - x[i2][2];
domain->minimum_image(delx2,dely2,delz2);
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
double e13=0.0;
if (repflag) {
// 1-3 LJ interaction.
double delx3 = x[i1][0] - x[i3][0];
double dely3 = x[i1][1] - x[i3][1];
double delz3 = x[i1][2] - x[i3][2];
domain->minimum_image(delx3,dely3,delz3);
const int type1 = atom->type[i1];
const int type3 = atom->type[i3];
const double rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3;
if (rsq3 < rminsq[type1][type3]) {
const int ljt = lj_type[type1][type3];
const double r2inv = 1.0/rsq3;
if (ljt == LJ12_4) {
const double r4inv=r2inv*r2inv;
e13 = r4inv*(lj3[type1][type3]*r4inv*r4inv - lj4[type1][type3]);
} else if (ljt == LJ9_6) {
const double r3inv = r2inv*sqrt(r2inv);
const double r6inv = r3inv*r3inv;
e13 = r6inv*(lj3[type1][type3]*r3inv - lj4[type1][type3]);
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
e13 = r6inv*(lj3[type1][type3]*r6inv - lj4[type1][type3]);
}
// make sure energy is 0.0 at the cutoff.
e13 -= emin[type1][type3];
}
}
double dtheta = acos(c) - theta0[type];
double tk = k[type] * dtheta;
return tk*dtheta + e13;
}

62
src/USER-CG-CMM/angle_sdk.h Normal file
View File

@ -0,0 +1,62 @@
/* -*- c++ -*- ----------------------------------------------------------
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 ANGLE_CLASS
AngleStyle(sdk,AngleSDK)
AngleStyle(cg/cmm,AngleSDK)
#else
#ifndef LMP_ANGLE_SDK_H
#define LMP_ANGLE_SDK_H
#include "stdio.h"
#include "angle.h"
namespace LAMMPS_NS {
class AngleSDK : public Angle {
public:
AngleSDK(class LAMMPS *);
virtual ~AngleSDK();
virtual void compute(int, int);
void coeff(int, char **);
void init_style();
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
protected:
double *k,*theta0;
// scaling factor for repulsive 1-3 interaction
double *repscale;
// parameters from SDK pair style
int **lj_type;
double **lj1,**lj2, **lj3, **lj4;
double **rminsq,**emin;
int repflag; // 1 if we have to handle 1-3 repulsion
void ev_tally13(int, int, int, int, double, double,
double, double, double);
void allocate();
};
}
#endif
#endif

44
src/USER-CG-CMM/lj_sdk_common.h Normal file
View File

@ -0,0 +1,44 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Common data for the Shinoda, DeVane, Klein (SDK) coarse grain model
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifndef LMP_LJ_SDK_COMMON_H
#define LMP_LJ_SDK_COMMON_H
#include "string.h"
namespace LAMMPS_NS {
namespace LJSDKParms {
// LJ type flags. list of supported LJ exponent combinations
enum {LJ_NOT_SET=0, LJ9_6, LJ12_4, LJ12_6, NUM_LJ_TYPES};
static int find_lj_type(const char *label,
const char * const * const list) {
for (int i=0; i < NUM_LJ_TYPES; ++i)
if (strcmp(label,list[i]) == 0) return i;
return LJ_NOT_SET;
}
static const char * const lj_type_list[] = {"none", "lj9_6", "lj12_4", "lj12_6"};
static const double lj_prefact[] = {0.0, 6.75, 2.59807621135332, 4.0};
static const double lj_pow1[] = {0.0, 9.00, 12.0, 12.0};
static const double lj_pow2[] = {0.0, 6.00, 4.0, 6.0};
}}
#endif

484
src/USER-CG-CMM/pair_lj_sdk.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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
This style is a simplified re-implementation of the CG/CMM pair style
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lj_sdk.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "lj_sdk_common.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace LJSDKParms;
/* ---------------------------------------------------------------------- */
PairLJSDK::PairLJSDK(LAMMPS *lmp) : Pair(lmp)
{
respa_enable = 0;
single_enable = 1;
}
/* ---------------------------------------------------------------------- */
PairLJSDK::~PairLJSDK()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(lj_type);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(offset);
memory->destroy(rminsq);
memory->destroy(emin);
allocated = 0;
}
}
/* ---------------------------------------------------------------------- */
void PairLJSDK::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 0;
if (evflag) {
if (eflag) {
if (force->newton_pair) eval<1,1,1>();
else eval<1,1,0>();
} else {
if (force->newton_pair) eval<1,0,1>();
else eval<1,0,0>();
}
} else {
if (force->newton_pair) eval<0,0,1>();
else eval<0,0,0>();
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairLJSDK::eval()
{
int i,j,ii,jj,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,forcelj,factor_lj;
evdwl = 0.0;
const double * const * const x = atom->x;
double * const * const f = atom->f;
const int * const type = atom->type;
const int nlocal = atom->nlocal;
const double * const special_lj = force->special_lj;
double fxtmp,fytmp,fztmp;
const int inum = list->inum;
const int * const ilist = list->ilist;
const int * const numneigh = list->numneigh;
const int * const * const firstneigh = list->firstneigh;
// 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];
fxtmp=fytmp=fztmp=0.0;
const int itype = type[i];
const int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
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]) {
r2inv = 1.0/rsq;
const int ljt = lj_type[itype][jtype];
if (ljt == LJ12_4) {
const double r4inv=r2inv*r2inv;
forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv
- lj2[itype][jtype]);
if (EFLAG)
evdwl = r4inv*(lj3[itype][jtype]*r4inv*r4inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ9_6) {
const double r3inv = r2inv*sqrt(r2inv);
const double r6inv = r3inv*r3inv;
forcelj = r6inv*(lj1[itype][jtype]*r3inv
- lj2[itype][jtype]);
if (EFLAG)
evdwl = r6inv*(lj3[itype][jtype]*r3inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv*(lj1[itype][jtype]*r6inv
- lj2[itype][jtype]);
if (EFLAG)
evdwl = r6inv*(lj3[itype][jtype]*r6inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else continue;
fpair = factor_lj*forcelj*r2inv;
fxtmp += delx*fpair;
fytmp += dely*fpair;
fztmp += delz*fpair;
if (NEWTON_PAIR || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (EFLAG) evdwl *= factor_lj;
if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR,
evdwl,0.0,fpair,delx,dely,delz);
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJSDK::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(lj_type,n+1,n+1,"pair:lj_type");
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
setflag[i][j] = 0;
lj_type[i][j] = LJ_NOT_SET;
}
}
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
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");
memory->create(offset,n+1,n+1,"pair:offset");
memory->create(rminsq,n+1,n+1,"pair:rminsq");
memory->create(emin,n+1,n+1,"pair:emin");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJSDK::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 PairLJSDK::coeff(int narg, char **arg)
{
if (narg < 5 || narg > 6) 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);
int lj_type_one = find_lj_type(arg[2],lj_type_list);
if (lj_type_one == LJ_NOT_SET)
error->all(FLERR,"Cannot parse LJ type flag.");
double epsilon_one = force->numeric(arg[3]);
double sigma_one = force->numeric(arg[4]);
double cut_one = cut_global;
if (narg == 6) cut_one = force->numeric(arg[5]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
lj_type[i][j] = lj_type_one;
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 for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJSDK::init_one(int i, int j)
{
if (setflag[i][j] == 0)
error->all(FLERR,"No mixing support for lj/sdk. "
"Coefficients for all pairs need to be set explicitly.");
const int ljt = lj_type[i][j];
if (ljt == LJ_NOT_SET)
error->all(FLERR,"unrecognized LJ parameter flag");
lj1[i][j] = lj_prefact[ljt] * lj_pow1[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow1[ljt]);
lj2[i][j] = lj_prefact[ljt] * lj_pow2[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow2[ljt]);
lj3[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow1[ljt]);
lj4[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow2[ljt]);
if (offset_flag) {
double ratio = sigma[i][j] / cut[i][j];
offset[i][j] = lj_prefact[ljt] * epsilon[i][j] * (pow(ratio,lj_pow1[ljt]) - pow(ratio,lj_pow2[ljt]));
} else offset[i][j] = 0.0;
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
cut[j][i] = cut[i][j];
cutsq[j][i] = cutsq[i][j];
offset[j][i] = offset[i][j];
lj_type[j][i] = lj_type[i][j];
// compute derived parameters for SDK angle potential
const double eps = epsilon[i][j];
const double sig = sigma[i][j];
const double rmin = sig*exp(1.0/(lj_pow1[ljt]-lj_pow2[ljt])
*log(lj_pow1[ljt]/lj_pow2[ljt]) );
rminsq[j][i] = rminsq[i][j] = rmin*rmin;
const double ratio = sig/rmin;
const double emin_one = lj_prefact[ljt] * eps * (pow(ratio,lj_pow1[ljt])
- pow(ratio,lj_pow2[ljt]));
emin[j][i] = emin[i][j] = emin_one;
// compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce
if (tail_flag)
error->all(FLERR,"Tail flag not supported by lj/sdk pair style");
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDK::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(&lj_type[i][j],sizeof(int),1,fp);
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDK::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) {
fread(&lj_type[i][j],sizeof(int),1,fp);
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&lj_type[i][j],1,MPI_INT,0,world);
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDK::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDK::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&cut_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairLJSDK::single(int, int, int itype, int jtype, double rsq,
double, double factor_lj, double &fforce)
{
if (rsq < cutsq[itype][jtype]) {
const int ljt = lj_type[itype][jtype];
const double ljpow1 = lj_pow1[ljt];
const double ljpow2 = lj_pow2[ljt];
const double ljpref = lj_prefact[ljt];
const double ratio = sigma[itype][jtype]/sqrt(rsq);
const double eps = epsilon[itype][jtype];
fforce = factor_lj * ljpref*eps * (ljpow1*pow(ratio,ljpow1)
- ljpow2*pow(ratio,ljpow2))/rsq;
return factor_lj * (ljpref*eps * (pow(ratio,ljpow1) - pow(ratio,ljpow2))
- offset[itype][jtype]);
} else fforce=0.0;
return 0.0;
}
/* ---------------------------------------------------------------------- */
void *PairLJSDK::extract(char *str, int &dim)
{
dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
if (strcmp(str,"sigma") == 0) return (void *) sigma;
if (strcmp(str,"lj_type") == 0) return (void *) lj_type;
if (strcmp(str,"lj1") == 0) return (void *) lj1;
if (strcmp(str,"lj2") == 0) return (void *) lj2;
if (strcmp(str,"lj3") == 0) return (void *) lj3;
if (strcmp(str,"lj4") == 0) return (void *) lj4;
if (strcmp(str,"rminsq") == 0) return (void *) rminsq;
if (strcmp(str,"emin") == 0) return (void *) emin;
return NULL;
}
/* ---------------------------------------------------------------------- */
double PairLJSDK::memory_usage()
{
double bytes = Pair::memory_usage();
int n = atom->ntypes;
// setflag/lj_type
bytes += 2 * (n+1)*(n+1)*sizeof(int);
// cut/cutsq/epsilon/sigma/offset/lj1/lj2/lj3/lj4/rminsq/emin
bytes += 11 * (n+1)*(n+1)*sizeof(double);
return bytes;
}

74
src/USER-CG-CMM/pair_lj_sdk.h Normal file
View File

@ -0,0 +1,74 @@
/* -*- c++ -*- ----------------------------------------------------------
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 author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lj/sdk,PairLJSDK)
PairStyle(cg/cmm,PairLJSDK)
#else
#ifndef LMP_PAIR_LJ_SDK_H
#define LMP_PAIR_LJ_SDK_H
#include "pair.h"
namespace LAMMPS_NS {
class LAMMPS;
class PairLJSDK : public Pair {
public:
PairLJSDK(LAMMPS *);
virtual ~PairLJSDK();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(char *, int &);
virtual double memory_usage();
protected:
int **lj_type; // type of lennard jones potential
double **cut;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset;
// cutoff and offset for minimum of LJ potential
// to be used in SDK angle potential, which
// uses only the repulsive part of the potential
double **rminsq, **emin;
double cut_global;
void allocate();
private:
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval();
};
}
#endif
#endif

765
src/USER-CG-CMM/pair_lj_sdk_coul_long.cpp Normal file
View File

@ -0,0 +1,765 @@
/* ----------------------------------------------------------------------
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 author: Axel Kohlmeyer (Temple U)
This style is a simplified re-implementation of the CG/CMM pair style
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lj_sdk_coul_long.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "kspace.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "lj_sdk_common.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace LJSDKParms;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
PairLJSDKCoulLong::PairLJSDKCoulLong(LAMMPS *lmp) : Pair(lmp)
{
respa_enable = 0;
ftable = NULL;
}
/* ---------------------------------------------------------------------- */
PairLJSDKCoulLong::~PairLJSDKCoulLong()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(lj_type);
memory->destroy(cutsq);
memory->destroy(cut_lj);
memory->destroy(cut_ljsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(offset);
memory->destroy(rminsq);
memory->destroy(emin);
allocated = 0;
}
if (ftable) free_tables();
}
/* ---------------------------------------------------------------------- */
void PairLJSDKCoulLong::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 0;
if (evflag) {
if (eflag) {
if (force->newton_pair) eval<1,1,1>();
else eval<1,1,0>();
} else {
if (force->newton_pair) eval<1,0,1>();
else eval<1,0,0>();
}
} else {
if (force->newton_pair) eval<0,0,1>();
else eval<0,0,0>();
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairLJSDKCoulLong::eval()
{
int i,ii,j,jj,jtype,itable;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double fraction,table;
double r,rsq,r2inv,forcecoul,forcelj,factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
evdwl = ecoul = 0.0;
const double * const * const x = atom->x;
double * const * const f = atom->f;
const double * const q = atom->q;
const int * const type = atom->type;
const int nlocal = atom->nlocal;
const double * const special_coul = force->special_coul;
const double * const special_lj = force->special_lj;
const double qqrd2e = force->qqrd2e;
double fxtmp,fytmp,fztmp;
const int inum = list->inum;
const int * const ilist = list->ilist;
const int * const numneigh = list->numneigh;
const int * const * const firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
fxtmp=fytmp=fztmp=0.0;
const int itype = type[i];
const int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
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]) {
r2inv = 1.0/rsq;
const int ljt = lj_type[itype][jtype];
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
}
} else {
forcecoul = 0.0;
ecoul = 0.0;
}
if (rsq < cut_ljsq[itype][jtype]) {
if (ljt == LJ12_4) {
const double r4inv=r2inv*r2inv;
forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv
- lj2[itype][jtype]);
if (EFLAG)
evdwl = r4inv*(lj3[itype][jtype]*r4inv*r4inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ9_6) {
const double r3inv = r2inv*sqrt(r2inv);
const double r6inv = r3inv*r3inv;
forcelj = r6inv*(lj1[itype][jtype]*r3inv
- lj2[itype][jtype]);
if (EFLAG)
evdwl = r6inv*(lj3[itype][jtype]*r3inv
- lj4[itype][jtype]) - offset[itype][jtype];
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv*(lj1[itype][jtype]*r6inv
- lj2[itype][jtype]);
if (EFLAG)
evdwl = r6inv*(lj3[itype][jtype]*r6inv
- lj4[itype][jtype]) - offset[itype][jtype];
}
} else {
forcelj=0.0;
evdwl = 0.0;
}
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
fxtmp += delx*fpair;
fytmp += dely*fpair;
fztmp += delz*fpair;
if (NEWTON_PAIR || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (EFLAG) {
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq)
ecoul = prefactor*erfc;
else {
table = etable[itable] + fraction*detable[itable];
ecoul = qtmp*q[j] * table;
}
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
} else ecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
evdwl *= factor_lj;
} else evdwl = 0.0;
}
if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR,
evdwl,ecoul,fpair,delx,dely,delz);
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(lj_type,n+1,n+1,"pair:lj_type");
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
setflag[i][j] = 0;
lj_type[i][j] = LJ_NOT_SET;
}
}
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut_lj,n+1,n+1,"pair:cut_lj");
memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
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");
memory->create(offset,n+1,n+1,"pair:offset");
memory->create(rminsq,n+1,n+1,"pair:rminsq");
memory->create(emin,n+1,n+1,"pair:emin");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
cut_lj_global = force->numeric(arg[0]);
if (narg == 1) cut_coul = cut_lj_global;
else cut_coul = force->numeric(arg[1]);
// 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_lj[i][j] = cut_lj_global;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::coeff(int narg, char **arg)
{
if (narg < 5 || narg > 6) 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);
int lj_type_one = find_lj_type(arg[2],lj_type_list);
if (lj_type_one == LJ_NOT_SET)
error->all(FLERR,"Cannot parse LJ type flag.");
double epsilon_one = force->numeric(arg[3]);
double sigma_one = force->numeric(arg[4]);
double cut_lj_one = cut_lj_global;
if (narg == 6) cut_lj_one = force->numeric(arg[5]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
lj_type[i][j] = lj_type_one;
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut_lj[i][j] = cut_lj_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q");
neighbor->request(this);
cut_coulsq = cut_coul * cut_coul;
// insure use of KSpace long-range solver, set g_ewald
if (force->kspace == NULL)
error->all(FLERR,"Pair style is incompatible with KSpace style");
g_ewald = force->kspace->g_ewald;
// setup force tables
if (ncoultablebits) init_tables();
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJSDKCoulLong::init_one(int i, int j)
{
if (setflag[i][j] == 0)
error->all(FLERR,"No mixing support for lj/sdk/coul/long. "
"Coefficients for all pairs need to be set explicitly.");
const int ljt = lj_type[i][j];
if (ljt == LJ_NOT_SET)
error->all(FLERR,"unrecognized LJ parameter flag");
double cut = MAX(cut_lj[i][j],cut_coul);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
lj1[i][j] = lj_prefact[ljt] * lj_pow1[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow1[ljt]);
lj2[i][j] = lj_prefact[ljt] * lj_pow2[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow2[ljt]);
lj3[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow1[ljt]);
lj4[i][j] = lj_prefact[ljt] * epsilon[i][j] * pow(sigma[i][j],lj_pow2[ljt]);
if (offset_flag) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = lj_prefact[ljt] * epsilon[i][j] * (pow(ratio,lj_pow1[ljt]) - pow(ratio,lj_pow2[ljt]));
} else offset[i][j] = 0.0;
cut_ljsq[j][i] = cut_ljsq[i][j];
cut_lj[j][i] = cut_lj[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];
offset[j][i] = offset[i][j];
lj_type[j][i] = lj_type[i][j];
// compute LJ derived parameters for SDK angle potential (LJ only!)
const double eps = epsilon[i][j];
const double sig = sigma[i][j];
const double rmin = sig*exp(1.0/(lj_pow1[ljt]-lj_pow2[ljt])
*log(lj_pow1[ljt]/lj_pow2[ljt]) );
rminsq[j][i] = rminsq[i][j] = rmin*rmin;
const double ratio = sig/rmin;
const double emin_one = lj_prefact[ljt] * eps * (pow(ratio,lj_pow1[ljt])
- pow(ratio,lj_pow2[ljt]));
emin[j][i] = emin[i][j] = emin_one;
// compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce
if (tail_flag)
error->all(FLERR,"Tail flag not supported by lj/sdk/coul/long pair style");
return cut;
}
/* ----------------------------------------------------------------------
setup force tables used in compute routines
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::init_tables()
{
int masklo,maskhi;
double r,grij,expm2,derfc;
double qqrd2e = force->qqrd2e;
tabinnersq = tabinner*tabinner;
init_bitmap(tabinner,cut_coul,ncoultablebits,
masklo,maskhi,ncoulmask,ncoulshiftbits);
int ntable = 1;
for (int i = 0; i < ncoultablebits; i++) ntable *= 2;
// linear lookup tables of length N = 2^ncoultablebits
// stored value = value at lower edge of bin
// d values = delta from lower edge to upper edge of bin
if (ftable) free_tables();
memory->create(rtable,ntable,"pair:rtable");
memory->create(ftable,ntable,"pair:ftable");
memory->create(ctable,ntable,"pair:ctable");
memory->create(etable,ntable,"pair:etable");
memory->create(drtable,ntable,"pair:drtable");
memory->create(dftable,ntable,"pair:dftable");
memory->create(dctable,ntable,"pair:dctable");
memory->create(detable,ntable,"pair:detable");
union_int_float_t rsq_lookup;
union_int_float_t minrsq_lookup;
int itablemin;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
for (int i = 0; i < ntablem1; i++) {
drtable[i] = 1.0/(rtable[i+1] - rtable[i]);
dftable[i] = ftable[i+1] - ftable[i];
dctable[i] = ctable[i+1] - ctable[i];
detable[i] = etable[i+1] - etable[i];
}
// get the delta values for the last table entries
// tables are connected periodically between 0 and ntablem1
drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]);
dftable[ntablem1] = ftable[0] - ftable[ntablem1];
dctable[ntablem1] = ctable[0] - ctable[ntablem1];
detable[ntablem1] = etable[0] - etable[ntablem1];
// get the correct delta values at itablemax
// smallest r is in bin itablemin
// largest r is in bin itablemax, which is itablemin-1,
// or ntablem1 if itablemin=0
// deltas at itablemax only needed if corresponding rsq < cut*cut
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq_lookup.f < cut_coulsq) {
rsq_lookup.f = cut_coulsq;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
c_tmp = qqrd2e/r;
e_tmp = qqrd2e/r * derfc;
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::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(&lj_type[i][j],sizeof(int),1,fp);
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::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) {
fread(&lj_type[i][j],sizeof(int),1,fp);
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut_lj[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&lj_type[i][j],1,MPI_INT,0,world);
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
free memory for tables used in pair computations
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::free_tables()
{
memory->destroy(rtable);
memory->destroy(drtable);
memory->destroy(ftable);
memory->destroy(dftable);
memory->destroy(ctable);
memory->destroy(dctable);
memory->destroy(etable);
memory->destroy(detable);
}
/* ---------------------------------------------------------------------- */
double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype,
double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,r,grij,expm2,t,erfc,prefactor;
double fraction,table,forcecoul,forcelj,phicoul,philj;
int itable;
r2inv = 1.0/rsq;
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
union_int_float_t rsq_lookup_single;
rsq_lookup_single.f = rsq;
itable = rsq_lookup_single.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = atom->q[i]*atom->q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = atom->q[i]*atom->q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
}
} else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
const int ljt = lj_type[itype][jtype];
const double ljpow1 = lj_pow1[ljt];
const double ljpow2 = lj_pow2[ljt];
const double ljpref = lj_prefact[ljt];
const double ratio = sigma[itype][jtype]/sqrt(rsq);
const double eps = epsilon[itype][jtype];
fforce = factor_lj * ljpref*eps * (ljpow1*pow(ratio,ljpow1)
- ljpow2*pow(ratio,ljpow2))/rsq;
philj = factor_lj * (ljpref*eps * (pow(ratio,ljpow1) - pow(ratio,ljpow2))
- offset[itype][jtype]);
} else fforce=0.0;
fforce = (forcecoul + factor_lj*forcelj) * r2inv;
double eng = 0.0;
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq)
phicoul = prefactor*erfc;
else {
table = etable[itable] + fraction*detable[itable];
phicoul = atom->q[i]*atom->q[j] * table;
}
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
eng += phicoul;
}
if (rsq < cut_ljsq[itype][jtype])
eng += philj;
return eng;
}
/* ---------------------------------------------------------------------- */
void *PairLJSDKCoulLong::extract(char *str, int &dim)
{
dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
if (strcmp(str,"sigma") == 0) return (void *) sigma;
if (strcmp(str,"lj_type") == 0) return (void *) lj_type;
if (strcmp(str,"lj1") == 0) return (void *) lj1;
if (strcmp(str,"lj2") == 0) return (void *) lj2;
if (strcmp(str,"lj3") == 0) return (void *) lj3;
if (strcmp(str,"lj4") == 0) return (void *) lj4;
if (strcmp(str,"rminsq") == 0) return (void *) rminsq;
if (strcmp(str,"emin") == 0) return (void *) emin;
dim = 0;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
return NULL;
}
/* ---------------------------------------------------------------------- */
double PairLJSDKCoulLong::memory_usage()
{
double bytes = Pair::memory_usage();
int n = atom->ntypes;
// setflag/lj_type
bytes += 2 * (n+1)*(n+1)*sizeof(int);
// lj_cut/lj_cutsq/epsilon/sigma/offset/lj1/lj2/lj3/lj4/rminsq/emin
bytes += 11 * (n+1)*(n+1)*sizeof(double);
if (ncoultablebits) {
int ntable = 1<<ncoultablebits;
bytes += 8 * ntable*sizeof(double);
}
return bytes;
}

82
src/USER-CG-CMM/pair_lj_sdk_coul_long.h Normal file
View File

@ -0,0 +1,82 @@
/* -*- c++ -*- ----------------------------------------------------------
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 author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lj/sdk/coul/long,PairLJSDKCoulLong)
PairStyle(cg/cmm/coul/long,PairLJSDKCoulLong)
#else
#ifndef LMP_PAIR_LJ_SDK_COUL_LONG_H
#define LMP_PAIR_LJ_SDK_COUL_LONG_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJSDKCoulLong : public Pair {
public:
PairLJSDKCoulLong(class LAMMPS *);
virtual ~PairLJSDKCoulLong();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
void *extract(char *, int &);
virtual double memory_usage();
protected:
double **cut_lj,**cut_ljsq;
double cut_coul,cut_coulsq;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset;
int **lj_type;
// cutoff and offset for minimum of LJ potential
// to be used in SDK angle potential, which
// uses only the repulsive part of the potential
double **rminsq, **emin;
double cut_lj_global;
double g_ewald;
double tabinnersq;
double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable;
double *etable,*detable;
int ncoulshiftbits,ncoulmask;
void allocate();
void init_tables();
void free_tables();
private:
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval();
};
}
#endif
#endif