diff --git a/src/USER-CG-CMM/pair_cmm_common.cpp b/src/USER-CG-CMM/pair_cmm_common.cpp index 97d8323c1d..8a09a16295 100644 --- a/src/USER-CG-CMM/pair_cmm_common.cpp +++ b/src/USER-CG-CMM/pair_cmm_common.cpp @@ -41,8 +41,7 @@ PairCMMCommon::PairCMMCommon(class LAMMPS *lmp) : Pair(lmp) * clean up common arrays * * ---------------------------------------------------------------------- */ -PairCMMCommon::~PairCMMCommon() -{ +PairCMMCommon::~PairCMMCommon() { if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_int_array(cg_type); @@ -447,16 +446,21 @@ double PairCMMCommon::eval_single(int coul_type, int i, int j, int itype, int jt 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_force = cgpref*eps * (cgpow1*pow(ratio,cgpow1) + - cgpow2*pow(ratio,cgpow2))/rsq; 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."); + error->all("single energy computation with long-range 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."); + const double r2inv = 1.0/rsq; + const double rinv = sqrt(r2inv); + const double qscreen=exp(-kappa*sqrt(rsq)); + coul_force = force->qqrd2e * atom->q[i]*atom->q[j]*rinv * qscreen * (kappa + rinv); + coul_erg = force->qqrd2e * atom->q[i]*atom->q[j]*rinv * qscreen; + // error->all("single energy computation with coulomb not supported by CG potentials."); } else if (coul_type == CG_COUL_NONE) { ; // do nothing } else { diff --git a/src/USER-CG-CMM/pair_cmm_common.h b/src/USER-CG-CMM/pair_cmm_common.h index 1635ba9cfd..46b3111579 100644 --- a/src/USER-CG-CMM/pair_cmm_common.h +++ b/src/USER-CG-CMM/pair_cmm_common.h @@ -10,7 +10,7 @@ /* ---------------------------------------------------------------------- Common functionality for the CMM coarse grained MD potentials. - Contributing author: Axel Kohlmeyer + Contributing author: Axel Kohlmeyer ------------------------------------------------------------------------- */ #ifndef PAIR_CMM_COMMON_H @@ -247,11 +247,11 @@ namespace LAMMPS_NS { 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); + double qscreen=exp(-kappa*r); forcecoul = factor_coul * qqrd2e - * qtmp * q[j] * screen * (kappa + 1.0/r); + * qtmp * q[j] * qscreen * (kappa + 1.0/r); if (EFLAG) ecoul=factor_coul*qqrd2e - * qtmp*q[j] * screen / r; + * qtmp*q[j] * qscreen / r; } } @@ -644,11 +644,11 @@ namespace LAMMPS_NS { 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); + double qscreen=exp(-kappa*r); forcecoul = factor_coul * qqrd2e - * qtmp * q[j] * screen * (kappa + 1.0/r); + * qtmp * q[j] * qscreen * (kappa + 1.0/r); if (EFLAG) ecoul=factor_coul*qqrd2e - * qtmp*q[j] * screen / r; + * qtmp*q[j] * qscreen / r; } } diff --git a/src/finish.cpp b/src/finish.cpp index d9708edb5b..72acf8497c 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -82,7 +82,6 @@ void Finish::end(int flag) // overall loop time // use actual natoms, in case atoms were lost - natoms; double rlocal = atom->nlocal; MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/pair.cpp b/src/pair.cpp index 1453cd701f..73daa5eeff 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -390,6 +390,44 @@ void Pair::ev_tally(int i, int j, int nlocal, int newton_pair, } } +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + can use this version with full neighbor lists +------------------------------------------------------------------------- */ + +void Pair::ev_tally_full(int i, double evdwl, double ecoul, double fpair, + double delx, double dely, double delz) +{ + + if (eflag_either) { + if (eflag_global) { + eng_vdwl += 0.5*evdwl; + eng_coul += 0.5*ecoul; + } + if (eflag_atom) eatom[i] += 0.5 * (evdwl + ecoul); + } + + if (vflag_either) { + if (vflag_global) { + virial[0] += 0.5*delx*delx*fpair; + virial[1] += 0.5*dely*dely*fpair; + virial[2] += 0.5*delz*delz*fpair; + virial[3] += 0.5*delx*dely*fpair; + virial[4] += 0.5*delx*delz*fpair; + virial[5] += 0.5*dely*delz*fpair; + } + } + + if (vflag_atom) { + vatom[i][0] += 0.5*delx*delx*fpair; + vatom[i][1] += 0.5*dely*dely*fpair; + vatom[i][2] += 0.5*delz*delz*fpair; + vatom[i][3] += 0.5*delx*dely*fpair; + vatom[i][4] += 0.5*delx*delz*fpair; + vatom[i][5] += 0.5*dely*delz*fpair; + } +} + /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz diff --git a/src/pair.h b/src/pair.h index e896d98abf..653e57119f 100644 --- a/src/pair.h +++ b/src/pair.h @@ -121,6 +121,7 @@ class Pair : protected Pointers { void ev_setup(int, int); void ev_tally(int, int, int, int, double, double, double, double, double, double); + void ev_tally_full(int, double, double, double, double, double, double); void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double, double); void ev_tally3(int, int, int, double, double,