git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3295 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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) 
+    lj_force = cgpref*eps * (cgpow1*pow(ratio,cgpow1) 
-                                       - cgpow2*pow(ratio,cgpow2));
+                            - 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 {

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 <akohlmey@cmm.chem.upenn.edu>
+   Contributing author: Axel Kohlmeyer <akohlmey@gmail.com>
 ------------------------------------------------------------------------- */
 
 #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;
               }
             }
 

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);
     

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

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,