From 96dae5a76c03ff330c366f2d6fd280ba98033355 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 15 May 2008 23:03:19 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1826 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-CG-CMM/angle_cg_cmm.cpp | 91 +++++++++++++++++++++++++++-- src/USER-CG-CMM/angle_cg_cmm.h | 4 ++ src/USER-CG-CMM/pair_cg_cmm.h | 2 +- src/USER-CG-CMM/pair_cmm_common.cpp | 5 +- src/USER-CG-CMM/pair_cmm_common.h | 1 + 5 files changed, 97 insertions(+), 6 deletions(-) diff --git a/src/USER-CG-CMM/angle_cg_cmm.cpp b/src/USER-CG-CMM/angle_cg_cmm.cpp index 9506813c4d..63129ce5a0 100644 --- a/src/USER-CG-CMM/angle_cg_cmm.cpp +++ b/src/USER-CG-CMM/angle_cg_cmm.cpp @@ -49,6 +49,88 @@ AngleCGCMM::~AngleCGCMM() /* ---------------------------------------------------------------------- */ +void AngleCGCMM::ev_tally_lj13(int i, int j, int nlocal, int newton_bond, + double evdwl, double fpair, + double delx, double dely, double delz) +{ + double v[6]; + + if (eflag_either) { + if (eflag_global) { + if (newton_bond) { + energy += evdwl; + } else { + if (i < nlocal) + energy += 0.5*evdwl; + if (j < nlocal) + energy += 0.5*evdwl; + } + } + if (eflag_atom) { + if (newton_bond || i < nlocal) eatom[i] += 0.5*evdwl; + if (newton_bond || j < nlocal) eatom[i] += 0.5*evdwl; + } + } + + if (vflag_either) { + v[0] = delx*delx*fpair; + v[1] = dely*dely*fpair; + v[2] = delz*delz*fpair; + v[3] = delx*dely*fpair; + v[4] = delx*delz*fpair; + v[5] = dely*delz*fpair; + + if (vflag_global) { + if (newton_bond) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } else { + if (i < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + if (j < nlocal) { + virial[0] += 0.5*v[0]; + virial[1] += 0.5*v[1]; + virial[2] += 0.5*v[2]; + virial[3] += 0.5*v[3]; + virial[4] += 0.5*v[4]; + virial[5] += 0.5*v[5]; + } + } + } + + if (vflag_atom) { + if (newton_bond || i < nlocal) { + vatom[i][0] += 0.5*v[0]; + vatom[i][1] += 0.5*v[1]; + vatom[i][2] += 0.5*v[2]; + vatom[i][3] += 0.5*v[3]; + vatom[i][4] += 0.5*v[4]; + vatom[i][5] += 0.5*v[5]; + } + if (newton_bond || j < nlocal) { + vatom[j][0] += 0.5*v[0]; + vatom[j][1] += 0.5*v[1]; + vatom[j][2] += 0.5*v[2]; + vatom[j][3] += 0.5*v[3]; + vatom[j][4] += 0.5*v[4]; + vatom[j][5] += 0.5*v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + void AngleCGCMM::compute(int eflag, int vflag) { int i1,i2,i3,n,type; @@ -142,7 +224,7 @@ void AngleCGCMM::compute(int eflag, int vflag) dtheta = acos(c) - theta0[type]; tk = k[type] * dtheta; - if (eflag) eangle = tk*dtheta + e13; + if (eflag) eangle = tk*dtheta; a = -2.0 * tk * s; a11 = a*c / rsq1; @@ -176,11 +258,12 @@ void AngleCGCMM::compute(int eflag, int vflag) 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); + + if (evflag) ev_tally_lj13(i1,i3,nlocal,newton_bond, + e13,f13,delx3,dely3,delz3); + } } diff --git a/src/USER-CG-CMM/angle_cg_cmm.h b/src/USER-CG-CMM/angle_cg_cmm.h index d60dc3243d..c34ac71f11 100644 --- a/src/USER-CG-CMM/angle_cg_cmm.h +++ b/src/USER-CG-CMM/angle_cg_cmm.h @@ -34,6 +34,10 @@ class AngleCGCMM : public Angle, public CGCMMParms { void read_restart(FILE *); double single(int, int, int, int); + protected: + void ev_tally_lj13(int, int, int, int, double, double, + double, double, double); + private: double *k,*theta0; int *cg_type; diff --git a/src/USER-CG-CMM/pair_cg_cmm.h b/src/USER-CG-CMM/pair_cg_cmm.h index 3ec2c0a96a..f3de701298 100644 --- a/src/USER-CG-CMM/pair_cg_cmm.h +++ b/src/USER-CG-CMM/pair_cg_cmm.h @@ -33,7 +33,7 @@ namespace LAMMPS_NS { double single(int, int, int, int, double, double, double, double &); private: - // disable default destructor + // disable default constructor PairCGCMM(); }; } diff --git a/src/USER-CG-CMM/pair_cmm_common.cpp b/src/USER-CG-CMM/pair_cmm_common.cpp index f91910c01b..ba9b55f7ed 100644 --- a/src/USER-CG-CMM/pair_cmm_common.cpp +++ b/src/USER-CG-CMM/pair_cmm_common.cpp @@ -259,6 +259,7 @@ double PairCMMCommon::init_one(int i, int j) if (allocated_coul) { mycut = MAX(cut_lj[i][j],cut_coul[i][j]); + cut[i][j] = mycut; 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) { @@ -274,7 +275,8 @@ double PairCMMCommon::init_one(int i, int j) lj4[j][i] = lj4[i][j]; offset[j][i] = offset[i][j]; cg_type[j][i] = cg_type[i][j]; - + cut[j][i] = mycut; + if (allocated_coul) { cut_lj[j][i]=cut_lj[i][j]; cut_ljsq[j][i]=cut_ljsq[i][j]; @@ -407,6 +409,7 @@ void PairCMMCommon::read_restart_settings(FILE *fp) 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); + cut_coulsq_global = cut_coul_global*cut_coul_global; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-CG-CMM/pair_cmm_common.h b/src/USER-CG-CMM/pair_cmm_common.h index f91bb34ad1..61af5ebfc4 100644 --- a/src/USER-CG-CMM/pair_cmm_common.h +++ b/src/USER-CG-CMM/pair_cmm_common.h @@ -201,6 +201,7 @@ namespace LAMMPS_NS { 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];