use CHARMM energy conversion factor with new CHARMM pair styles

src/KSPACE/pair_lj_charmmfsw_coul_long.cpp

@@ -25,6 +25,7 @@
 #include <string.h>
 #include "pair_lj_charmmfsw_coul_long.h"
 #include "atom.h"
+#include "update.h"
 #include "comm.h"
 #include "force.h"
 #include "kspace.h"
@@ -61,6 +62,11 @@ PairLJCharmmfswCoulLong::PairLJCharmmfswCoulLong(LAMMPS *lmp) : Pair(lmp)
   // short-range/long-range flag accessed by DihedralCharmmfsw
 
   dihedflag = 1;
+
+  // switch qqr2e from LAMMPS value to CHARMM value
+
+  if (strcmp(update->unit_style,"real") == 0)
+    force->qqr2e = force->qqr2e_charmm_real;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -87,6 +93,11 @@ PairLJCharmmfswCoulLong::~PairLJCharmmfswCoulLong()
     }
     if (ftable) free_tables();
   }
+
+  // switch qqr2e back from CHARMM value to LAMMPS value
+
+  if (strcmp(update->unit_style,"real") == 0)
+    force->qqr2e = force->qqr2e_lammps_real;
 }
 
 /* ---------------------------------------------------------------------- */

src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp

@@ -25,6 +25,7 @@
 #include <string.h>
 #include "pair_lj_charmmfsw_coul_charmmfsh.h"
 #include "atom.h"
+#include "update.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
@@ -46,6 +47,11 @@ PairLJCharmmfswCoulCharmmfsh::PairLJCharmmfswCoulCharmmfsh(LAMMPS *lmp) :
   // short-range/long-range flag accessed by DihedralCharmmfsw
 
   dihedflag = 0;
+
+  // switch qqr2e from LAMMPS value to CHARMM value
+
+  if (strcmp(update->unit_style,"real") == 0)
+    force->qqr2e = force->qqr2e_charmm_real;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -71,6 +77,11 @@ PairLJCharmmfswCoulCharmmfsh::~PairLJCharmmfswCoulCharmmfsh()
       memory->destroy(lj14_4);
     }
   }
+
+  // switch qqr2e back from CHARMM value to LAMMPS value
+
+  if (strcmp(update->unit_style,"real") == 0)
+    force->qqr2e = force->qqr2e_lammps_real;
 }
 
 /* ---------------------------------------------------------------------- */

src/force.cpp

@@ -53,6 +53,8 @@ Force::Force(LAMMPS *lmp) : Pointers(lmp)
   special_extra = 0;
 
   dielectric = 1.0;
+  qqr2e_lammps_real = 332.06371;          // these constants are toggled
+  qqr2e_charmm_real = 332.0716;           // by new CHARMM pair styles
 
   pair = NULL;
   bond = NULL;

src/force.h

@@ -43,6 +43,9 @@ class Force : protected Pointers {
   double femtosecond;                // 1 femtosecond in native units
   double qelectron;                  // 1 electron charge abs() in native units
 
+  double qqr2e_lammps_real;          // different versions of this constant
+  double qqr2e_charmm_real;          // used by new CHARMM pair styles
+
   int newton,newton_pair,newton_bond;   // Newton's 3rd law settings
 
   class Pair *pair;

src/update.cpp

@@ -154,7 +154,7 @@ void Update::set_units(const char *style)
     force->ftm2v = 1.0 / 48.88821291 / 48.88821291;
     force->mv2d = 1.0 / 0.602214129;
     force->nktv2p = 68568.415;
-    force->qqr2e = 332.06371;
+    force->qqr2e = 332.06371;     // see also force->qqr2d_lammps_real
     force->qe2f = 23.060549;
     force->vxmu2f = 1.4393264316e4;
     force->xxt2kmu = 0.1;