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

Browse Source
This commit is contained in:
sjplimp
2008-02-08 21:57:06 +00:00
parent bd89b8239f
commit 0a5ca8126b
17 changed files with 2768 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/Makefile
View File

@ -16,7 +16,7 @@ OBJ = $(SRC:.cpp=.o)
PACKAGE = asphere class2 colloid dipole dpd granular \
kspace manybody meam molecule opt poems xtc
PACKUSER = user-ackland user-ewaldn
PACKUSER = user-ackland user-cg-cmm user-ewaldn
PACKALL = $(PACKAGE) $(PACKUSER)

43
src/USER-CG-CMM/Install.csh Normal file
View File

@ -0,0 +1,43 @@
# Install/unInstall package classes in LAMMPS
if ($1 == 1) then
cp style_user_cg_cmm.h ..
cp angle_cg_cmm.h ..
cp angle_cg_cmm.cpp ..
cp cg_cmm_parms.h ..
cp cg_cmm_parms.cpp ..
cp pair_cmm_common.h ..
cp pair_cmm_common.cpp ..
cp pair_cg_cmm.cpp ..
cp pair_cg_cmm.h ..
cp pair_cg_cmm_coul_cut.cpp ..
cp pair_cg_cmm_coul_cut.h ..
cp pair_cg_cmm_coul_long.cpp ..
cp pair_cg_cmm_coul_long.h ..
else if ($1 == 0) then
rm ../style_user_cg_cmm.h
touch ../style_user_cg_cmm.h
rm ../angle_cg_cmm.h
rm ../angle_cg_cmm.cpp
rm ../cg_cmm_parms.h
rm ../cg_cmm_parms.cpp
rm ../pair_cmm_common.h
rm ../pair_cmm_common.cpp
rm ../pair_cg_cmm.cpp
rm ../pair_cg_cmm.h
rm ../pair_cg_cmm_coul_cut.cpp
rm ../pair_cg_cmm_coul_cut.h
rm ../pair_cg_cmm_coul_long.cpp
rm ../pair_cg_cmm_coul_long.h
endif

37
src/USER-CG-CMM/README Normal file
View File

@ -0,0 +1,37 @@
The files in this directory are a user-contributed package for LAMMPS.
The person who created these files is Axel Kohlmeyer
(akohlmey@cmm.chem.upenn.edu). Contact him directly if you have
questions.
The current version of this package should be considered beta
quality. The CG potentials work correctly and well, but there will be
optimizations, cleanups and additional tools to aid in setting up and
analyzing simulations with this package added in the next months.
This package implements 4 commands which can be used in a LAMMPS input
script:
pair_style cg/cmm
pair_style cg/cmm/coul/cut
pair_style cg/cmm/coul/long
and angle_style cg/cmm.
See the documentation files for these commands for details.
There is also an cg-cmm example directory under examples/USER with
sample inputs and outputs.
These styles allow coarse grained MD simulations with the
parametrization of Shinoda, DeVane, Klein, Mol Sim, 33, 27 (2007)
(cg/cmm), with extensions to simulate ionic liquids, electrolytes,
lipids and charged amino acids (to be published soon).
Thanks for contributions, support and testing goes to
Wataru Shinoda (AIST, Tsukuba)
Russell DeVane (CMM / U Penn, Philadelphia)
Balasubramanian Sundaram (JNCASR, Bangalore)
Michael L. Klein (CMM / U Penn, Philadelphia)
version: 0.98 / 2008-01-31

349
src/USER-CG-CMM/angle_cg_cmm.cpp Normal file
View File

