From 60c3f3d64c5872febe9def2f0b94651b5304e415 Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Thu, 22 Jun 2017 09:15:15 -0600
Subject: [PATCH] use CHARMM energy conversion factor with new CHARMM pair
styles
---
src/KSPACE/pair_lj_charmmfsw_coul_long.cpp | 11 +++++++++++
src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp | 11 +++++++++++
src/force.cpp | 2 ++
src/force.h | 3 +++
src/update.cpp | 2 +-
5 files changed, 28 insertions(+), 1 deletion(-)
diff --git a/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp b/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
index 6e17a9bbd7..9ba59b6995 100644
--- a/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
+++ b/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;
}
/* ---------------------------------------------------------------------- */
diff --git a/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp b/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp
index 1e34b06478..f1366a379d 100644
--- a/src/MOLECULE/pair_lj_charmmfsw_coul_charmmfsh.cpp
+++ b/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;
}
/* ---------------------------------------------------------------------- */
diff --git a/src/force.cpp b/src/force.cpp
index 3dd28cc710..33e6630406 100644
--- a/src/force.cpp
+++ b/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;
diff --git a/src/force.h b/src/force.h
index f2d9abc7dd..edaac1b527 100644
--- a/src/force.h
+++ b/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;
diff --git a/src/update.cpp b/src/update.cpp
index 5599dc6c88..e4c85dde73 100644
--- a/src/update.cpp
+++ b/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;
--
GitLab