diff --git a/src/USER-CG-CMM/Install.sh b/src/USER-CG-CMM/Install.sh
index 1562a9896..ac4705f01 100644
--- a/src/USER-CG-CMM/Install.sh
+++ b/src/USER-CG-CMM/Install.sh
@@ -1,62 +1,69 @@
# Install/unInstall package files in LAMMPS
# do not install child files if parent does not exist
if (test $1 = 1) then
if (test -e ../angle_harmonic.cpp) then
cp angle_cg_cmm.h ..
cp angle_cg_cmm.cpp ..
cp angle_sdk.h ..
cp angle_sdk.cpp ..
fi
if (test -e ../pppm.cpp) then
cp pair_lj_sdk_coul_long.cpp ..
cp pair_cg_cmm_coul_long.cpp ..
cp pair_lj_sdk_coul_long.h ..
cp pair_cg_cmm_coul_long.h ..
fi
+ if (test -e ../msm.cpp) then
+ cp pair_lj_sdk_coul_msm.cpp ..
+ cp pair_lj_sdk_coul_msm.h ..
+ fi
+
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_lj_sdk.cpp ..
cp pair_lj_sdk.h ..
cp lj_sdk_common.h ..
elif (test $1 = 0) then
rm -f ../angle_cg_cmm.h
rm -f ../angle_cg_cmm.cpp
rm -f ../cg_cmm_parms.h
rm -f ../cg_cmm_parms.cpp
rm -f ../pair_cmm_common.h
rm -f ../pair_cmm_common.cpp
rm -f ../pair_cg_cmm.cpp
rm -f ../pair_cg_cmm.h
rm -f ../pair_cg_cmm_coul_cut.cpp
rm -f ../pair_cg_cmm_coul_cut.h
rm -f ../pair_cg_cmm_coul_long.cpp
rm -f ../pair_cg_cmm_coul_long.h
rm -f ../lj_sdk_common.h
rm -f ../angle_sdk.h
rm -f ../angle_sdk.cpp
rm -f ../pair_lj_sdk.cpp
rm -f ../pair_lj_sdk.h
rm -f ../pair_lj_sdk_coul_long.cpp
rm -f ../pair_lj_sdk_coul_long.h
+ rm -f ../pair_lj_sdk_coul_msm.cpp
+ rm -f ../pair_lj_sdk_coul_msm.h
fi
diff --git a/src/USER-CG-CMM/Package.sh b/src/USER-CG-CMM/Package.sh
index 317648236..6f481d663 100644
--- a/src/USER-CG-CMM/Package.sh
+++ b/src/USER-CG-CMM/Package.sh
@@ -1,38 +1,44 @@
# Update package files in LAMMPS
# cp package file to src if doesn't exist or is different
# do not copy molecular and kspace files if corresponding versions do not exist
for file in *.cpp *.h; do
if (test $file = angle_cg_cmm.cpp -a ! -e ../angle_harmonic.cpp) then
continue
fi
if (test $file = angle_cg_cmm.h -a ! -e ../angle_harmonic.h) then
continue
fi
if (test $file = pair_cg_cmm_coul_long.cpp -a ! -e ../pair_lj_cut_coul_long.cpp) then
continue
fi
if (test $file = pair_cg_cmm_coul_long.h -a ! -e ../pair_lj_cut_coul_long.h) then
continue
fi
if (test $file = angle_sdk.cpp -a ! -e ../angle_harmonic.cpp) then
continue
fi
if (test $file = angle_sdk.h -a ! -e ../angle_harmonic.h) then
continue
fi
if (test $file = pair_lj_sdk_coul_long.cpp -a ! -e ../pair_lj_cut_coul_long.cpp) then
continue
fi
if (test $file = pair_lj_sdk_coul_long.h -a ! -e ../pair_lj_cut_coul_long.h) then
continue
fi
+ if (test $file = pair_lj_sdk_coul_msm.cpp -a ! -e ../pair_lj_cut_coul_msm.cpp) then
+ continue
+ fi
+ if (test $file = pair_lj_sdk_coul_msm.h -a ! -e ../pair_lj_cut_coul_msm.h) then
+ continue
+ fi
if (test ! -e ../$file) then
echo " creating src/$file"
cp $file ..
elif (test "`diff --brief $file ../$file`" != "") then
echo " updating src/$file"
cp $file ..
fi
done
diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_long.h b/src/USER-CG-CMM/pair_lj_sdk_coul_long.h
index c929712ec..8000507e2 100644
--- a/src/USER-CG-CMM/pair_lj_sdk_coul_long.h
+++ b/src/USER-CG-CMM/pair_lj_sdk_coul_long.h
@@ -1,82 +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(const char *, int &);
+ virtual void *extract(const 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
diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_msm.cpp b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.cpp
new file mode 100644
index 000000000..fe8e8ce9b
--- /dev/null
+++ b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.cpp
@@ -0,0 +1,307 @@
+/* ----------------------------------------------------------------------
+ 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_msm.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;
+
+/* ---------------------------------------------------------------------- */
+
+PairLJSDKCoulMSM::PairLJSDKCoulMSM(LAMMPS *lmp) : PairLJSDKCoulLong(lmp)
+{
+ ewaldflag = pppmflag = 0;
+ msmflag = 1;
+ respa_enable = 0;
+ ftable = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJSDKCoulMSM::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_msm<1,1,1>();
+ else eval_msm<1,1,0>();
+ } else {
+ if (force->newton_pair) eval_msm<1,0,1>();
+ else eval_msm<1,0,0>();
+ }
+ } else {
+ if (force->newton_pair) eval_msm<0,0,1>();
+ else eval_msm<0,0,0>();
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
+void PairLJSDKCoulMSM::eval_msm()
+{
+ 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 egamma,fgamma,prefactor;
+
+ 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++) {
+ forcecoul = forcelj = evdwl = ecoul = 0.0;
+
+ 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);
+ prefactor = qqrd2e * qtmp*q[j]/r;
+ fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
+ forcecoul = prefactor * fgamma;
+ if (EFLAG) {
+ egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
+ ecoul = prefactor*egamma;
+ }
+ if (factor_coul < 1.0) {
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ if (EFLAG) ecoul -= (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 (EFLAG) ecoul = qtmp*q[j] *
+ (etable[itable] + fraction*detable[itable]);
+ if (factor_coul < 1.0) {
+ table = ctable[itable] + fraction*dctable[itable];
+ prefactor = qtmp*q[j] * table;
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ } else forcecoul = 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];
+ }
+ forcelj *= factor_lj;
+ if (EFLAG) evdwl *= factor_lj;
+ }
+
+ fpair = (forcecoul + 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 (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;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJSDKCoulMSM::single(int i, int j, int itype, int jtype,
+ double rsq,
+ double factor_coul, double factor_lj,
+ double &fforce)
+{
+ double r2inv,r,fgamma,egamma,prefactor;
+ double fraction,table,forcecoul,forcelj,phicoul,philj;
+ int itable;
+
+ forcecoul = forcelj = phicoul = philj = 0.0;
+
+ r2inv = 1.0/rsq;
+ if (rsq < cut_coulsq) {
+ if (!ncoultablebits || rsq <= tabinnersq) {
+ r = sqrt(rsq);
+ prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
+ egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
+ fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
+ forcecoul = prefactor * fgamma;
+ phicoul = prefactor * egamma;
+ if (factor_coul < 1.0) {
+ forcecoul -= (1.0-factor_coul)*prefactor;
+ phicoul -= (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;
+ table = etable[itable] + fraction*detable[itable];
+ phicoul = 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;
+ phicoul -= (1.0-factor_coul)*prefactor;
+ }
+ }
+ }
+
+ 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];
+
+ forcelj = 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]);
+ }
+
+ fforce = (forcecoul + forcelj) * r2inv;
+
+ return phicoul + philj;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairLJSDKCoulMSM::extract(const 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;
+ if (strcmp(str,"cut_msm") == 0) return (void *) &cut_coul;
+ return NULL;
+}
+
diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_msm.h b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.h
new file mode 100644
index 000000000..af79b6e0b
--- /dev/null
+++ b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.h
@@ -0,0 +1,48 @@
+/* -*- 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/msm,PairLJSDKCoulMSM)
+PairStyle(cg/cmm/coul/msm,PairLJSDKCoulMSM)
+
+#else
+
+#ifndef LMP_PAIR_LJ_SDK_COUL_MSM_H
+#define LMP_PAIR_LJ_SDK_COUL_MSM_H
+
+#include "pair_lj_sdk_coul_long.h"
+
+namespace LAMMPS_NS {
+
+class PairLJSDKCoulMSM : public PairLJSDKCoulLong {
+ public:
+ PairLJSDKCoulMSM(class LAMMPS *);
+ virtual ~PairLJSDKCoulMSM() {};
+ virtual void compute(int, int);
+ virtual double single(int, int, int, int, double, double, double, double &);
+ virtual void *extract(const char *, int &);
+
+private:
+ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval_msm();
+
+};
+
+}
+
+#endif
+#endif