@ -0,0 +1,349 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Special Angle Potential for the CMM coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "angle_cg_cmm.h"
#include "atom.h"
#include "neighbor.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define SMALL 0.001
/* ---------------------------------------------------------------------- */
AngleCGCMM::AngleCGCMM(LAMMPS *lmp) : Angle(lmp) {}
/* ---------------------------------------------------------------------- */
AngleCGCMM::~AngleCGCMM()
{
if (allocated) {
memory->sfree(setflag);
memory->sfree(k);
memory->sfree(theta0);
memory->sfree(cg_type);
memory->sfree(epsilon);
memory->sfree(sigma);
memory->sfree(rcut);
}
}
/* ---------------------------------------------------------------------- */
void AngleCGCMM::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,r3,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,
// so this has to be done here and not in the
// general non-bonded code.
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;
r3 = sqrt(rsq3);
f13=0.0;
e13=0.0;
if (r3 < rcut[type]) {
const int cgt = cg_type[type];
const double cgpow1 = cg_pow1[cgt];
const double cgpow2 = cg_pow2[cgt];
const double cgpref = cg_prefact[cgt];
const double ratio = sigma[type]/r3;
const double eps = epsilon[type];
f13 = cgpref*eps / rsq3 * (cgpow1*pow(ratio,cgpow1)
- cgpow2*pow(ratio,cgpow2));
if (eflag) e13 = eps + cgpref*eps * (pow(ratio,cgpow1)
- pow(ratio,cgpow2));
}
// force & energy
dtheta = acos(c) - theta0[type];
tk = k[type] * dtheta;
if (eflag) eangle = tk*dtheta + e13;
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 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;
}
// FIXME: we'll need to override this one for accurate virial.
// for the time being, we hope that this contribution
// is very small (it should be).
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
delx1,dely1,delz1,delx2,dely2,delz2);
}
}
/* ---------------------------------------------------------------------- */
void AngleCGCMM::allocate()
{
allocated = 1;
int n = atom->nangletypes;
k = (double *) memory->smalloc((n+1)*sizeof(double),"angle:k");
theta0 = (double *) memory->smalloc((n+1)*sizeof(double),"angle:theta0");
epsilon = (double *) memory->smalloc((n+1)*sizeof(double),"angle:epsilon");
sigma = (double *) memory->smalloc((n+1)*sizeof(double),"angle:sigma");
rcut = (double *) memory->smalloc((n+1)*sizeof(double),"angle:rcut");
cg_type = (int *) memory->smalloc((n+1)*sizeof(int),"angle:cg_type");
setflag = (int *) memory->smalloc((n+1)*sizeof(int),"angle:setflag");
for (int i = 1; i <= n; i++) {
cg_type[i] = CG_NOT_SET;
setflag[i] = 0;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
void AngleCGCMM::coeff(int which, int narg, char **arg)
{
if (which > 0) return;
if (narg != 6) error->all("Incorrect args for angle coefficients");
if (!allocated) allocate();
int ilo,ihi;
force->bounds(arg[0],atom->nangletypes,ilo,ihi);
double k_one = atof(arg[1]);
double theta0_one = atof(arg[2]);
int cg_type_one=find_cg_type(arg[3]);
if (cg_type_one == CG_NOT_SET) error->all("Error reading CG type flag.");
double epsilon_one = atof(arg[4]);
double sigma_one = atof(arg[5]);
// find minimum of LJ potential. we only want to include
// the repulsive part of the 1-3 LJ.
double rcut_one = sigma_one*exp(
1.0/(cg_pow1[cg_type_one]-cg_pow2[cg_type_one])
*log(cg_pow1[cg_type_one]/cg_pow2[cg_type_one])
);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
// convert theta0 from degrees to radians
theta0[i] = theta0_one/180.0 * PI;
epsilon[i] = epsilon_one;
sigma[i] = sigma_one;
rcut[i] = rcut_one;
cg_type[i] = cg_type_one;
setflag[i] = 1;
count++;
}
if (count == 0) error->all("Incorrect args for angle coefficients");
}
/* ---------------------------------------------------------------------- */
double AngleCGCMM::equilibrium_angle(int i)
{
return theta0[i];
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void AngleCGCMM::write_restart(FILE *fp)
{
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
fwrite(&epsilon[1],sizeof(double),atom->nangletypes,fp);
fwrite(&sigma[1],sizeof(double),atom->nangletypes,fp);
fwrite(&rcut[1],sizeof(double),atom->nangletypes,fp);
fwrite(&cg_type[1],sizeof(int),atom->nangletypes,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void AngleCGCMM::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(&epsilon[1],sizeof(double),atom->nangletypes,fp);
fread(&sigma[1],sizeof(double),atom->nangletypes,fp);
fread(&rcut[1],sizeof(double),atom->nangletypes,fp);
fread(&cg_type[1],sizeof(int),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(&epsilon[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&rcut[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&cg_type[1],atom->nangletypes,MPI_INT,0,world);
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
}
/* ---------------------------------------------------------------------- */
double AngleCGCMM::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;
// 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 double r3 = sqrt(delx3*delx3 + dely3*dely3 + delz3*delz3);
double e13=0.0;
if (r3 < rcut[type]) {
const int cgt = cg_type[type];
const double cgpow1 = cg_pow1[cgt];
const double cgpow2 = cg_pow2[cgt];
const double cgpref = cg_prefact[cgt];
const double ratio = sigma[type]/r3;
const double eps = epsilon[type];
e13 = eps + cgpref*eps * (pow(ratio,cgpow1)
- pow(ratio,cgpow2));
}
double dtheta = acos(c) - theta0[type];
double tk = k[type] * dtheta;
return tk*dtheta + e13;
}

47
src/USER-CG-CMM/angle_cg_cmm.h Normal file
View File

@ -0,0 +1,47 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Special Angle Potential for the CMM coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#ifndef ANGLE_CG_CMM_H
#define ANGLE_CG_CMM_H
#include "stdio.h"
#include "angle.h"
#include "cg_cmm_parms.h"
namespace LAMMPS_NS {
class AngleCGCMM : public Angle, public CGCMMParms {
public:
AngleCGCMM(class LAMMPS *);
~AngleCGCMM();
void compute(int, int);
void coeff(int, int, char **);
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
private:
double *k,*theta0;
int *cg_type;
double *epsilon, *sigma, *rcut;
void allocate();
};
}
#endif

40
src/USER-CG-CMM/cg_cmm_parms.cpp Normal file
View File

@ -0,0 +1,40 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Common parameters for the CMM coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#include "cg_cmm_parms.h"
#include "string.h"
using namespace LAMMPS_NS;
/* static constant class members */
const char * const CGCMMParms::cg_type_list[] = {"none", "lj9_6", "lj12_4", "lj12_6"};
const double CGCMMParms::cg_prefact[] = {0.0, 6.75, 2.59807621135332, 4.0};
const double CGCMMParms::cg_pow1[] = {0.0, 9.0, 12.0, 12.0};
const double CGCMMParms::cg_pow2[] = {0.0, 6.0, 4.0, 6.0};
/* ---------------------------------------------------------------------- */
int CGCMMParms::find_cg_type(const char *label)
{
for (int i=0; i < NUM_CG_TYPES; ++i) {
if (strcmp(label,cg_type_list[i]) == 0) {
return i;
}
}
return CG_NOT_SET;
}
/* ---------------------------------------------------------------------- */

40
src/USER-CG-CMM/cg_cmm_parms.h Normal file
View File

@ -0,0 +1,40 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Common parameters for the CMM coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#ifndef CG_CMM_PARMS_H
#define CG_CMM_PARMS_H
namespace LAMMPS_NS {
class CGCMMParms {
public:
// CG type flags. list of supported LJ exponent combinations
enum {CG_NOT_SET=0, CG_LJ9_6, CG_LJ12_4, CG_LJ12_6, NUM_CG_TYPES,
CG_COUL_NONE, CG_COUL_CUT, CG_COUL_DEBYE, CG_COUL_LONG};
int find_cg_type(const char *);
protected:
// coarse grain flags
static const char * const cg_type_list[NUM_CG_TYPES];
static const double cg_prefact[NUM_CG_TYPES];
static const double cg_pow1[NUM_CG_TYPES] ;
static const double cg_pow2[NUM_CG_TYPES] ;
};
}
#endif

171
src/USER-CG-CMM/pair_cg_cmm.cpp Normal file
View File

@ -0,0 +1,171 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
CMM coarse grained MD potentials. Plain version w/o charges.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#include "pair_cg_cmm.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairCGCMM::PairCGCMM(LAMMPS *lmp) : PairCMMCommon(lmp)
{
respa_enable = 0;
single_enable = 0;
}
/* ---------------------------------------------------------------------- */
PairCGCMM::~PairCGCMM()
{
/* empty */ ;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- *
* the real compute work is done in the PairCMMCommon::eval_XXX<>() templates
* in the common PairCG class. Through using templates we can have one
* implementation for all CG varieties _and_ gain speed through having
* the compiler optimize away conditionals within the innerloops that
* can be predetermined outside the loop through instantiation of the
* different combination of template flags.
* ---------------------------------------------------------------------- */
void PairCGCMM::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) {
return eval_verlet<1,1,1,CG_COUL_NONE>();
} else {
return eval_verlet<1,1,0,CG_COUL_NONE>();
}
} else {
if (force->newton_pair) {
return eval_verlet<1,0,1,CG_COUL_NONE>();
} else {
return eval_verlet<1,0,0,CG_COUL_NONE>();
}
}
} else {
if (force->newton_pair) {
return eval_verlet<0,0,1,CG_COUL_NONE>();
} else {
return eval_verlet<0,0,0,CG_COUL_NONE>();
}
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMM::compute_inner()
{
if (force->newton_pair) {
return eval_inner<1,CG_COUL_NONE>();
} else {
return eval_inner<0,CG_COUL_NONE>();
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMM::compute_middle()
{
if (force->newton_pair) {
return eval_middle<1,CG_COUL_NONE>();
} else {
return eval_middle<0,CG_COUL_NONE>();
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMM::compute_outer(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else {
evflag = 0;
}
if (evflag) {
if (eflag) {
if (vflag) {
if (force->newton_pair) {
return eval_outer<1,1,1,1,CG_COUL_NONE>();
} else {
return eval_outer<1,1,1,0,CG_COUL_NONE>();
}
} else {
if (force->newton_pair) {
return eval_outer<1,1,0,1,CG_COUL_NONE>();
} else {
return eval_outer<1,1,0,0,CG_COUL_NONE>();
}
}
} else {
if (vflag) {
if (force->newton_pair) {
return eval_outer<1,0,1,1,CG_COUL_NONE>();
} else {
return eval_outer<1,0,1,0,CG_COUL_NONE>();
}
} else {
if (force->newton_pair) {
return eval_outer<1,0,0,1,CG_COUL_NONE>();
} else {
return eval_outer<1,0,0,0,CG_COUL_NONE>();
}
}
}
} else {
if (force->newton_pair) {
return eval_outer<0,0,0,1,CG_COUL_NONE>();
} else {
return eval_outer<0,0,0,0,CG_COUL_NONE>();
}
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMM::write_restart(FILE *fp)
{
write_restart_settings(fp);
PairCMMCommon::write_restart(fp);
}
/* ---------------------------------------------------------------------- */
void PairCGCMM::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
PairCMMCommon::read_restart(fp);
}
/* ---------------------------------------------------------------------- */
double PairCGCMM::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj, double &fforce)
{
return eval_single(CG_COUL_NONE,i,j,itype,jtype,rsq,factor_coul,factor_lj,fforce);
}

41
src/USER-CG-CMM/pair_cg_cmm.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
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef PAIR_CG_CMM_H
#define PAIR_CG_CMM_H
#include "pair_cmm_common.h"
namespace LAMMPS_NS {
class PairCGCMM : public PairCMMCommon {
public:
PairCGCMM(class LAMMPS *);
virtual ~PairCGCMM();
void compute(int, int);
void compute_inner();
void compute_middle();
void compute_outer(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int, double, double, double, double &);
private:
// disable default destructor
PairCGCMM();
};
}
#endif

234
src/USER-CG-CMM/pair_cg_cmm_coul_cut.cpp Normal file
View File

@ -0,0 +1,234 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
CMM coarse grained MD potentials. Coulomb with cutoff version.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#include "pair_cg_cmm_coul_cut.h"
#include "memory.h"
#include "atom.h"
#include "string.h"
#define MIN(a,b) ((a) < (b) ? (a) : (b))
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairCGCMMCoulCut::PairCGCMMCoulCut(LAMMPS *lmp) : PairCMMCommon(lmp)
{
respa_enable = 0;
single_enable = 0;
}
/* ---------------------------------------------------------------------- */
PairCGCMMCoulCut::~PairCGCMMCoulCut()
{
if (allocated_coul) {
memory->destroy_2d_double_array(cut_lj);
memory->destroy_2d_double_array(cut_ljsq);
memory->destroy_2d_double_array(cut_coul);
memory->destroy_2d_double_array(cut_coulsq);
allocated_coul=0;
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::allocate()
{
PairCMMCommon::allocate();
allocated_coul = 1;
int n = atom->ntypes;
cut_lj = memory->create_2d_double_array(n+1,n+1,"paircg:cut_lj");
cut_ljsq = memory->create_2d_double_array(n+1,n+1,"paircg:cut_ljsq");
cut_coul = memory->create_2d_double_array(n+1,n+1,"paircg:cut_coul");
cut_coulsq = memory->create_2d_double_array(n+1,n+1,"paircg:cut_coulsq");
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::init_style()
{
if (!atom->q_flag)
error->all("Pair style cg/cut/coul/cut requires atom attribute q");
PairCMMCommon::init_style();
// set & error check interior rRESPA cutoffs
if (strcmp(update->integrate_style,"respa") == 0) {
if (((Respa *) update->integrate)->level_inner >= 0) {
cut_respa = ((Respa *) update->integrate)->cutoff;
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
if (MIN(cut_lj[i][j],cut_coul[i][j]) < cut_respa[3])
error->all("Pair cutoff < Respa interior cutoff");
}
} else cut_respa = NULL;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- *
* the real compute work is done in the PairCMMCommon::eval_XXX<>() templates
* in the common PairCG class. Through using templates we can have one
* implementation for all CG varieties _and_ gain speed through having
* the compiler optimize away conditionals within the innerloops that
* can be predetermined outside the loop through instantiation of the
* different combination of template flags.
* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::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) {
return eval_verlet<1,1,1,CG_COUL_CUT>();
} else {
return eval_verlet<1,1,0,CG_COUL_CUT>();
}
} else {
if (force->newton_pair) {
return eval_verlet<1,0,1,CG_COUL_CUT>();
} else {
return eval_verlet<1,0,0,CG_COUL_CUT>();
}
}
} else {
if (force->newton_pair) {
return eval_verlet<0,0,1,CG_COUL_CUT>();
} else {
return eval_verlet<0,0,0,CG_COUL_CUT>();
}
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::compute_inner()
{
if (force->newton_pair) {
return eval_inner<1,CG_COUL_CUT>();
} else {
return eval_inner<0,CG_COUL_CUT>();
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::compute_middle()
{
if (force->newton_pair) {
return eval_middle<1,CG_COUL_CUT>();
} else {
return eval_middle<0,CG_COUL_CUT>();
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::compute_outer(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else {
evflag = 0;
}
if (evflag) {
if (eflag) {
if (vflag) {
if (force->newton_pair) {
return eval_outer<1,1,1,1,CG_COUL_CUT>();
} else {
return eval_outer<1,1,1,0,CG_COUL_CUT>();
}
} else {
if (force->newton_pair) {
return eval_outer<1,1,0,1,CG_COUL_CUT>();
} else {
return eval_outer<1,1,0,0,CG_COUL_CUT>();
}
}
} else {
if (vflag) {
if (force->newton_pair) {
return eval_outer<1,0,1,1,CG_COUL_CUT>();
} else {
return eval_outer<1,0,1,0,CG_COUL_CUT>();
}
} else {
if (force->newton_pair) {
return eval_outer<1,0,0,1,CG_COUL_CUT>();
} else {
return eval_outer<1,0,0,0,CG_COUL_CUT>();
}
}
}
} else {
if (force->newton_pair) {
return eval_outer<0,0,0,1,CG_COUL_CUT>();
} else {
return eval_outer<0,0,0,0,CG_COUL_CUT>();
}
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::write_restart(FILE *fp)
{
write_restart_settings(fp);
PairCMMCommon::write_restart(fp);
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulCut::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
PairCMMCommon::read_restart(fp);
}
/* ---------------------------------------------------------------------- */
double PairCGCMMCoulCut::memory_usage()
{
double bytes=PairCMMCommon::memory_usage();
int n = atom->ntypes;
// cut_coul/cut_coulsq/cut_lj/cut_ljsq;
bytes += (n+1)*(n+1)*sizeof(double)*4;
return bytes;
}
/* ---------------------------------------------------------------------- */
double PairCGCMMCoulCut::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj, double &fforce)
{
return eval_single(CG_COUL_CUT,i,j,itype,jtype,rsq,factor_coul,factor_lj,fforce);
}

42
src/USER-CG-CMM/pair_cg_cmm_coul_cut.h Normal file
View File

@ -0,0 +1,42 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef PAIR_CG_CMM_COUL_CUT_H
#define PAIR_CG_CMM_COUL_CUT_H
#include "pair_cmm_common.h"
namespace LAMMPS_NS {
class PairCGCMMCoulCut : public PairCMMCommon {
public:
PairCGCMMCoulCut(class LAMMPS *);
~PairCGCMMCoulCut();
void compute(int, int);
void compute_inner();
void compute_middle();
void compute_outer(int, int);
void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
double memory_usage();
double single(int, int, int, int, double, double, double, double &);
protected:
void allocate();
};
}
#endif

446
src/USER-CG-CMM/pair_cg_cmm_coul_long.cpp Normal file
View File

@ -0,0 +1,446 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
CMM coarse grained MD potentials. Coulomb with k-space version.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#include "pair_cg_cmm_coul_long.h"
#include "memory.h"
#include "atom.h"
#include "force.h"
#include "kspace.h"
#include "string.h"
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define EWALD_F 1.12837917
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairCGCMMCoulLong::PairCGCMMCoulLong(LAMMPS *lmp) : PairCMMCommon(lmp)
{
respa_enable = 0;
single_enable = 0;
}
/* ---------------------------------------------------------------------- */
PairCGCMMCoulLong::~PairCGCMMCoulLong()
{
if (allocated_coul) {
memory->destroy_2d_double_array(cut_lj);
memory->destroy_2d_double_array(cut_ljsq);
memory->destroy_2d_double_array(cut_coul);
memory->destroy_2d_double_array(cut_coulsq);
allocated_coul=0;
}
if (ftable) free_tables();
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::allocate()
{
PairCMMCommon::allocate();
allocated_coul = 1;
int n = atom->ntypes;
cut_lj = memory->create_2d_double_array(n+1,n+1,"paircg:cut_lj");
cut_ljsq = memory->create_2d_double_array(n+1,n+1,"paircg:cut_ljsq");
cut_coul = memory->create_2d_double_array(n+1,n+1,"paircg:cut_coul");
cut_coulsq = memory->create_2d_double_array(n+1,n+1,"paircg:cut_coulsq");
}
/* ----------------------------------------------------------------------
free memory for tables used in pair computations
------------------------------------------------------------------------- */
void PairCGCMMCoulLong::free_tables()
{
memory->sfree(rtable);
memory->sfree(drtable);
memory->sfree(ftable);
memory->sfree(dftable);
memory->sfree(ctable);
memory->sfree(dctable);
memory->sfree(etable);
memory->sfree(detable);
memory->sfree(vtable);
memory->sfree(dvtable);
memory->sfree(ptable);
memory->sfree(dptable);
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::init_style()
{
if (!atom->q_flag)
error->all("Pair style cg/cut/coul/long requires atom attribute q");
PairCMMCommon::init_style();
// set & error check interior rRESPA cutoffs
if (strcmp(update->integrate_style,"respa") == 0) {
if (((Respa *) update->integrate)->level_inner >= 0) {
cut_respa = ((Respa *) update->integrate)->cutoff;
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
if (MIN(cut_lj[i][j],cut_coul_global) < cut_respa[3])
error->all("Pair cutoff < Respa interior cutoff");
}
} else cut_respa = NULL;
// ensure use of KSpace long-range solver, set g_ewald
if (force->kspace == NULL)
error->all("Pair style is incompatible with KSpace style");
g_ewald = force->kspace->g_ewald;
// setup force tables
if (ncoultablebits) init_tables();
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::init_tables()
{
int masklo,maskhi;
double r,grij,expm2,derfc,rsw;
double qqrd2e = force->qqrd2e;
tabinnersq = tabinner*tabinner;
init_bitmap(tabinner,cut_coul_global,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();
rtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:rtable");
ftable = (double *) memory->smalloc(ntable*sizeof(double),"pair:ftable");
ctable = (double *) memory->smalloc(ntable*sizeof(double),"pair:ctable");
etable = (double *) memory->smalloc(ntable*sizeof(double),"pair:etable");
drtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:drtable");
dftable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dftable");
dctable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dctable");
detable = (double *) memory->smalloc(ntable*sizeof(double),"pair:detable");
if (cut_respa == NULL) {
vtable = ptable = dvtable = dptable = NULL;
} else {
vtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:vtable");
ptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:ptable");
dvtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dvtable");
dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable");
}
float rsq;
int *int_rsq = (int *) &rsq;
float minrsq;
int *int_minrsq = (int *) &minrsq;
int itablemin;
*int_minrsq = 0 << ncoulshiftbits;
*int_minrsq = *int_minrsq | maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tabinnersq) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
}
r = sqrtf(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
} else {
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
}
}
}
minrsq = MIN(minrsq,rsq);
}
tabinnersq = minrsq;
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];
}
if (cut_respa) {
for (int i = 0; i < ntablem1; i++) {
dvtable[i] = vtable[i+1] - vtable[i];
dptable[i] = ptable[i+1] - ptable[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];
if (cut_respa) {
dvtable[ntablem1] = vtable[0] - vtable[ntablem1];
dptable[ntablem1] = ptable[0] - ptable[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,p_tmp,v_tmp;
itablemin = *int_minrsq & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
*int_rsq = itablemax << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
if (rsq < cut_coulsq_global) {
rsq = cut_coulsq_global;
r = sqrtf(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
c_tmp = qqrd2e/r;
e_tmp = qqrd2e/r * derfc;
} else {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
c_tmp = 0.0;
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
} else {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
c_tmp = qqrd2e/r;
}
}
}
drtable[itablemax] = 1.0/(rsq - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
if (cut_respa) {
dvtable[itablemax] = v_tmp - vtable[itablemax];
dptable[itablemax] = p_tmp - ptable[itablemax];
}
}
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- *
* the real compute work is done in the PairCMMCommon::eval_XXX<>() templates
* in the common PairCG class. Through using templates we can have one
* implementation for all CG varieties _and_ gain speed through having
* the compiler optimize away conditionals within the innerloops that
* can be predetermined outside the loop through instantiation of the
* different combination of template flags.
* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::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) {
return eval_verlet<1,1,1,CG_COUL_LONG>();
} else {
return eval_verlet<1,1,0,CG_COUL_LONG>();
}
} else {
if (force->newton_pair) {
return eval_verlet<1,0,1,CG_COUL_LONG>();
} else {
return eval_verlet<1,0,0,CG_COUL_LONG>();
}
}
} else {
if (force->newton_pair) {
return eval_verlet<0,0,1,CG_COUL_LONG>();
} else {
return eval_verlet<0,0,0,CG_COUL_LONG>();
}
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::compute_inner()
{
if (force->newton_pair) {
return eval_inner<1,CG_COUL_LONG>();
} else {
return eval_inner<0,CG_COUL_LONG>();
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::compute_middle()
{
if (force->newton_pair) {
return eval_middle<1,CG_COUL_LONG>();
} else {
return eval_middle<0,CG_COUL_LONG>();
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::compute_outer(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else {
evflag = 0;
}
if (evflag) {
if (eflag) {
if (vflag) {
if (force->newton_pair) {
return eval_outer<1,1,1,1,CG_COUL_LONG>();
} else {
return eval_outer<1,1,1,0,CG_COUL_LONG>();
}
} else {
if (force->newton_pair) {
return eval_outer<1,1,0,1,CG_COUL_LONG>();
} else {
return eval_outer<1,1,0,0,CG_COUL_LONG>();
}
}
} else {
if (vflag) {
if (force->newton_pair) {
return eval_outer<1,0,1,1,CG_COUL_LONG>();
} else {
return eval_outer<1,0,1,0,CG_COUL_LONG>();
}
} else {
if (force->newton_pair) {
return eval_outer<1,0,0,1,CG_COUL_LONG>();
} else {
return eval_outer<1,0,0,0,CG_COUL_LONG>();
}
}
}
} else {
if (force->newton_pair) {
return eval_outer<0,0,0,1,CG_COUL_LONG>();
} else {
return eval_outer<0,0,0,0,CG_COUL_LONG>();
}
}
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::write_restart(FILE *fp)
{
write_restart_settings(fp);
PairCMMCommon::write_restart(fp);
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
PairCMMCommon::read_restart(fp);
}
/* ---------------------------------------------------------------------- */
double PairCGCMMCoulLong::memory_usage()
{
double bytes=PairCMMCommon::memory_usage();
int n = atom->ntypes;
// cut_coul/cut_coulsq/cut_ljsq
bytes += (n+1)*(n+1)*sizeof(double)*4;
return bytes;
}
/* ---------------------------------------------------------------------- */
double PairCGCMMCoulLong::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj, double &fforce)
{
return eval_single(CG_COUL_LONG,i,j,itype,jtype,rsq,factor_coul,factor_lj,fforce);
}
/* ---------------------------------------------------------------------- */
void *PairCGCMMCoulLong::extract(char *str)
{
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul_global;
return NULL;
}
/* ---------------------------------------------------------------------- */

45
src/USER-CG-CMM/pair_cg_cmm_coul_long.h Normal file
View File

@ -0,0 +1,45 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef PAIR_CG_CMM_COUL_LONG_H
#define PAIR_CG_CMM_COUL_LONG_H
#include "pair_cmm_common.h"
namespace LAMMPS_NS {
class PairCGCMMCoulLong : public PairCMMCommon {
public:
PairCGCMMCoulLong(class LAMMPS *);
~PairCGCMMCoulLong();
void compute(int, int);
void compute_inner();
void compute_middle();
void compute_outer(int, int);
void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
double memory_usage();
double single(int, int, int, int, double, double, double, double &);
void *extract(char *str);
protected:
void allocate();
void init_tables();
void free_tables();
};
}
#endif

477
src/USER-CG-CMM/pair_cmm_common.cpp Normal file
View File

@ -0,0 +1,477 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Common functionality for the CMM coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#include "pair_cmm_common.h"
#include "memory.h"
#include "stdlib.h"
#include "string.h"
#include "ctype.h"
#include "math.h"
using namespace LAMMPS_NS;
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#define SMALL 1.0e-6
/* ---------------------------------------------------------------------- */
PairCMMCommon::PairCMMCommon(class LAMMPS *lmp) : Pair(lmp)
{
ftable = NULL;
allocated_coul = 0;
kappa = 0.0;
respa_enable = 0;
single_enable = 0;
}
/* ---------------------------------------------------------------------- *
* clean up common arrays *
* ---------------------------------------------------------------------- */
PairCMMCommon::~PairCMMCommon()
{
if (allocated) {
memory->destroy_2d_int_array(setflag);
memory->destroy_2d_int_array(cg_type);
memory->destroy_2d_double_array(cut);
memory->destroy_2d_double_array(cutsq);
memory->destroy_2d_double_array(epsilon);
memory->destroy_2d_double_array(sigma);
memory->destroy_2d_double_array(offset);
memory->destroy_2d_double_array(lj1);
memory->destroy_2d_double_array(lj2);
memory->destroy_2d_double_array(lj3);
memory->destroy_2d_double_array(lj4);
allocated = 0;
}
}
/* ---------------------------------------------------------------------- *
* allocate common arrays *
* ---------------------------------------------------------------------- */
void PairCMMCommon::allocate()
{
allocated = 1;
int n = atom->ntypes;
setflag = memory->create_2d_int_array(n+1,n+1,"paircg:setflag");
cg_type = memory->create_2d_int_array(n+1,n+1,"paircg:cg_type");
for (int i = 1; i <= n; i++) {
for (int j = i; j <= n; j++) {
setflag[i][j] = 0;
cg_type[i][j] = CG_NOT_SET;
}
}
cut = memory->create_2d_double_array(n+1,n+1,"paircg:cut");
cutsq = memory->create_2d_double_array(n+1,n+1,"paircg:cutsq");
epsilon = memory->create_2d_double_array(n+1,n+1,"paircg:epsilon");
sigma = memory->create_2d_double_array(n+1,n+1,"paircg:sigma");
offset = memory->create_2d_double_array(n+1,n+1,"paircg:offset");
lj1 = memory->create_2d_double_array(n+1,n+1,"paircg:lj1");
lj2 = memory->create_2d_double_array(n+1,n+1,"paircg:lj2");
lj3 = memory->create_2d_double_array(n+1,n+1,"paircg:lj3");
lj4 = memory->create_2d_double_array(n+1,n+1,"paircg:lj4");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
// arguments to the pair_style command (global version)
// args = cutoff (cutoff2 (kappa))
void PairCMMCommon::settings(int narg, char **arg)
{
if ((narg < 1) || (narg > 3)) error->all("Illegal pair_style command");
cut_lj_global = atof(arg[0]);
if (narg == 1) cut_coul_global = cut_lj_global;
else cut_coul_global = atof(arg[1]);
cut_coulsq_global = cut_coul_global*cut_coul_global;
// exponential coulomb screening (optional)
if (narg == 3) kappa = atof(arg[2]);
if (fabs(kappa) < SMALL) kappa=0.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_lj_global;
if (allocated_coul) {
cut[i][j] = MAX(cut_lj_global,cut_coul_global);
cut_lj[i][j] = cut_lj_global;
cut_coul[i][j] = cut_coul_global;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairCMMCommon::coeff(int narg, char **arg)
{
if (narg < 5 || narg > 7) error->all("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 cg_type_one=find_cg_type(arg[2]);
if (cg_type_one == CG_NOT_SET) error->all("Error reading CG type flag.");
double epsilon_one = atof(arg[3]);
double sigma_one = atof(arg[4]);
double cut_lj_one = cut_lj_global;
double cut_coul_one = cut_coul_global;
if (narg >= 6) cut_lj_one = atof(arg[5]);
if (narg == 7) cut_coul_one = atof(arg[6]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cg_type[i][j] = cg_type_one;
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
setflag[i][j] = 1;
if (allocated_coul) {
cut_lj[i][j] = cut_lj_one;
cut_coul[i][j] = cut_coul_one;
} else {
cut[i][j] = cut_lj_one;
}
count++;
}
}
if (count == 0) error->all("Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairCMMCommon::init_style()
{
// request regular or rRESPA neighbor lists
int irequest;
if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) {
int respa = 0;
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
if (respa == 0) irequest = neighbor->request(this);
else if (respa == 1) {
irequest = neighbor->request(this);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
} else {
irequest = neighbor->request(this);
neighbor->requests[irequest]->id = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respainner = 1;
irequest = neighbor->request(this);
neighbor->requests[irequest]->id = 2;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respamiddle = 1;
irequest = neighbor->request(this);
neighbor->requests[irequest]->id = 3;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->respaouter = 1;
}
} else irequest = neighbor->request(this);
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
regular or rRESPA
------------------------------------------------------------------------- */
void PairCMMCommon::init_list(int id, NeighList *ptr)
{
if (id == 0) list = ptr;
else if (id == 1) listinner = ptr;
else if (id == 2) listmiddle = ptr;
else if (id == 3) listouter = ptr;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairCMMCommon::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
error->all("for CG styles, epsilon and sigma need to be set explicitly for all pairs.");
}
const int cgt = cg_type[i][j];
if (cgt == CG_NOT_SET)
error->all("unrecognized LJ parameter flag");
lj1[i][j] = cg_prefact[cgt] * cg_pow1[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow1[cgt]);
lj2[i][j] = cg_prefact[cgt] * cg_pow2[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow2[cgt]);
lj3[i][j] = cg_prefact[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow1[cgt]);
lj4[i][j] = cg_prefact[cgt] * epsilon[i][j] * pow(sigma[i][j],cg_pow2[cgt]);
double mycut = cut[i][j];
if (offset_flag) {
double ratio = sigma[i][j] / mycut;
offset[i][j] = cg_prefact[cgt] * epsilon[i][j] * (pow(ratio,cg_pow1[cgt]) - pow(ratio,cg_pow2[cgt]));
} else offset[i][j] = 0.0;
if (allocated_coul) {
mycut = MAX(cut_lj[i][j],cut_coul[i][j]);
cut_ljsq[i][j]=cut_lj[i][j]*cut_lj[i][j];
cut_coulsq[i][j]=cut_coul[i][j]*cut_coul[i][j];
if (offset_flag) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = cg_prefact[cgt] * epsilon[i][j] * (pow(ratio,cg_pow1[cgt]) - pow(ratio,cg_pow2[cgt]));
} else offset[i][j] = 0.0;
}
// make sure data is stored symmetrically
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];
cg_type[j][i] = cg_type[i][j];
if (allocated_coul) {
cut_lj[j][i]=cut_lj[i][j];
cut_ljsq[j][i]=cut_ljsq[i][j];
cut_coul[j][i]=cut_coul[i][j];
cut_coulsq[j][i]=cut_coulsq[i][j];
}
// set & error check interior rRESPA cutoff
if (strcmp(update->integrate_style,"respa") == 0) {
if (((Respa *) update->integrate)->level_inner >= 0) {
cut_respa = ((Respa *) update->integrate)->cutoff;
if (cut[i][j] < cut_respa[3])
error->all("Pair cutoff < Respa interior cutoff");
}
} else cut_respa = NULL;
// compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce
if (tail_flag) {
#if 1
error->all("tail correction not (yet) supported by CG potentials.");
#else
int *type = atom->type;
int nlocal = atom->nlocal;
double count[2],all[2];
count[0] = count[1] = 0.0;
for (int k = 0; k < nlocal; k++) {
if (type[k] == i) count[0] += 1.0;
if (type[k] == j) count[1] += 1.0;
}
MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double PI = 4.0*atan(1.0);
double sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
double rc6 = rc3*rc3;
double rc9 = rc3*rc6;
etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
ptail_ij = 16.0*PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
#endif
}
return mycut;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCMMCommon::write_restart(FILE *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(&cg_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);
if (allocated_coul) {
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
fwrite(&cut_coul[i][j],sizeof(double),1,fp);
}
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCMMCommon::read_restart(FILE *fp)
{
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(&cg_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);
if(allocated_coul) {
fread(&cut_lj[i][j],sizeof(double),1,fp);
fread(&cut_coul[i][j],sizeof(double),1,fp);
}
}
MPI_Bcast(&cg_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);
if (allocated_coul) {
MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);
}
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCMMCommon::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul_global,sizeof(double),1,fp);
fwrite(&kappa,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 PairCMMCommon::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul_global,sizeof(double),1,fp);
fread(&kappa,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_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&kappa,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairCMMCommon::memory_usage()
{
double bytes=Pair::memory_usage();
int n = atom->ntypes;
// setflag/cg_type
bytes += (n+1)*(n+1)*sizeof(int)*2;
// cut/cutsq/epsilon/sigma/offset/lj1/lj2/lj3/lj4
bytes += (n+1)*(n+1)*sizeof(double)*9;
return bytes;
}
/* ------------------------------------------------------------------------ */
double PairCMMCommon::eval_single(int coul_type, int i, int j, int itype, int jtype,
double rsq, double factor_coul, double factor_lj,
double &fforce)
{
double lj_force, lj_erg, coul_force, coul_erg;
lj_force=lj_erg=coul_force=coul_erg=0.0;
if (rsq < cut_ljsq[itype][jtype]) {
const int cgt = cg_type[itype][jtype];
const double cgpow1 = cg_pow1[cgt];
const double cgpow2 = cg_pow2[cgt];
const double cgpref = cg_prefact[cgt];
const double ratio = sigma[itype][jtype]/sqrt(rsq);
const double eps = epsilon[itype][jtype];
lj_force = cgpref*eps * rsq * (cgpow1*pow(ratio,cgpow1)
- cgpow2*pow(ratio,cgpow2));
lj_erg = cgpref*eps * (pow(ratio,cgpow1) - pow(ratio,cgpow2));
}
if (rsq < cut_coul[itype][jtype]) {
if(coul_type == CG_COUL_LONG) {
error->all("single energy computation with coulomb not supported by CG potentials.");
} else if ((coul_type == CG_COUL_CUT) || (coul_type == CG_COUL_DEBYE)) {
error->all("single energy computation with coulomb not supported by CG potentials.");
} else if (coul_type == CG_COUL_NONE) {
; // do nothing
} else {
error->all("unknown coulomb type with CG potentials.");
}
}
fforce = factor_lj*lj_force + factor_coul*coul_force;
return factor_lj*lj_erg + factor_coul*coul_erg;
}
/* ------------------------------------------------------------------------ */

718
src/USER-CG-CMM/pair_cmm_common.h Normal file
View File

@ -0,0 +1,718 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Common functionality for the CMM coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@cmm.chem.upenn.edu>
------------------------------------------------------------------------- */
#ifndef PAIR_CMM_COMMON_H
#define PAIR_CMM_COMMON_H
#include "pair.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "respa.h"
#include "update.h"
#include "cg_cmm_parms.h"
#include "math.h"
namespace LAMMPS_NS {
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define EWALD_A1 0.254829592
#define EWALD_A2 -0.284496736
#define EWALD_A3 1.421413741
#define EWALD_A4 -1.453152027
#define EWALD_A5 1.061405429
class PairCMMCommon : public Pair , public CGCMMParms {
public:
PairCMMCommon(class LAMMPS *);
virtual ~PairCMMCommon();
virtual void settings(int, char **);
virtual void coeff(int, char **);
virtual void init_style();
virtual void init_list(int, class NeighList *);
virtual double init_one(int, int);
virtual void write_restart(FILE *);
virtual void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual double memory_usage();
protected:
// coarse grain flags
int **cg_type;
// lennard jones parameters
double cut_lj_global, **cut, **cut_lj, **cut_ljsq;
double **epsilon, **sigma;
double **lj1, **lj2, **lj3, **lj4, **offset;
// coulomb parameters
int allocated_coul; // 0/1 = whether coulomb arrays are allocated
double cut_coul_global, cut_coulsq_global, kappa, g_ewald;
double **cut_coul, **cut_coulsq;
// tables
double tabinnersq;
double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable;
double *etable,*detable,*ptable,*dptable,*vtable,*dvtable;
int ncoulshiftbits,ncoulmask;
// r-RESPA parameters
double *cut_respa;
// methods
virtual void allocate();
private:
// disable default constructor
PairCMMCommon();
protected:
// general optimizeable real space loops
template < const int EVFLAG, const int EFLAG,
const int NEWTON_PAIR, const int COUL_TYPE >
void eval_verlet();
template < const int NEWTON_PAIR, const int COUL_TYPE >
void eval_inner();
template < const int NEWTON_PAIR, const int COUL_TYPE >
void eval_middle();
template < const int EVFLAG, const int EFLAG, const int VFLAG,
const int NEWTON_PAIR, const int COUL_TYPE >
void eval_outer();
// this one is not performance critical... no template needed.
double eval_single(int, int, int, int, int,
double, double, double, double &);
};
/* ---------------------------------------------------------------------- */
/* this is the inner heart of the CG potentials. */
#define CG_LJ_INNER(eflag,fvar) \
fvar=factor_lj; \
if (eflag) evdwl=factor_lj; \
\
if (cgt == CG_LJ12_4) { \
const double r4inv=r2inv*r2inv; \
\
fvar *= 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 (cgt == CG_LJ9_6) { \
const double r3inv = r2inv*sqrt(r2inv); \
const double r6inv = r3inv*r3inv; \
fvar *= r6inv*(lj1[itype][jtype]*r3inv \
- lj2[itype][jtype]); \
if (eflag) { \
evdwl *= r6inv*(lj3[itype][jtype]*r3inv \
- lj4[itype][jtype]) - offset[itype][jtype]; \
} \
} else if (cgt == CG_LJ12_6) { \
const double r6inv = r2inv*r2inv*r2inv; \
fvar *= r6inv*(lj1[itype][jtype]*r6inv \
- lj2[itype][jtype]); \
if (eflag) { \
evdwl *= r6inv*(lj3[itype][jtype]*r6inv \
- lj4[itype][jtype]) - offset[itype][jtype]; \
} \
} else { \
/* do nothing. this is a "cannot happen(TM)" case */ \
; \
}
#define CG_LJ_ENERGY(eflag) \
if (eflag) { \
evdwl=factor_lj; \
\
if (cgt == CG_LJ12_4) { \
const double r4inv=r2inv*r2inv; \
evdwl *= r4inv*(lj3[itype][jtype]*r4inv*r4inv \
- lj4[itype][jtype]) - offset[itype][jtype]; \
} else if (cgt == CG_LJ9_6) { \
const double r3inv = r2inv*sqrt(r2inv); \
const double r6inv = r3inv*r3inv; \
evdwl *= r6inv*(lj3[itype][jtype]*r3inv \
- lj4[itype][jtype]) - offset[itype][jtype]; \
} else if (cgt == CG_LJ12_6) { \
const double r6inv = r2inv*r2inv*r2inv; \
evdwl *= r6inv*(lj3[itype][jtype]*r6inv \
- lj4[itype][jtype]) - offset[itype][jtype]; \
} else { \
/* do nothing. this is a "cannot happen(TM)" case */ \
; \
} \
} \
template < const int EVFLAG, const int EFLAG,
const int NEWTON_PAIR, const int COUL_TYPE >
void PairCMMCommon::eval_verlet()
{
double ** const x = atom->x;
double ** const f = atom->f;
const double * const q = atom->q;
const int * const type = atom->type;
const int nlocal = atom->nlocal;
const int nall = nlocal + atom->nghost;
const double * const special_lj = force->special_lj;
const double * const special_coul = force->special_coul;
const double qqrd2e = force->qqrd2e;
const int inum = list->inum;
const int * const ilist = list->ilist;
const int * const numneigh = list->numneigh;
int * const * const firstneigh = list->firstneigh;
// loop over neighbors of my atoms
int ii,jj;
for (ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
double qtmp = (COUL_TYPE != CG_COUL_NONE) ? q[i] : 0.0;
const int itype = type[i];
const int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
int j2 = jlist[jj];
double factor_lj = 1.0;
double factor_coul = 1.0;
if (j2 >= nall) {
factor_lj = special_lj[j2/nall];
if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall];
j2 %= nall;
}
const int j = j2;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
const int jtype = type[j];
double evdwl = 0.0;
double ecoul = 0.0;
double fpair = 0.0;
const double r2inv = 1.0/rsq;
const int cgt=cg_type[itype][jtype];
if (rsq < cutsq[itype][jtype]) {
if (COUL_TYPE == CG_COUL_NONE) {
CG_LJ_INNER(EFLAG,fpair);
fpair *= r2inv;
} else {
double forcelj = 0.0;
double forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
CG_LJ_INNER(EFLAG,forcelj);
}
// coulomb with cutoff and screening
if ((COUL_TYPE == CG_COUL_CUT) || (COUL_TYPE == CG_COUL_DEBYE)) {
if (rsq < cut_coulsq[itype][jtype]) {
double r=sqrt(rsq);
double screen=exp(-kappa*r);
forcecoul = factor_coul * qqrd2e
* qtmp * q[j] * screen * (kappa + 1.0/r);
if (EFLAG) ecoul=factor_coul*qqrd2e
* qtmp*q[j] * screen / r;
}
}
if (COUL_TYPE == CG_COUL_LONG) {
if (rsq < cut_coulsq_global) {
const float rsqf = rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
const float r = sqrtf(rsq);
const double grij = g_ewald * r;
const double expm2 = exp(-grij*grij);
const double t = 1.0 / (1.0 + EWALD_P*grij);
const double erfc = t * (EWALD_A1+t*(EWALD_A2+t*(EWALD_A3+t*(EWALD_A4+t*EWALD_A5)))) * expm2;
const double prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (EFLAG) ecoul = prefactor*erfc;
if (factor_coul < 1.0) {
forcecoul -= (1.0-factor_coul)*prefactor;
if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
}
} else {
const int *int_rsq = (int *) &rsqf;
int itable = *int_rsq & ncoulmask;
itable >>= ncoulshiftbits;
const double fraction = (rsq - rtable[itable]) * drtable[itable];
const double table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (EFLAG) {
const double table2 = etable[itable] + fraction*detable[itable];
ecoul = qtmp*q[j] * table2;
}
if (factor_coul < 1.0) {
const double table2 = ctable[itable] + fraction*dctable[itable];
const double prefactor = qtmp*q[j] * table2;
forcecoul -= (1.0-factor_coul)*prefactor;
if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
}
}
}
}
fpair = (forcecoul + forcelj) * r2inv;
}
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,ecoul,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_compute();
}
/* ---------------------------------------------------------------------- */
template < const int NEWTON_PAIR, const int COUL_TYPE >
void PairCMMCommon::eval_inner()
{
double ** const x = atom->x;
double ** const f = atom->f;
const double * const q = atom->q;
const int * const type = atom->type;
const int nlocal = atom->nlocal;
const int nall = nlocal + atom->nghost;
const double * const special_lj = force->special_lj;
const double * const special_coul = force->special_coul;
const double qqrd2e = force->qqrd2e;
const int inum = listinner->inum;
const int * const ilist = listinner->ilist;
const int * const numneigh = listinner->numneigh;
int * const * const firstneigh = listinner->firstneigh;
const double cut_out_on = cut_respa[0];
const double cut_out_off = cut_respa[1];
const double cut_out_diff = cut_out_off - cut_out_on;
const double cut_out_on_sq = cut_out_on*cut_out_on;
const double cut_out_off_sq = cut_out_off*cut_out_off;
// loop over neighbors of my atoms
int ii,jj;
for (ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
double qtmp = (COUL_TYPE != CG_COUL_NONE) ? q[i] : 0.0;
const int itype = type[i];
const int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
int j2 = jlist[jj];
double factor_lj = 1.0;
double factor_coul = 1.0;
if (j2 >= nall) {
factor_lj = special_lj[j2/nall];
if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall];
j2 %= nall;
}
const int j = j2;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
const int jtype = type[j];
double evdwl = 0.0;
double ecoul = 0.0;
double fpair = 0.0;
const double r2inv = 1.0/rsq;
const int cgt=cg_type[itype][jtype];
if (rsq < cut_out_off_sq) {
if (COUL_TYPE == CG_COUL_NONE) {
CG_LJ_INNER(0,fpair);
fpair *= r2inv;
if (rsq > cut_out_on_sq) {
const double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw);
}
} else {
double forcelj = 0.0;
double forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
CG_LJ_INNER(0,forcelj);
}
forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
if (factor_coul < 1.0) forcecoul -= (1.0 -factor_coul)*forcecoul;
fpair = (forcecoul + forcelj) * r2inv;
if (rsq > cut_out_on_sq) {
const double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw);
}
}
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;
}
}
}
}
}
/* ---------------------------------------------------------------------- */
template < const int NEWTON_PAIR, const int COUL_TYPE >
void PairCMMCommon::eval_middle()
{
double ** const x = atom->x;
double ** const f = atom->f;
const double * const q = atom->q;
const int * const type = atom->type;
const int nlocal = atom->nlocal;
const int nall = nlocal + atom->nghost;
const double * const special_lj = force->special_lj;
const double * const special_coul = force->special_coul;
const double qqrd2e = force->qqrd2e;
const int inum = listmiddle->inum;
const int * const ilist = listmiddle->ilist;
const int * const numneigh = listmiddle->numneigh;
int * const * const firstneigh = listmiddle->firstneigh;
const double cut_in_off = cut_respa[0];
const double cut_in_on = cut_respa[1];
const double cut_out_on = cut_respa[2];
const double cut_out_off = cut_respa[3];
const double cut_in_diff = cut_in_on - cut_in_off;
const double cut_out_diff = cut_out_off - cut_out_on;
const double cut_in_off_sq = cut_in_off*cut_in_off;
const double cut_in_on_sq = cut_in_on*cut_in_on;
const double cut_out_on_sq = cut_out_on*cut_out_on;
const double cut_out_off_sq = cut_out_off*cut_out_off;
// loop over neighbors of my atoms
int ii,jj;
for (ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
double qtmp = (COUL_TYPE != CG_COUL_NONE) ? q[i] : 0.0;
const int itype = type[i];
const int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
int j2 = jlist[jj];
double factor_lj = 1.0;
double factor_coul = 1.0;
if (j2 >= nall) {
factor_lj = special_lj[j2/nall];
if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall];
j2 %= nall;
}
const int j = j2;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
const int jtype = type[j];
double evdwl = 0.0;
double ecoul = 0.0;
double fpair = 0.0;
const double r2inv = 1.0/rsq;
const int cgt=cg_type[itype][jtype];
if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
if (COUL_TYPE == CG_COUL_NONE) {
CG_LJ_INNER(0,fpair);
fpair *= r2inv;
if (rsq < cut_in_on_sq) {
const double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
}
if (rsq > cut_out_on_sq) {
const double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
}
} else {
double forcelj = 0.0;
double forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
CG_LJ_INNER(0,forcelj);
}
forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
if (factor_coul < 1.0) forcecoul -= (1.0 -factor_coul)*forcecoul;
fpair = (forcecoul + forcelj) * r2inv;
if (rsq < cut_in_on_sq) {
const double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
}
if (rsq > cut_out_on_sq) {
const double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
}
}
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;
}
}
}
}
}
/* ---------------------------------------------------------------------- */
template < const int EVFLAG, const int EFLAG, const int VFLAG,
const int NEWTON_PAIR, const int COUL_TYPE >
void PairCMMCommon::eval_outer()
{
double ** const x = atom->x;
double ** const f = atom->f;
const double * const q = atom->q;
const int * const type = atom->type;
const int nlocal = atom->nlocal;
const int nall = nlocal + atom->nghost;
const double * const special_lj = force->special_lj;
const double * const special_coul = force->special_coul;
const double qqrd2e = force->qqrd2e;
const int inum = listouter->inum;
const int * const ilist = listouter->ilist;
const int * const numneigh = listouter->numneigh;
int * const * const firstneigh = listouter->firstneigh;
const double cut_in_off = cut_respa[2];
const double cut_in_on = cut_respa[3];
const double cut_in_diff = cut_in_on - cut_in_off;
const double cut_in_off_sq = cut_in_off*cut_in_off;
const double cut_in_on_sq = cut_in_on*cut_in_on;
// loop over neighbors of my atoms
int ii,jj;
for (ii = 0; ii < inum; ii++) {
const int i = ilist[ii];
const double xtmp = x[i][0];
const double ytmp = x[i][1];
const double ztmp = x[i][2];
double qtmp = (COUL_TYPE != CG_COUL_NONE) ? q[i] : 0.0;
const int itype = type[i];
const int * const jlist = firstneigh[i];
const int jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
int j2 = jlist[jj];
double factor_lj = 1.0;
double factor_coul = 1.0;
if (j2 >= nall) {
factor_lj = special_lj[j2/nall];
if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall];
j2 %= nall;
}
const int j = j2;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
const int jtype = type[j];
double evdwl = 0.0;
double ecoul = 0.0;
double fpair = 0.0;
const double r2inv = 1.0/rsq;
const int cgt=cg_type[itype][jtype];
if (rsq < cutsq[itype][jtype]) {
if (COUL_TYPE == CG_COUL_NONE) {
double forcelj=0.0;
if (rsq > cut_in_off_sq) {
CG_LJ_INNER(0,forcelj);
fpair = forcelj*r2inv;
if (rsq < cut_in_on_sq) {
const double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
}
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;
}
}
CG_LJ_ENERGY(EFLAG);
if (VFLAG) {
if (rsq <= cut_in_off_sq) {
CG_LJ_INNER(0,fpair);
fpair *= r2inv;
} else if (rsq < cut_in_on_sq) {
fpair = forcelj*r2inv;
}
}
if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR,
evdwl,ecoul,fpair,delx,dely,delz);
} else {
double forcelj = 0.0;
double forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
CG_LJ_INNER(EFLAG,forcelj);
}
// coulomb with cutoff and screening
if ((COUL_TYPE == CG_COUL_CUT) || (COUL_TYPE == CG_COUL_DEBYE)) {
if (rsq < cut_coulsq[itype][jtype]) {
double r=sqrt(rsq);
double screen=exp(-kappa*r);
forcecoul = factor_coul * qqrd2e
* qtmp * q[j] * screen * (kappa + 1.0/r);
if (EFLAG) ecoul=factor_coul*qqrd2e
* qtmp*q[j] * screen / r;
}
}
if (COUL_TYPE == CG_COUL_LONG) {
if (rsq < cut_coulsq_global) {
const float rsqf = rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
const float r = sqrtf(rsq);
const double grij = g_ewald * r;
const double expm2 = exp(-grij*grij);
const double t = 1.0 / (1.0 + EWALD_P*grij);
const double erfc = t * (EWALD_A1+t*(EWALD_A2+t*(EWALD_A3+t*(EWALD_A4+t*EWALD_A5)))) * expm2;
const double prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (EFLAG) ecoul = prefactor*erfc;
if (factor_coul < 1.0) {
forcecoul -= (1.0-factor_coul)*prefactor;
if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
}
} else {
const int *int_rsq = (int *) &rsqf;
int itable = *int_rsq & ncoulmask;
itable >>= ncoulshiftbits;
const double fraction = (rsq - rtable[itable]) * drtable[itable];
const double table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (EFLAG) {
const double table2 = etable[itable] + fraction*detable[itable];
ecoul = qtmp*q[j] * table2;
}
if (factor_coul < 1.0) {
const double table2 = ctable[itable] + fraction*dctable[itable];
const double prefactor = qtmp*q[j] * table2;
forcecoul -= (1.0-factor_coul)*prefactor;
if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
}
}
}
}
fpair = (forcecoul + forcelj) * r2inv;
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,ecoul,fpair,delx,dely,delz);
}
}
}
}
}
/* ------------------------------------------------------------------------ */
}
#undef EWALD_F
#undef EWALD_P
#undef EWALD_A1
#undef EWALD_A2
#undef EWALD_A3
#undef EWALD_A4
#undef EWALD_A5
#endif

36
src/USER-CG-CMM/style_user_cg_cmm.h Normal file
View File

@ -0,0 +1,36 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
// add new include files in appropriate Include ifdef
// add new style keywords and class names in appropriate Class ifdef
// see style.h for examples
#ifdef AngleInclude
#include "angle_cg_cmm.h"
#endif
#ifdef AngleClass
AngleStyle(cg/cmm,AngleCGCMM)
#endif
#ifdef PairInclude
#include "pair_cg_cmm.h"
#include "pair_cg_cmm_coul_cut.h"
#include "pair_cg_cmm_coul_long.h"
#endif
#ifdef PairClass
PairStyle(cg/cmm,PairCGCMM)
PairStyle(cg/cmm/coul/cut,PairCGCMMCoulCut)
PairStyle(cg/cmm/coul/long,PairCGCMMCoulLong)
#endif

1
src/style_user_packages.h
View File

@ -15,4 +15,5 @@
// see the README files in individual user-package directories for details
#include "style_user_ackland.h"
#include "style_user_cg_cmm.h"
#include "style_user_ewaldn.h"