ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/branches/development/src/nonbonded/Electrostatic.cpp

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1535 by gezelter, Fri Dec 31 18:31:56 2010 UTC vs.
Revision 1877 by gezelter, Thu Jun 6 15:43:35 2013 UTC

#	Line 34 | Line 34
34		* work.  Good starting points are:
35		*
36		* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	<	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	<	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	<	* [4]  Vardeman & Gezelter, in progress (2009).
37	>	* [2]  Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).
38	>	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	>	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	>	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		#include <stdio.h>
#	Line 46 | Line 47
47		#include "nonbonded/Electrostatic.hpp"
48		#include "utils/simError.h"
49		#include "types/NonBondedInteractionType.hpp"
50	<	#include "types/DirectionalAtomType.hpp"
50	>	#include "types/FixedChargeAdapter.hpp"
51	>	#include "types/FluctuatingChargeAdapter.hpp"
52	>	#include "types/MultipoleAdapter.hpp"
53		#include "io/Globals.hpp"
54	+	#include "nonbonded/SlaterIntegrals.hpp"
55	+	#include "utils/PhysicalConstants.hpp"
56	+	#include "math/erfc.hpp"
57	+	#include "math/SquareMatrix.hpp"
58
59		namespace OpenMD {
60
61		Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
62	<	forceField_(NULL) {}
62	>	forceField_(NULL), info_(NULL),
63	>	haveCutoffRadius_(false),
64	>	haveDampingAlpha_(false),
65	>	haveDielectric_(false),
66	>	haveElectroSplines_(false)
67	>	{}
68
69		void Electrostatic::initialize() {
70	+
71	+	Globals* simParams_ = info_->getSimParams();
72
59	–	Globals* simParams_;
60	–
73		summationMap_["HARD"]               = esm_HARD;
74	+	summationMap_["NONE"]               = esm_HARD;
75		summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
76		summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
77		summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;
78	+	summationMap_["TAYLOR_SHIFTED"]     = esm_TAYLOR_SHIFTED;
79		summationMap_["REACTION_FIELD"]     = esm_REACTION_FIELD;
80		summationMap_["EWALD_FULL"]         = esm_EWALD_FULL;
81		summationMap_["EWALD_PME"]          = esm_EWALD_PME;
#	Line 76 | Line 90 | namespace OpenMD {
90		// Charge-Dipole, assuming charges are measured in electrons, and
91		// dipoles are measured in debyes
92		pre12_ = 69.13373;
93	<	// Dipole-Dipole, assuming dipoles are measured in debyes
93	>	// Dipole-Dipole, assuming dipoles are measured in Debye
94		pre22_ = 14.39325;
95		// Charge-Quadrupole, assuming charges are measured in electrons, and
96		// quadrupoles are measured in 10^-26 esu cm^2
97	<	// This unit is also known affectionately as an esu centi-barn.
97	>	// This unit is also known affectionately as an esu centibarn.
98		pre14_ = 69.13373;
99	<
99	>	// Dipole-Quadrupole, assuming dipoles are measured in debyes and
100	>	// quadrupoles in esu centibarns:
101	>	pre24_ = 14.39325;
102	>	// Quadrupole-Quadrupole, assuming esu centibarns:
103	>	pre44_ = 14.39325;
104	>
105		// conversions for the simulation box dipole moment
106		chargeToC_ = 1.60217733e-19;
107		angstromToM_ = 1.0e-10;
108		debyeToCm_ = 3.33564095198e-30;
109
110	<	// number of points for electrostatic splines
110	>	// Default number of points for electrostatic splines
111		np_ = 100;
112
113		// variables to handle different summation methods for long-range
#	Line 96 | Line 115 | namespace OpenMD {
115		summationMethod_ = esm_HARD;
116		screeningMethod_ = UNDAMPED;
117		dielectric_ = 1.0;
99	–	one_third_ = 1.0 / 3.0;
100	–	haveCutoffRadius_ = false;
101	–	haveDampingAlpha_ = false;
102	–	haveDielectric_ = false;
103	–	haveElectroSpline_ = false;
118
119		// check the summation method:
120		if (simParams_->haveElectrostaticSummationMethod()) {
#	Line 113 | Line 127 | namespace OpenMD {
127		} else {
128		// throw error
129		sprintf( painCave.errMsg,
130	<	"SimInfo error: Unknown electrostaticSummationMethod.\n"
130	>	"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
131		"\t(Input file specified %s .)\n"
132	<	"\telectrostaticSummationMethod must be one of: \"none\",\n"
133	<	"\t\"shifted_potential\", \"shifted_force\", or \n"
134	<	"\t\"reaction_field\".\n", myMethod.c_str() );
132	>	"\telectrostaticSummationMethod must be one of: \"hard\",\n"
133	>	"\t\"shifted_potential\", \"shifted_force\",\n"
134	>	"\t\"taylor_shifted\", or \"reaction_field\".\n",
135	>	myMethod.c_str() );
136		painCave.isFatal = 1;
137		simError();
138		}
#	Line 185 | Line 200 | namespace OpenMD {
200
201		// throw warning
202		sprintf( painCave.errMsg,
203	<	"Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
204	<	"\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
203	>	"Electrostatic::initialize: dampingAlpha was not specified in the\n"
204	>	"\tinput file.  A default value of %f (1/ang) will be used for the\n"
205	>	"\tcutoff of %f (ang).\n",
206		dampingAlpha_, cutoffRadius_);
207		painCave.severity = OPENMD_INFO;
208		painCave.isFatal = 0;
#	Line 207 | Line 223 | namespace OpenMD {
223
224		if (at->isElectrostatic())
225		addType(at);
226	<	}
226	>	}
227
212	–
213	–	cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
214	–	rcuti_ = 1.0 / cutoffRadius_;
215	–	rcuti2_ = rcuti_ * rcuti_;
216	–	rcuti3_ = rcuti2_ * rcuti_;
217	–	rcuti4_ = rcuti2_ * rcuti2_;
218	–
219	–	if (screeningMethod_ == DAMPED) {
220	–
221	–	alpha2_ = dampingAlpha_ * dampingAlpha_;
222	–	alpha4_ = alpha2_ * alpha2_;
223	–	alpha6_ = alpha4_ * alpha2_;
224	–	alpha8_ = alpha4_ * alpha4_;
225	–
226	–	constEXP_ = exp(-alpha2_ * cutoffRadius2_);
227	–	invRootPi_ = 0.56418958354775628695;
228	–	alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
229	–
230	–	c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
231	–	c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
232	–	c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
233	–	c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_;
234	–	c5c_ = 8.0 * alphaPi_ * alpha6_ + 7.0 * c4c_ * rcuti2_;
235	–	c6c_ = 16.0 * alphaPi_ * alpha8_ + 9.0 * c5c_ * rcuti2_;
236	–	} else {
237	–	c1c_ = rcuti_;
238	–	c2c_ = c1c_ * rcuti_;
239	–	c3c_ = 3.0 * c2c_ * rcuti_;
240	–	c4c_ = 5.0 * c3c_ * rcuti2_;
241	–	c5c_ = 7.0 * c4c_ * rcuti2_;
242	–	c6c_ = 9.0 * c5c_ * rcuti2_;
243	–	}
244	–
228		if (summationMethod_ == esm_REACTION_FIELD) {
229		preRF_ = (dielectric_ - 1.0) /
230	<	((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
248	<	preRF2_ = 2.0 * preRF_;
230	>	((2.0 * dielectric_ + 1.0) * pow(cutoffRadius_,3) );
231		}
232
233	<	RealType dx = cutoffRadius_ / RealType(np_ - 1);
234	<	RealType rval;
235	<	vector<RealType> rvals;
236	<	vector<RealType> yvals;
237	<	for (int i = 0; i < np_; i++) {
238	<	rval = RealType(i) * dx;
239	<	rvals.push_back(rval);
240	<	yvals.push_back(erfc(dampingAlpha_ * rval));
259	<	}
260	<	erfcSpline_ = new CubicSpline();
261	<	erfcSpline_->addPoints(rvals, yvals);
262	<	haveElectroSpline_ = true;
263	<
264	<	initialized_ = true;
265	<	}
266	<
267	<	void Electrostatic::addType(AtomType* atomType){
233	>	RealType b0c, b1c, b2c, b3c, b4c, b5c;
234	>	RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5;
235	>	RealType a2, expTerm, invArootPi;
236	>
237	>	RealType r = cutoffRadius_;
238	>	RealType r2 = r * r;
239	>	RealType ric = 1.0 / r;
240	>	RealType ric2 = ric * ric;
241
242	<	ElectrostaticAtomData electrostaticAtomData;
243	<	electrostaticAtomData.is_Charge = false;
244	<	electrostaticAtomData.is_Dipole = false;
245	<	electrostaticAtomData.is_SplitDipole = false;
246	<	electrostaticAtomData.is_Quadrupole = false;
242	>	if (screeningMethod_ == DAMPED) {
243	>	a2 = dampingAlpha_ * dampingAlpha_;
244	>	invArootPi = 1.0 / (dampingAlpha_ * sqrt(M_PI));
245	>	expTerm = exp(-a2 * r2);
246	>	// values of Smith's B_l functions at the cutoff radius:
247	>	b0c = erfc(dampingAlpha_ * r) / r;
248	>	b1c = (      b0c     + 2.0*a2     * expTerm * invArootPi) / r2;
249	>	b2c = (3.0 * b1c + pow(2.0*a2, 2) * expTerm * invArootPi) / r2;
250	>	b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
251	>	b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
252	>	b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
253	>	selfMult_ = b0c + a2 * invArootPi;
254	>	} else {
255	>	a2 = 0.0;
256	>	b0c = 1.0 / r;
257	>	b1c = (      b0c) / r2;
258	>	b2c = (3.0 * b1c) / r2;
259	>	b3c = (5.0 * b2c) / r2;
260	>	b4c = (7.0 * b3c) / r2;
261	>	b5c = (9.0 * b4c) / r2;
262	>	selfMult_ = b0c;
263	>	}
264
265	<	if (atomType->isCharge()) {
266	<	GenericData* data = atomType->getPropertyByName("Charge");
265	>	// higher derivatives of B_0 at the cutoff radius:
266	>	db0c_1 = -r * b1c;
267	>	db0c_2 =     -b1c + r2 * b2c;
268	>	db0c_3 =          3.0*r*b2c  - r2*r*b3c;
269	>	db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
270	>	db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
271	>
272
273	<	if (data == NULL) {
274	<	sprintf( painCave.errMsg, "Electrostatic::addType could not find "
275	<	"Charge\n"
276	<	"\tparameters for atomType %s.\n",
277	<	atomType->getName().c_str());
278	<	painCave.severity = OPENMD_ERROR;
279	<	painCave.isFatal = 1;
280	<	simError();
281	<	}
282	<
288	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
289	<	if (doubleData == NULL) {
290	<	sprintf( painCave.errMsg,
291	<	"Electrostatic::addType could not convert GenericData to "
292	<	"Charge for\n"
293	<	"\tatom type %s\n", atomType->getName().c_str());
294	<	painCave.severity = OPENMD_ERROR;
295	<	painCave.isFatal = 1;
296	<	simError();
297	<	}
298	<	electrostaticAtomData.is_Charge = true;
299	<	electrostaticAtomData.charge = doubleData->getData();
300	<	}
273	>	// working variables for the splines:
274	>	RealType ri, ri2;
275	>	RealType b0, b1, b2, b3, b4, b5;
276	>	RealType db0_1, db0_2, db0_3, db0_4, db0_5;
277	>	RealType f, fc, f0;
278	>	RealType g, gc, g0, g1, g2, g3, g4;
279	>	RealType h, hc, h1, h2, h3, h4;
280	>	RealType s, sc, s2, s3, s4;
281	>	RealType t, tc, t3, t4;
282	>	RealType u, uc, u4;
283
284	<	if (atomType->isDirectional()) {
285	<	DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
286	<
287	<	if (daType->isDipole()) {
288	<	GenericData* data = daType->getPropertyByName("Dipole");
289	<
290	<	if (data == NULL) {
291	<	sprintf( painCave.errMsg,
292	<	"Electrostatic::addType could not find Dipole\n"
293	<	"\tparameters for atomType %s.\n",
294	<	daType->getName().c_str());
295	<	painCave.severity = OPENMD_ERROR;
296	<	painCave.isFatal = 1;
297	<	simError();
298	<	}
299	<
300	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
301	<	if (doubleData == NULL) {
302	<	sprintf( painCave.errMsg,
303	<	"Electrostatic::addType could not convert GenericData to "
304	<	"Dipole Moment\n"
305	<	"\tfor atom type %s\n", daType->getName().c_str());
306	<	painCave.severity = OPENMD_ERROR;
307	<	painCave.isFatal = 1;
308	<	simError();
309	<	}
310	<	electrostaticAtomData.is_Dipole = true;
311	<	electrostaticAtomData.dipole_moment = doubleData->getData();
284	>	// working variables for Taylor expansion:
285	>	RealType rmRc, rmRc2, rmRc3, rmRc4;
286	>
287	>	// Approximate using splines using a maximum of 0.1 Angstroms
288	>	// between points.
289	>	int nptest = int((cutoffRadius_ + 2.0) / 0.1);
290	>	np_ = (np_ > nptest) ? np_ : nptest;
291	>
292	>	// Add a 2 angstrom safety window to deal with cutoffGroups that
293	>	// have charged atoms longer than the cutoffRadius away from each
294	>	// other.  Splining is almost certainly the best choice here.
295	>	// Direct calls to erfc would be preferrable if it is a very fast
296	>	// implementation.
297	>
298	>	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_);
299	>
300	>	// Storage vectors for the computed functions
301	>	vector<RealType> rv;
302	>	vector<RealType> v01v;
303	>	vector<RealType> v11v;
304	>	vector<RealType> v21v, v22v;
305	>	vector<RealType> v31v, v32v;
306	>	vector<RealType> v41v, v42v, v43v;
307	>
308	>	/*
309	>	vector<RealType> dv01v;
310	>	vector<RealType> dv11v;
311	>	vector<RealType> dv21v, dv22v;
312	>	vector<RealType> dv31v, dv32v;
313	>	vector<RealType> dv41v, dv42v, dv43v;
314	>	*/
315	>
316	>	for (int i = 1; i < np_ + 1; i++) {
317	>	r = RealType(i) * dx;
318	>	rv.push_back(r);
319	>
320	>	ri = 1.0 / r;
321	>	ri2 = ri * ri;
322	>
323	>	r2 = r * r;
324	>	expTerm = exp(-a2 * r2);
325	>
326	>	// Taylor expansion factors (no need for factorials this way):
327	>	rmRc = r - cutoffRadius_;
328	>	rmRc2 = rmRc  * rmRc / 2.0;
329	>	rmRc3 = rmRc2 * rmRc / 3.0;
330	>	rmRc4 = rmRc3 * rmRc / 4.0;
331	>
332	>	// values of Smith's B_l functions at r:
333	>	if (screeningMethod_ == DAMPED) {
334	>	b0 = erfc(dampingAlpha_ * r) * ri;
335	>	b1 = (      b0 +     2.0*a2     * expTerm * invArootPi) * ri2;
336	>	b2 = (3.0 * b1 + pow(2.0*a2, 2) * expTerm * invArootPi) * ri2;
337	>	b3 = (5.0 * b2 + pow(2.0*a2, 3) * expTerm * invArootPi) * ri2;
338	>	b4 = (7.0 * b3 + pow(2.0*a2, 4) * expTerm * invArootPi) * ri2;
339	>	b5 = (9.0 * b4 + pow(2.0*a2, 5) * expTerm * invArootPi) * ri2;
340	>	} else {
341	>	b0 = ri;
342	>	b1 = (      b0) * ri2;
343	>	b2 = (3.0 * b1) * ri2;
344	>	b3 = (5.0 * b2) * ri2;
345	>	b4 = (7.0 * b3) * ri2;
346	>	b5 = (9.0 * b4) * ri2;
347		}
348	+
349	+	// higher derivatives of B_0 at r:
350	+	db0_1 = -r * b1;
351	+	db0_2 =     -b1 + r2 * b2;
352	+	db0_3 =          3.0*r*b2   - r2*r*b3;
353	+	db0_4 =          3.0*b2   - 6.0*r2*b3     + r2*r2*b4;
354	+	db0_5 =                    -15.0*r*b3 + 10.0*r2*r*b4 - r2*r2*r*b5;
355
356	<	if (daType->isSplitDipole()) {
357	<	GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
358	<
359	<	if (data == NULL) {
360	<	sprintf(painCave.errMsg,
361	<	"Electrostatic::addType could not find SplitDipoleDistance\n"
362	<	"\tparameter for atomType %s.\n",
363	<	daType->getName().c_str());
364	<	painCave.severity = OPENMD_ERROR;
365	<	painCave.isFatal = 1;
366	<	simError();
367	<	}
356	>	f = b0;
357	>	fc = b0c;
358	>	f0 = f - fc - rmRc*db0c_1;
359	>
360	>	g = db0_1;
361	>	gc = db0c_1;
362	>	g0 = g - gc;
363	>	g1 = g0 - rmRc *db0c_2;
364	>	g2 = g1 - rmRc2*db0c_3;
365	>	g3 = g2 - rmRc3*db0c_4;
366	>	g4 = g3 - rmRc4*db0c_5;
367	>
368	>	h = db0_2;
369	>	hc = db0c_2;
370	>	h1 = h - hc;
371	>	h2 = h1 - rmRc *db0c_3;
372	>	h3 = h2 - rmRc2*db0c_4;
373	>	h4 = h3 - rmRc3*db0c_5;
374	>
375	>	s = db0_3;
376	>	sc = db0c_3;
377	>	s2 = s - sc;
378	>	s3 = s2 - rmRc *db0c_4;
379	>	s4 = s3 - rmRc2*db0c_5;
380	>
381	>	t = db0_4;
382	>	tc = db0c_4;
383	>	t3 = t - tc;
384	>	t4 = t3 - rmRc *db0c_5;
385
386	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
387	<	if (doubleData == NULL) {
388	<	sprintf( painCave.errMsg,
348	<	"Electrostatic::addType could not convert GenericData to "
349	<	"SplitDipoleDistance for\n"
350	<	"\tatom type %s\n", daType->getName().c_str());
351	<	painCave.severity = OPENMD_ERROR;
352	<	painCave.isFatal = 1;
353	<	simError();
354	<	}
355	<	electrostaticAtomData.is_SplitDipole = true;
356	<	electrostaticAtomData.split_dipole_distance = doubleData->getData();
357	<	}
386	>	u = db0_5;
387	>	uc = db0c_5;
388	>	u4 = u - uc;
389
390	<	if (daType->isQuadrupole()) {
391	<	GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
390	>	// in what follows below, the various v functions are used for
391	>	// potentials and torques, while the w functions show up in the
392	>	// forces.
393	>
394	>	switch (summationMethod_) {
395	>	case esm_SHIFTED_FORCE:
396	>
397	>	v01 = f - fc - rmRc*gc;
398	>	v11 = g - gc - rmRc*hc;
399	>	v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
400	>	v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
401	>	v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
402	>	v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
403	>	- rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
404	>	v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
405	>	- rmRc*(sc - 3.0*(hc-gc*ric)*ric)*ric2;
406	>	v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric
407	>	- rmRc*(tc - (4.0*sc - 9.0*(hc - gc*ric)*ric)*ric)*ric;
408
409	<	if (data == NULL) {
410	<	sprintf( painCave.errMsg,
411	<	"Electrostatic::addType could not find QuadrupoleMoments\n"
412	<	"\tparameter for atomType %s.\n",
413	<	daType->getName().c_str());
414	<	painCave.severity = OPENMD_ERROR;
415	<	painCave.isFatal = 1;
416	<	simError();
417	<	}
409	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
410	>	- (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric)
411	>	- rmRc*(uc-3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
412	>
413	>	dv01 = g - gc;
414	>	dv11 = h - hc;
415	>	dv21 = (h - g*ri)*ri - (hc - gc*ric)*ric;
416	>	dv22 = (s - (h - g*ri)*ri) - (sc - (hc - gc*ric)*ric);
417	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri - (sc - 2.0*(hc-gc*ric)*ric)*ric;
418	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri)
419	>	- (tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
420	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2 - (sc - 3.0*(hc - gc*ric)*ric)*ric2;
421	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri
422	>	- (tc - (4.0*sc - 9.0*(hc-gc*ric)*ric)*ric)*ric;
423	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri)
424	>	- (uc - 3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
425
426	<	// Quadrupoles in OpenMD are set as the diagonal elements
373	<	// of the diagonalized traceless quadrupole moment tensor.
374	<	// The column vectors of the unitary matrix that diagonalizes
375	<	// the quadrupole moment tensor become the eFrame (or the
376	<	// electrostatic version of the body-fixed frame.
426	>	break;
427
428	<	Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
429	<	if (v3dData == NULL) {
430	<	sprintf( painCave.errMsg,
431	<	"Electrostatic::addType could not convert GenericData to "
432	<	"Quadrupole Moments for\n"
433	<	"\tatom type %s\n", daType->getName().c_str());
434	<	painCave.severity = OPENMD_ERROR;
435	<	painCave.isFatal = 1;
436	<	simError();
428	>	case esm_TAYLOR_SHIFTED:
429	>
430	>	v01 = f0;
431	>	v11 = g1;
432	>	v21 = g2 * ri;
433	>	v22 = h2 - v21;
434	>	v31 = (h3 - g3 * ri) * ri;
435	>	v32 = s3 - 3.0*v31;
436	>	v41 = (h4 - g4 * ri) * ri2;
437	>	v42 = s4 * ri - 3.0*v41;
438	>	v43 = t4 - 6.0*v42 - 3.0*v41;
439	>
440	>	dv01 = g0;
441	>	dv11 = h1;
442	>	dv21 = (h2 - g2*ri)*ri;
443	>	dv22 = (s2 - (h2 - g2*ri)*ri);
444	>	dv31 = (s3 - 2.0*(h3-g3*ri)*ri)*ri;
445	>	dv32 = (t3 - 3.0*(s3-2.0*(h3-g3*ri)*ri)*ri);
446	>	dv41 = (s4 - 3.0*(h4 - g4*ri)*ri)*ri2;
447	>	dv42 = (t4 - (4.0*s4 - 9.0*(h4-g4*ri)*ri)*ri)*ri;
448	>	dv43 = (u4 - 3.0*(2.0*t4 - (7.0*s4 - 15.0*(h4 - g4*ri)*ri)*ri)*ri);
449	>
450	>	break;
451	>
452	>	case esm_SHIFTED_POTENTIAL:
453	>
454	>	v01 = f - fc;
455	>	v11 = g - gc;
456	>	v21 = g*ri - gc*ric;
457	>	v22 = h - g*ri - (hc - gc*ric);
458	>	v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
459	>	v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
460	>	v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
461	>	v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;
462	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
463	>	- (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
464	>
465	>	dv01 = g;
466	>	dv11 = h;
467	>	dv21 = (h - g*ri)*ri;
468	>	dv22 = (s - (h - g*ri)*ri);
469	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
470	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
471	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
472	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
473	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
474	>
475	>	break;
476	>
477	>	case esm_SWITCHING_FUNCTION:
478	>	case esm_HARD:
479	>
480	>	v01 = f;
481	>	v11 = g;
482	>	v21 = g*ri;
483	>	v22 = h - g*ri;
484	>	v31 = (h-g*ri)*ri;
485	>	v32 = (s - 3.0*(h-g*ri)*ri);
486	>	v41 = (h - g*ri)*ri2;
487	>	v42 = (s-3.0*(h-g*ri)*ri)*ri;
488	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri);
489	>
490	>	dv01 = g;
491	>	dv11 = h;
492	>	dv21 = (h - g*ri)*ri;
493	>	dv22 = (s - (h - g*ri)*ri);
494	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
495	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
496	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
497	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
498	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
499	>
500	>	break;
501	>
502	>	case esm_REACTION_FIELD:
503	>
504	>	// following DL_POLY's lead for shifting the image charge potential:
505	>	f = b0 + preRF_ * r2;
506	>	fc = b0c + preRF_ * cutoffRadius_ * cutoffRadius_;
507	>
508	>	g = db0_1 + preRF_ * 2.0 * r;
509	>	gc = db0c_1 + preRF_ * 2.0 * cutoffRadius_;
510	>
511	>	h = db0_2 + preRF_ * 2.0;
512	>	hc = db0c_2 + preRF_ * 2.0;
513	>
514	>	v01 = f - fc;
515	>	v11 = g - gc;
516	>	v21 = g*ri - gc*ric;
517	>	v22 = h - g*ri - (hc - gc*ric);
518	>	v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
519	>	v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
520	>	v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
521	>	v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;
522	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
523	>	- (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
524	>
525	>	dv01 = g;
526	>	dv11 = h;
527	>	dv21 = (h - g*ri)*ri;
528	>	dv22 = (s - (h - g*ri)*ri);
529	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
530	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
531	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
532	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
533	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
534	>
535	>	break;
536	>
537	>	case esm_EWALD_FULL:
538	>	case esm_EWALD_PME:
539	>	case esm_EWALD_SPME:
540	>	default :
541	>	map<string, ElectrostaticSummationMethod>::iterator i;
542	>	std::string meth;
543	>	for (i = summationMap_.begin(); i != summationMap_.end(); ++i) {
544	>	if ((*i).second == summationMethod_) meth = (*i).first;
545		}
546	+	sprintf( painCave.errMsg,
547	+	"Electrostatic::initialize: electrostaticSummationMethod %s \n"
548	+	"\thas not been implemented yet. Please select one of:\n"
549	+	"\t\"hard\", \"shifted_potential\", or \"shifted_force\"\n",
550	+	meth.c_str() );
551	+	painCave.isFatal = 1;
552	+	simError();
553	+	break;
554	+	}
555	+
556	+	// Add these computed values to the storage vectors for spline creation:
557	+	v01v.push_back(v01);
558	+	v11v.push_back(v11);
559	+	v21v.push_back(v21);
560	+	v22v.push_back(v22);
561	+	v31v.push_back(v31);
562	+	v32v.push_back(v32);
563	+	v41v.push_back(v41);
564	+	v42v.push_back(v42);
565	+	v43v.push_back(v43);
566	+	/*
567	+	dv01v.push_back(dv01);
568	+	dv11v.push_back(dv11);
569	+	dv21v.push_back(dv21);
570	+	dv22v.push_back(dv22);
571	+	dv31v.push_back(dv31);
572	+	dv32v.push_back(dv32);
573	+	dv41v.push_back(dv41);
574	+	dv42v.push_back(dv42);
575	+	dv43v.push_back(dv43);
576	+	*/
577	+	}
578	+
579	+	// construct the spline structures and fill them with the values we've
580	+	// computed:
581	+
582	+	v01s = new CubicSpline();
583	+	v01s->addPoints(rv, v01v);
584	+	v11s = new CubicSpline();
585	+	v11s->addPoints(rv, v11v);
586	+	v21s = new CubicSpline();
587	+	v21s->addPoints(rv, v21v);
588	+	v22s = new CubicSpline();
589	+	v22s->addPoints(rv, v22v);
590	+	v31s = new CubicSpline();
591	+	v31s->addPoints(rv, v31v);
592	+	v32s = new CubicSpline();
593	+	v32s->addPoints(rv, v32v);
594	+	v41s = new CubicSpline();
595	+	v41s->addPoints(rv, v41v);
596	+	v42s = new CubicSpline();
597	+	v42s->addPoints(rv, v42v);
598	+	v43s = new CubicSpline();
599	+	v43s->addPoints(rv, v43v);
600	+
601	+	/*
602	+	dv01s = new CubicSpline();
603	+	dv01s->addPoints(rv, dv01v);
604	+	dv11s = new CubicSpline();
605	+	dv11s->addPoints(rv, dv11v);
606	+	dv21s = new CubicSpline();
607	+	dv21s->addPoints(rv, dv21v);
608	+	dv22s = new CubicSpline();
609	+	dv22s->addPoints(rv, dv22v);
610	+	dv31s = new CubicSpline();
611	+	dv31s->addPoints(rv, dv31v);
612	+	dv32s = new CubicSpline();
613	+	dv32s->addPoints(rv, dv32v);
614	+	dv41s = new CubicSpline();
615	+	dv41s->addPoints(rv, dv41v);
616	+	dv42s = new CubicSpline();
617	+	dv42s->addPoints(rv, dv42v);
618	+	dv43s = new CubicSpline();
619	+	dv43s->addPoints(rv, dv43v);
620	+	*/
621	+
622	+	haveElectroSplines_ = true;
623	+
624	+	initialized_ = true;
625	+	}
626	+
627	+	void Electrostatic::addType(AtomType* atomType){
628	+
629	+	ElectrostaticAtomData electrostaticAtomData;
630	+	electrostaticAtomData.is_Charge = false;
631	+	electrostaticAtomData.is_Dipole = false;
632	+	electrostaticAtomData.is_Quadrupole = false;
633	+	electrostaticAtomData.is_Fluctuating = false;
634	+
635	+	FixedChargeAdapter fca = FixedChargeAdapter(atomType);
636	+
637	+	if (fca.isFixedCharge()) {
638	+	electrostaticAtomData.is_Charge = true;
639	+	electrostaticAtomData.fixedCharge = fca.getCharge();
640	+	}
641	+
642	+	MultipoleAdapter ma = MultipoleAdapter(atomType);
643	+	if (ma.isMultipole()) {
644	+	if (ma.isDipole()) {
645	+	electrostaticAtomData.is_Dipole = true;
646	+	electrostaticAtomData.dipole = ma.getDipole();
647	+	}
648	+	if (ma.isQuadrupole()) {
649		electrostaticAtomData.is_Quadrupole = true;
650	<	electrostaticAtomData.quadrupole_moments = v3dData->getData();
650	>	electrostaticAtomData.quadrupole = ma.getQuadrupole();
651		}
652		}
653
654	<	AtomTypeProperties atp = atomType->getATP();
654	>	FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
655
656	+	if (fqa.isFluctuatingCharge()) {
657	+	electrostaticAtomData.is_Fluctuating = true;
658	+	electrostaticAtomData.electronegativity = fqa.getElectronegativity();
659	+	electrostaticAtomData.hardness = fqa.getHardness();
660	+	electrostaticAtomData.slaterN = fqa.getSlaterN();
661	+	electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
662	+	}
663	+
664		pair<map<int,AtomType*>::iterator,bool> ret;
665	<	ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
665	>	ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
666	>	atomType) );
667		if (ret.second == false) {
668		sprintf( painCave.errMsg,
669		"Electrostatic already had a previous entry with ident %d\n",
670	<	atp.ident);
670	>	atomType->getIdent() );
671		painCave.severity = OPENMD_INFO;
672		painCave.isFatal = 0;
673		simError();
674		}
675
676	<	ElectrostaticMap[atomType] = electrostaticAtomData;
676	>	ElectrostaticMap[atomType] = electrostaticAtomData;
677	>
678	>	// Now, iterate over all known types and add to the mixing map:
679	>
680	>	map<AtomType*, ElectrostaticAtomData>::iterator it;
681	>	for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
682	>	AtomType* atype2 = (*it).first;
683	>	ElectrostaticAtomData eaData2 = (*it).second;
684	>	if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
685	>
686	>	RealType a = electrostaticAtomData.slaterZeta;
687	>	RealType b = eaData2.slaterZeta;
688	>	int m = electrostaticAtomData.slaterN;
689	>	int n = eaData2.slaterN;
690	>
691	>	// Create the spline of the coulombic integral for s-type
692	>	// Slater orbitals.  Add a 2 angstrom safety window to deal
693	>	// with cutoffGroups that have charged atoms longer than the
694	>	// cutoffRadius away from each other.
695	>
696	>	RealType rval;
697	>	RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
698	>	vector<RealType> rvals;
699	>	vector<RealType> Jvals;
700	>	// don't start at i = 0, as rval = 0 is undefined for the
701	>	// slater overlap integrals.
702	>	for (int i = 1; i < np_+1; i++) {
703	>	rval = RealType(i) * dr;
704	>	rvals.push_back(rval);
705	>	Jvals.push_back(sSTOCoulInt( a, b, m, n, rval *
706	>	PhysicalConstants::angstromToBohr ) *
707	>	PhysicalConstants::hartreeToKcal );
708	>	}
709	>
710	>	CubicSpline* J = new CubicSpline();
711	>	J->addPoints(rvals, Jvals);
712	>
713	>	pair<AtomType*, AtomType*> key1, key2;
714	>	key1 = make_pair(atomType, atype2);
715	>	key2 = make_pair(atype2, atomType);
716	>
717	>	Jij[key1] = J;
718	>	Jij[key2] = J;
719	>	}
720	>	}
721	>
722		return;
723		}
724
725	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
726	<	RealType theRSW ) {
412	<	cutoffRadius_ = theECR;
413	<	rrf_ = cutoffRadius_;
414	<	rt_ = theRSW;
725	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
726	>	cutoffRadius_ = rCut;
727		haveCutoffRadius_ = true;
728		}
729	+
730		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
731		summationMethod_ = esm;
732		}
#	Line 429 | Line 742 | namespace OpenMD {
742		haveDielectric_ = true;
743		}
744
745	<	void Electrostatic::calcForce(InteractionData idat) {
745	>	void Electrostatic::calcForce(InteractionData &idat) {
746
747	<	// utility variables.  Should clean these up and use the Vector3d and
748	<	// Mat3x3d to replace as many as we can in future versions:
747	>	RealType C_a, C_b;  // Charges
748	>	Vector3d D_a, D_b;  // Dipoles (space-fixed)
749	>	Mat3x3d  Q_a, Q_b;  // Quadrupoles (space-fixed)
750
751	<	RealType q_i, q_j, mu_i, mu_j, d_i, d_j;
752	<	RealType qxx_i, qyy_i, qzz_i;
753	<	RealType qxx_j, qyy_j, qzz_j;
754	<	RealType cx_i, cy_i, cz_i;
755	<	RealType cx_j, cy_j, cz_j;
756	<	RealType cx2, cy2, cz2;
443	<	RealType ct_i, ct_j, ct_ij, a1;
444	<	RealType riji, ri, ri2, ri3, ri4;
445	<	RealType pref, vterm, epot, dudr;
446	<	RealType scale, sc2;
447	<	RealType pot_term, preVal, rfVal;
448	<	RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
449	<	RealType preSw, preSwSc;
450	<	RealType c1, c2, c3, c4;
451	<	RealType erfcVal, derfcVal;
452	<	RealType BigR;
751	>	RealType ri;                                 // Distance utility scalar
752	>	RealType rdDa, rdDb;                         // Dipole utility scalars
753	>	Vector3d rxDa, rxDb;                         // Dipole utility vectors
754	>	RealType rdQar, rdQbr, trQa, trQb;           // Quadrupole utility scalars
755	>	Vector3d Qar, Qbr, rQa, rQb, rxQar, rxQbr;   // Quadrupole utility vectors
756	>	RealType pref;
757
758	<	Vector3d Q_i, Q_j;
759	<	Vector3d ux_i, uy_i, uz_i;
760	<	Vector3d ux_j, uy_j, uz_j;
761	<	Vector3d dudux_i, duduy_i, duduz_i;
762	<	Vector3d dudux_j, duduy_j, duduz_j;
459	<	Vector3d rhatdot2, rhatc4;
460	<	Vector3d dVdr;
758	>	RealType DadDb, trQaQb, DadQbr, DbdQar;       // Cross-interaction scalars
759	>	RealType rQaQbr;
760	>	Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors
761	>	Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr;
762	>	Mat3x3d  QaQb;                                // Cross-interaction matrices
763
764	<	pair<RealType, RealType> res;
764	>	RealType U(0.0);  // Potential
765	>	Vector3d F(0.0);  // Force
766	>	Vector3d Ta(0.0); // Torque on site a
767	>	Vector3d Tb(0.0); // Torque on site b
768	>	Vector3d Ea(0.0); // Electric field at site a
769	>	Vector3d Eb(0.0); // Electric field at site b
770	>	RealType dUdCa(0.0); // fluctuating charge force at site a
771	>	RealType dUdCb(0.0); // fluctuating charge force at site a
772
773	+	// Indirect interactions mediated by the reaction field.
774	+	RealType indirect_Pot(0.0);  // Potential
775	+	Vector3d indirect_F(0.0);    // Force
776	+	Vector3d indirect_Ta(0.0);   // Torque on site a
777	+	Vector3d indirect_Tb(0.0);   // Torque on site b
778	+
779	+	// Excluded potential that is still computed for fluctuating charges
780	+	RealType excluded_Pot(0.0);
781	+
782	+	RealType rfContrib, coulInt;
783	+
784	+	// spline for coulomb integral
785	+	CubicSpline* J;
786	+
787		if (!initialized_) initialize();
788
789	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
790	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
789	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
790	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
791
792		// some variables we'll need independent of electrostatic type:
793
794	<	riji = 1.0 / idat.rij;
795	<	Vector3d rhat = idat.d  * riji;
796	<
794	>	ri = 1.0 /  *(idat.rij);
795	>	Vector3d rhat =  *(idat.d)  * ri;
796	>
797		// logicals
798
799	<	bool i_is_Charge = data1.is_Charge;
800	<	bool i_is_Dipole = data1.is_Dipole;
801	<	bool i_is_SplitDipole = data1.is_SplitDipole;
802	<	bool i_is_Quadrupole = data1.is_Quadrupole;
799	>	bool a_is_Charge = data1.is_Charge;
800	>	bool a_is_Dipole = data1.is_Dipole;
801	>	bool a_is_Quadrupole = data1.is_Quadrupole;
802	>	bool a_is_Fluctuating = data1.is_Fluctuating;
803
804	<	bool j_is_Charge = data2.is_Charge;
805	<	bool j_is_Dipole = data2.is_Dipole;
806	<	bool j_is_SplitDipole = data2.is_SplitDipole;
807	<	bool j_is_Quadrupole = data2.is_Quadrupole;
485	<
486	<	if (i_is_Charge)
487	<	q_i = data1.charge;
804	>	bool b_is_Charge = data2.is_Charge;
805	>	bool b_is_Dipole = data2.is_Dipole;
806	>	bool b_is_Quadrupole = data2.is_Quadrupole;
807	>	bool b_is_Fluctuating = data2.is_Fluctuating;
808
809	<	if (i_is_Dipole) {
810	<	mu_i = data1.dipole_moment;
491	<	uz_i = idat.eFrame1.getColumn(2);
492	<
493	<	ct_i = dot(uz_i, rhat);
494	<
495	<	if (i_is_SplitDipole)
496	<	d_i = data1.split_dipole_distance;
497	<
498	<	duduz_i = V3Zero;
499	<	}
809	>	// Obtain all of the required radial function values from the
810	>	// spline structures:
811
812	<	if (i_is_Quadrupole) {
813	<	Q_i = data1.quadrupole_moments;
814	<	qxx_i = Q_i.x();
504	<	qyy_i = Q_i.y();
505	<	qzz_i = Q_i.z();
506	<
507	<	ux_i = idat.eFrame1.getColumn(0);
508	<	uy_i = idat.eFrame1.getColumn(1);
509	<	uz_i = idat.eFrame1.getColumn(2);
510	<
511	<	cx_i = dot(ux_i, rhat);
512	<	cy_i = dot(uy_i, rhat);
513	<	cz_i = dot(uz_i, rhat);
514	<
515	<	dudux_i = V3Zero;
516	<	duduy_i = V3Zero;
517	<	duduz_i = V3Zero;
812	>	// needed for fields (and forces):
813	>	if (a_is_Charge || b_is_Charge) {
814	>	v01s->getValueAndDerivativeAt( *(idat.rij), v01, dv01);
815		}
816	+	if (a_is_Dipole || b_is_Dipole) {
817	+	v11s->getValueAndDerivativeAt( *(idat.rij), v11, dv11);
818	+	v11or = ri * v11;
819	+	}
820	+	if (a_is_Quadrupole || b_is_Quadrupole ||  (a_is_Dipole && b_is_Dipole)) {
821	+	v21s->getValueAndDerivativeAt( *(idat.rij), v21, dv21);
822	+	v22s->getValueAndDerivativeAt( *(idat.rij), v22, dv22);
823	+	v22or = ri * v22;
824	+	}
825
826	<	if (j_is_Charge)
827	<	q_j = data2.charge;
828	<
829	<	if (j_is_Dipole) {
830	<	mu_j = data2.dipole_moment;
831	<	uz_j = idat.eFrame2.getColumn(2);
832	<
527	<	ct_j = dot(uz_j, rhat);
528	<
529	<	if (j_is_SplitDipole)
530	<	d_j = data2.split_dipole_distance;
531	<
532	<	duduz_j = V3Zero;
826	>	// needed for potentials (and forces and torques):
827	>	if ((a_is_Dipole && b_is_Quadrupole) ||
828	>	(b_is_Dipole && a_is_Quadrupole)) {
829	>	v31s->getValueAndDerivativeAt( *(idat.rij), v31, dv31);
830	>	v32s->getValueAndDerivativeAt( *(idat.rij), v32, dv32);
831	>	v31or = v31 * ri;
832	>	v32or = v32 * ri;
833		}
834	<
835	<	if (j_is_Quadrupole) {
836	<	Q_j = data2.quadrupole_moments;
837	<	qxx_j = Q_j.x();
838	<	qyy_j = Q_j.y();
839	<	qzz_j = Q_j.z();
840	<
541	<	ux_j = idat.eFrame2.getColumn(0);
542	<	uy_j = idat.eFrame2.getColumn(1);
543	<	uz_j = idat.eFrame2.getColumn(2);
834	>	if (a_is_Quadrupole && b_is_Quadrupole) {
835	>	v41s->getValueAndDerivativeAt( *(idat.rij), v41, dv41);
836	>	v42s->getValueAndDerivativeAt( *(idat.rij), v42, dv42);
837	>	v43s->getValueAndDerivativeAt( *(idat.rij), v43, dv43);
838	>	v42or = v42 * ri;
839	>	v43or = v43 * ri;
840	>	}
841
842	<	cx_j = dot(ux_j, rhat);
843	<	cy_j = dot(uy_j, rhat);
844	<	cz_j = dot(uz_j, rhat);
845	<
846	<	dudux_j = V3Zero;
847	<	duduy_j = V3Zero;
848	<	duduz_j = V3Zero;
842	>	// calculate the single-site contributions (fields, etc).
843	>
844	>	if (a_is_Charge) {
845	>	C_a = data1.fixedCharge;
846	>
847	>	if (a_is_Fluctuating) {
848	>	C_a += *(idat.flucQ1);
849	>	}
850	>
851	>	if (idat.excluded) {
852	>	*(idat.skippedCharge2) += C_a;
853	>	} else {
854	>	// only do the field if we're not excluded:
855	>	Eb -= C_a *  pre11_ * dv01 * rhat;
856	>	}
857		}
858
859	<	epot = 0.0;
860	<	dVdr = V3Zero;
859	>	if (a_is_Dipole) {
860	>	D_a = *(idat.dipole1);
861	>	rdDa = dot(rhat, D_a);
862	>	rxDa = cross(rhat, D_a);
863	>	if (!idat.excluded)
864	>	Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
865	>	}
866
867	<	if (i_is_Charge) {
867	>	if (a_is_Quadrupole) {
868	>	Q_a = *(idat.quadrupole1);
869	>	trQa =  Q_a.trace();
870	>	Qar =   Q_a * rhat;
871	>	rQa = rhat * Q_a;
872	>	rdQar = dot(rhat, Qar);
873	>	rxQar = cross(rhat, Qar);
874	>	if (!idat.excluded)
875	>	Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
876	>	+ rdQar * rhat * (dv22 - 2.0*v22or));
877	>	}
878	>
879	>	if (b_is_Charge) {
880	>	C_b = data2.fixedCharge;
881
882	<	if (j_is_Charge) {
883	<	if (screeningMethod_ == DAMPED) {
884	<	// assemble the damping variables
885	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
886	<	erfcVal = res.first;
887	<	derfcVal = res.second;
888	<	c1 = erfcVal * riji;
889	<	c2 = (-derfcVal + c1) * riji;
890	<	} else {
891	<	c1 = riji;
892	<	c2 = c1 * riji;
893	<	}
882	>	if (b_is_Fluctuating)
883	>	C_b += *(idat.flucQ2);
884	>
885	>	if (idat.excluded) {
886	>	*(idat.skippedCharge1) += C_b;
887	>	} else {
888	>	// only do the field if we're not excluded:
889	>	Ea += C_b *  pre11_ * dv01 * rhat;
890	>	}
891	>	}
892	>
893	>	if (b_is_Dipole) {
894	>	D_b = *(idat.dipole2);
895	>	rdDb = dot(rhat, D_b);
896	>	rxDb = cross(rhat, D_b);
897	>	if (!idat.excluded)
898	>	Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
899	>	}
900	>
901	>	if (b_is_Quadrupole) {
902	>	Q_b = *(idat.quadrupole2);
903	>	trQb =  Q_b.trace();
904	>	Qbr =   Q_b * rhat;
905	>	rQb = rhat * Q_b;
906	>	rdQbr = dot(rhat, Qbr);
907	>	rxQbr = cross(rhat, Qbr);
908	>	if (!idat.excluded)
909	>	Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
910	>	+ rdQbr * rhat * (dv22 - 2.0*v22or));
911	>	}
912	>
913	>	if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
914	>	J = Jij[idat.atypes];
915	>	}
916	>
917	>	if (a_is_Charge) {
918	>
919	>	if (b_is_Charge) {
920	>	pref =  pre11_ * *(idat.electroMult);
921	>	U  += C_a * C_b * pref * v01;
922	>	F  += C_a * C_b * pref * dv01 * rhat;
923	>
924	>	// If this is an excluded pair, there are still indirect
925	>	// interactions via the reaction field we must worry about:
926
927	<	preVal = idat.electroMult * pre11_ * q_i * q_j;
927	>	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
928	>	rfContrib = preRF_ * pref * C_a * C_b * *(idat.r2);
929	>	indirect_Pot += rfContrib;
930	>	indirect_F   += rfContrib * 2.0 * ri * rhat;
931	>	}
932
933	<	if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
934	<	vterm = preVal * (c1 - c1c_);
935	<	dudr  = -idat.sw * preVal * c2;
933	>	// Fluctuating charge forces are handled via Coulomb integrals
934	>	// for excluded pairs (i.e. those connected via bonds) and
935	>	// with the standard charge-charge interaction otherwise.
936
937	<	} else if (summationMethod_ == esm_SHIFTED_FORCE)  {
938	<	vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
939	<	dudr  = idat.sw * preVal * (c2c_ - c2);
937	>	if (idat.excluded) {
938	>	if (a_is_Fluctuating || b_is_Fluctuating) {
939	>	coulInt = J->getValueAt( *(idat.rij) );
940	>	if (a_is_Fluctuating)  dUdCa += coulInt * C_b;
941	>	if (b_is_Fluctuating)  dUdCb += coulInt * C_a;
942	>	excluded_Pot += C_a * C_b * coulInt;
943	>	}
944	>	} else {
945	>	if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
946	>	if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
947	>	}
948	>	}
949
950	<	} else if (summationMethod_ == esm_REACTION_FIELD) {
951	<	rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
952	<	vterm = preVal * ( riji + rfVal );
953	<	dudr  = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
950	>	if (b_is_Dipole) {
951	>	pref =  pre12_ * *(idat.electroMult);
952	>	U  += C_a * pref * v11 * rdDb;
953	>	F  += C_a * pref * ((dv11 - v11or) * rdDb * rhat + v11or * D_b);
954	>	Tb += C_a * pref * v11 * rxDb;
955
956	<	} else {
588	<	vterm = preVal * riji * erfcVal;
956	>	if (a_is_Fluctuating) dUdCa += pref * v11 * rdDb;
957
958	<	dudr  = - idat.sw * preVal * c2;
958	>	// Even if we excluded this pair from direct interactions, we
959	>	// still have the reaction-field-mediated charge-dipole
960	>	// interaction:
961
962	+	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
963	+	rfContrib = C_a * pref * preRF_ * 2.0 * *(idat.rij);
964	+	indirect_Pot += rfContrib * rdDb;
965	+	indirect_F   += rfContrib * D_b / (*idat.rij);
966	+	indirect_Tb  += C_a * pref * preRF_ * rxDb;
967		}
968	<
594	<	idat.vpair += vterm;
595	<	epot += idat.sw * vterm;
968	>	}
969
970	<	dVdr += dudr * rhat;
970	>	if (b_is_Quadrupole) {
971	>	pref = pre14_ * *(idat.electroMult);
972	>	U  +=  C_a * pref * (v21 * trQb + v22 * rdQbr);
973	>	F  +=  C_a * pref * (trQb * dv21 * rhat + 2.0 * Qbr * v22or);
974	>	F  +=  C_a * pref * rdQbr * rhat * (dv22 - 2.0*v22or);
975	>	Tb +=  C_a * pref * 2.0 * rxQbr * v22;
976	>
977	>	if (a_is_Fluctuating) dUdCa += pref * (v21 * trQb + v22 * rdQbr);
978		}
979	+	}
980
981	<	if (j_is_Dipole) {
601	<	// pref is used by all the possible methods
602	<	pref = idat.electroMult * pre12_ * q_i * mu_j;
603	<	preSw = idat.sw * pref;
981	>	if (a_is_Dipole) {
982
983	<	if (summationMethod_ == esm_REACTION_FIELD) {
984	<	ri2 = riji * riji;
607	<	ri3 = ri2 * riji;
608	<
609	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
610	<	idat.vpair += vterm;
611	<	epot += idat.sw * vterm;
983	>	if (b_is_Charge) {
984	>	pref = pre12_ * *(idat.electroMult);
985
986	<	dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
987	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);
986	>	U  -= C_b * pref * v11 * rdDa;
987	>	F  -= C_b * pref * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
988	>	Ta -= C_b * pref * v11 * rxDa;
989
990	<	} else {
617	<	// determine the inverse r used if we have split dipoles
618	<	if (j_is_SplitDipole) {
619	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
620	<	ri = 1.0 / BigR;
621	<	scale = idat.rij * ri;
622	<	} else {
623	<	ri = riji;
624	<	scale = 1.0;
625	<	}
626	<
627	<	sc2 = scale * scale;
990	>	if (b_is_Fluctuating) dUdCb -= pref * v11 * rdDa;
991
992	<	if (screeningMethod_ == DAMPED) {
993	<	// assemble the damping variables
994	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
995	<	erfcVal = res.first;
996	<	derfcVal = res.second;
997	<	c1 = erfcVal * ri;
998	<	c2 = (-derfcVal + c1) * ri;
999	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
1000	<	} else {
1001	<	c1 = ri;
639	<	c2 = c1 * ri;
640	<	c3 = 3.0 * c2 * ri;
641	<	}
642	<
643	<	c2ri = c2 * ri;
992	>	// Even if we excluded this pair from direct interactions,
993	>	// we still have the reaction-field-mediated charge-dipole
994	>	// interaction:
995	>	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
996	>	rfContrib = C_b * pref * preRF_ * 2.0 * *(idat.rij);
997	>	indirect_Pot -= rfContrib * rdDa;
998	>	indirect_F   -= rfContrib * D_a / (*idat.rij);
999	>	indirect_Ta  -= C_b * pref * preRF_ * rxDa;
1000	>	}
1001	>	}
1002
1003	<	// calculate the potential
1004	<	pot_term =  scale * c2;
1005	<	vterm = -pref * ct_j * pot_term;
1006	<	idat.vpair += vterm;
649	<	epot += idat.sw * vterm;
650	<
651	<	// calculate derivatives for forces and torques
1003	>	if (b_is_Dipole) {
1004	>	pref = pre22_ * *(idat.electroMult);
1005	>	DadDb = dot(D_a, D_b);
1006	>	DaxDb = cross(D_a, D_b);
1007
1008	<	dVdr += -preSw * (uz_j * c2ri - ct_j * rhat * sc2 * c3);
1009	<	duduz_j += -preSw * pot_term * rhat;
1008	>	U  -= pref * (DadDb * v21 + rdDa * rdDb * v22);
1009	>	F  -= pref * (dv21 * DadDb * rhat + v22or * (rdDb * D_a + rdDa * D_b));
1010	>	F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1011	>	Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1012	>	Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
1013
1014	+	// Even if we excluded this pair from direct interactions, we
1015	+	// still have the reaction-field-mediated dipole-dipole
1016	+	// interaction:
1017	+	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
1018	+	rfContrib = -pref * preRF_ * 2.0;
1019	+	indirect_Pot += rfContrib * DadDb;
1020	+	indirect_Ta  += rfContrib * DaxDb;
1021	+	indirect_Tb  -= rfContrib * DaxDb;
1022		}
1023		}
1024
1025	<	if (j_is_Quadrupole) {
1026	<	// first precalculate some necessary variables
1027	<	cx2 = cx_j * cx_j;
1028	<	cy2 = cy_j * cy_j;
1029	<	cz2 = cz_j * cz_j;
664	<	pref =  idat.electroMult * pre14_ * q_i * one_third_;
665	<
666	<	if (screeningMethod_ == DAMPED) {
667	<	// assemble the damping variables
668	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
669	<	erfcVal = res.first;
670	<	derfcVal = res.second;
671	<	c1 = erfcVal * riji;
672	<	c2 = (-derfcVal + c1) * riji;
673	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
674	<	c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * riji * riji;
675	<	} else {
676	<	c1 = riji;
677	<	c2 = c1 * riji;
678	<	c3 = 3.0 * c2 * riji;
679	<	c4 = 5.0 * c3 * riji * riji;
680	<	}
1025	>	if (b_is_Quadrupole) {
1026	>	pref = pre24_ * *(idat.electroMult);
1027	>	DadQb = D_a * Q_b;
1028	>	DadQbr = dot(D_a, Qbr);
1029	>	DaxQbr = cross(D_a, Qbr);
1030
1031	<	// precompute variables for convenience
1032	<	preSw = idat.sw * pref;
1033	<	c2ri = c2 * riji;
1034	<	c3ri = c3 * riji;
1035	<	c4rij = c4 * idat.rij;
1036	<	rhatdot2 = 2.0 * rhat * c3;
1037	<	rhatc4 = rhat * c4rij;
1031	>	U  -= pref * ((trQb*rdDa + 2.0*DadQbr)*v31 + rdDa*rdQbr*v32);
1032	>	F  -= pref * (trQb*D_a + 2.0*DadQb) * v31or;
1033	>	F  -= pref * (trQb*rdDa + 2.0*DadQbr) * (dv31-v31or) * rhat;
1034	>	F  -= pref * (D_a*rdQbr + 2.0*rdDa*rQb) * v32or;
1035	>	F  -= pref * (rdDa * rdQbr * rhat * (dv32-3.0*v32or));
1036	>	Ta += pref * ((-trQb*rxDa + 2.0 * DaxQbr)*v31 - rxDa*rdQbr*v32);
1037	>	Tb += pref * ((2.0*cross(DadQb, rhat) - 2.0*DaxQbr)*v31
1038	>	- 2.0*rdDa*rxQbr*v32);
1039	>	}
1040	>	}
1041
1042	<	// calculate the potential
1043	<	pot_term = ( qxx_j * (cx2*c3 - c2ri) +
1044	<	qyy_j * (cy2*c3 - c2ri) +
1045	<	qzz_j * (cz2*c3 - c2ri) );
1046	<	vterm = pref * pot_term;
1047	<	idat.vpair += vterm;
1048	<	epot += idat.sw * vterm;
697	<
698	<	// calculate derivatives for the forces and torques
1042	>	if (a_is_Quadrupole) {
1043	>	if (b_is_Charge) {
1044	>	pref = pre14_ * *(idat.electroMult);
1045	>	U  += C_b * pref * (v21 * trQa + v22 * rdQar);
1046	>	F  += C_b * pref * (trQa * rhat * dv21 + 2.0 * Qar * v22or);
1047	>	F  += C_b * pref * rdQar * rhat * (dv22 - 2.0*v22or);
1048	>	Ta += C_b * pref * 2.0 * rxQar * v22;
1049
1050	<	dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
701	<	qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
702	<	qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
703	<
704	<	dudux_j += preSw * qxx_j * cx_j * rhatdot2;
705	<	duduy_j += preSw * qyy_j * cy_j * rhatdot2;
706	<	duduz_j += preSw * qzz_j * cz_j * rhatdot2;
1050	>	if (b_is_Fluctuating) dUdCb += pref * (v21 * trQa + v22 * rdQar);
1051		}
1052	<	}
1053	<
1054	<	if (i_is_Dipole) {
1052	>	if (b_is_Dipole) {
1053	>	pref = pre24_ * *(idat.electroMult);
1054	>	DbdQa = D_b * Q_a;
1055	>	DbdQar = dot(D_b, Qar);
1056	>	DbxQar = cross(D_b, Qar);
1057
1058	<	if (j_is_Charge) {
1059	<	// variables used by all the methods
1060	<	pref = idat.electroMult * pre12_ * q_j * mu_i;
1061	<	preSw = idat.sw * pref;
1058	>	U  += pref * ((trQa*rdDb + 2.0*DbdQar)*v31 + rdDb*rdQar*v32);
1059	>	F  += pref * (trQa*D_b + 2.0*DbdQa) * v31or;
1060	>	F  += pref * (trQa*rdDb + 2.0*DbdQar) * (dv31-v31or) * rhat;
1061	>	F  += pref * (D_b*rdQar + 2.0*rdDb*rQa) * v32or;
1062	>	F  += pref * (rdDb * rdQar * rhat * (dv32-3.0*v32or));
1063	>	Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31
1064	>	+ 2.0*rdDb*rxQar*v32);
1065	>	Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1066	>	}
1067	>	if (b_is_Quadrupole) {
1068	>	pref = pre44_ * *(idat.electroMult);  // yes
1069	>	QaQb = Q_a * Q_b;
1070	>	trQaQb = QaQb.trace();
1071	>	rQaQb = rhat * QaQb;
1072	>	QaQbr = QaQb * rhat;
1073	>	QaxQb = cross(Q_a, Q_b);
1074	>	rQaQbr = dot(rQa, Qbr);
1075	>	rQaxQbr = cross(rQa, Qbr);
1076	>
1077	>	U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1078	>	U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
1079	>	U  += pref * (rdQar * rdQbr) * v43;
1080
1081	<	if (summationMethod_ == esm_REACTION_FIELD) {
1081	>	F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*dv41;
1082	>	F  += pref*rhat*(trQa*rdQbr + trQb*rdQar + 4.0*rQaQbr)*(dv42-2.0*v42or);
1083	>	F  += pref * rhat * (rdQar * rdQbr)*(dv43 - 4.0*v43or);
1084
1085	<	ri2 = riji * riji;
1086	<	ri3 = ri2 * riji;
1085	>	F  += pref * 2.0 * (trQb*rQa + trQa*rQb) * v42or;
1086	>	F  += pref * 4.0 * (rQaQb + QaQbr) * v42or;
1087	>	F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb) * v43or;
1088
1089	<	vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
1090	<	idat.vpair += vterm;
1091	<	epot += idat.sw * vterm;
1092	<
1093	<	dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
727	<
728	<	duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
729	<
730	<	} else {
731	<
732	<	// determine inverse r if we are using split dipoles
733	<	if (i_is_SplitDipole) {
734	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
735	<	ri = 1.0 / BigR;
736	<	scale = idat.rij * ri;
737	<	} else {
738	<	ri = riji;
739	<	scale = 1.0;
740	<	}
741	<
742	<	sc2 = scale * scale;
743	<
744	<	if (screeningMethod_ == DAMPED) {
745	<	// assemble the damping variables
746	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
747	<	erfcVal = res.first;
748	<	derfcVal = res.second;
749	<	c1 = erfcVal * ri;
750	<	c2 = (-derfcVal + c1) * ri;
751	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
752	<	} else {
753	<	c1 = ri;
754	<	c2 = c1 * ri;
755	<	c3 = 3.0 * c2 * ri;
756	<	}
757	<
758	<	c2ri = c2 * ri;
759	<
760	<	// calculate the potential
761	<	pot_term = c2 * scale;
762	<	vterm = pref * ct_i * pot_term;
763	<	idat.vpair += vterm;
764	<	epot += idat.sw * vterm;
1089	>	Ta += pref * (- 4.0 * QaxQb * v41);
1090	>	Ta += pref * (- 2.0 * trQb * cross(rQa, rhat)
1091	>	+ 4.0 * cross(rhat, QaQbr)
1092	>	- 4.0 * rQaxQbr) * v42;
1093	>	Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
1094
766	–	// calculate derivatives for the forces and torques
767	–	dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
768	–	duduz_i += preSw * pot_term * rhat;
769	–	}
770	–	}
1095
1096	<	if (j_is_Dipole) {
1097	<	// variables used by all methods
1098	<	ct_ij = dot(uz_i, uz_j);
1096	>	Tb += pref * (+ 4.0 * QaxQb * v41);
1097	>	Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1098	>	- 4.0 * cross(rQaQb, rhat)
1099	>	+ 4.0 * rQaxQbr) * v42;
1100	>	// Possible replacement for line 2 above:
1101	>	//             + 4.0 * cross(rhat, QbQar)
1102
1103	<	pref = idat.electroMult * pre22_ * mu_i * mu_j;
777	<	preSw = idat.sw * pref;
1103	>	Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1104
1105	<	if (summationMethod_ == esm_REACTION_FIELD) {
780	<	ri2 = riji * riji;
781	<	ri3 = ri2 * riji;
782	<	ri4 = ri2 * ri2;
783	<
784	<	vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
785	<	preRF2_ * ct_ij );
786	<	idat.vpair += vterm;
787	<	epot += idat.sw * vterm;
788	<
789	<	a1 = 5.0 * ct_i * ct_j - ct_ij;
790	<
791	<	dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
792	<
793	<	duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
794	<	duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
795	<
796	<	} else {
797	<
798	<	if (i_is_SplitDipole) {
799	<	if (j_is_SplitDipole) {
800	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
801	<	} else {
802	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
803	<	}
804	<	ri = 1.0 / BigR;
805	<	scale = idat.rij * ri;
806	<	} else {
807	<	if (j_is_SplitDipole) {
808	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
809	<	ri = 1.0 / BigR;
810	<	scale = idat.rij * ri;
811	<	} else {
812	<	ri = riji;
813	<	scale = 1.0;
814	<	}
815	<	}
816	<	if (screeningMethod_ == DAMPED) {
817	<	// assemble damping variables
818	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
819	<	erfcVal = res.first;
820	<	derfcVal = res.second;
821	<	c1 = erfcVal * ri;
822	<	c2 = (-derfcVal + c1) * ri;
823	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
824	<	c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * ri * ri;
825	<	} else {
826	<	c1 = ri;
827	<	c2 = c1 * ri;
828	<	c3 = 3.0 * c2 * ri;
829	<	c4 = 5.0 * c3 * ri * ri;
830	<	}
831	<
832	<	// precompute variables for convenience
833	<	sc2 = scale * scale;
834	<	cti3 = ct_i * sc2 * c3;
835	<	ctj3 = ct_j * sc2 * c3;
836	<	ctidotj = ct_i * ct_j * sc2;
837	<	preSwSc = preSw * scale;
838	<	c2ri = c2 * ri;
839	<	c3ri = c3 * ri;
840	<	c4rij = c4 * idat.rij;
841	<
842	<	// calculate the potential
843	<	pot_term = (ct_ij * c2ri - ctidotj * c3);
844	<	vterm = pref * pot_term;
845	<	idat.vpair += vterm;
846	<	epot += idat.sw * vterm;
847	<
848	<	// calculate derivatives for the forces and torques
849	<	dVdr += preSwSc * ( ctidotj * rhat * c4rij  -
850	<	(ct_i*uz_j + ct_j*uz_i + ct_ij*rhat) * c3ri);
851	<
852	<	duduz_i += preSw * (uz_j * c2ri - ctj3 * rhat);
853	<	duduz_j += preSw * (uz_i * c2ri - cti3 * rhat);
854	<	}
1105	>	//  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1106		}
1107		}
1108
1109	<	if (i_is_Quadrupole) {
1110	<	if (j_is_Charge) {
1111	<	// precompute some necessary variables
1112	<	cx2 = cx_i * cx_i;
862	<	cy2 = cy_i * cy_i;
863	<	cz2 = cz_i * cz_i;
1109	>	if (idat.doElectricField) {
1110	>	*(idat.eField1) += Ea * *(idat.electroMult);
1111	>	*(idat.eField2) += Eb * *(idat.electroMult);
1112	>	}
1113
1114	<	pref = idat.electroMult * pre14_ * q_j * one_third_;
1114	>	if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
1115	>	if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw);
1116
1117	<	if (screeningMethod_ == DAMPED) {
1118	<	// assemble the damping variables
1119	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
1120	<	erfcVal = res.first;
1121	<	derfcVal = res.second;
1122	<	c1 = erfcVal * riji;
1123	<	c2 = (-derfcVal + c1) * riji;
1124	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
875	<	c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * riji * riji;
876	<	} else {
877	<	c1 = riji;
878	<	c2 = c1 * riji;
879	<	c3 = 3.0 * c2 * riji;
880	<	c4 = 5.0 * c3 * riji * riji;
881	<	}
882	<
883	<	// precompute some variables for convenience
884	<	preSw = idat.sw * pref;
885	<	c2ri = c2 * riji;
886	<	c3ri = c3 * riji;
887	<	c4rij = c4 * idat.rij;
888	<	rhatdot2 = 2.0 * rhat * c3;
889	<	rhatc4 = rhat * c4rij;
1117	>	if (!idat.excluded) {
1118	>
1119	>	*(idat.vpair) += U;
1120	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += U * *(idat.sw);
1121	>	*(idat.f1) += F * *(idat.sw);
1122	>
1123	>	if (a_is_Dipole || a_is_Quadrupole)
1124	>	*(idat.t1) += Ta * *(idat.sw);
1125
1126	<	// calculate the potential
1127	<	pot_term = ( qxx_i * (cx2 * c3 - c2ri) +
1128	<	qyy_i * (cy2 * c3 - c2ri) +
1129	<	qzz_i * (cz2 * c3 - c2ri) );
895	<
896	<	vterm = pref * pot_term;
897	<	idat.vpair += vterm;
898	<	epot += idat.sw * vterm;
1126	>	if (b_is_Dipole || b_is_Quadrupole)
1127	>	*(idat.t2) += Tb * *(idat.sw);
1128	>
1129	>	} else {
1130
1131	<	// calculate the derivatives for the forces and torques
1131	>	// only accumulate the forces and torques resulting from the
1132	>	// indirect reaction field terms.
1133
1134	<	dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
1135	<	qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
1136	<	qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
1137	<
1138	<	dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1139	<	duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1140	<	duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1141	<	}
1134	>	*(idat.vpair) += indirect_Pot;
1135	>	(*(idat.excludedPot))[ELECTROSTATIC_FAMILY] +=  excluded_Pot;
1136	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += *(idat.sw) * indirect_Pot;
1137	>	*(idat.f1) += *(idat.sw) * indirect_F;
1138	>
1139	>	if (a_is_Dipole || a_is_Quadrupole)
1140	>	*(idat.t1) += *(idat.sw) * indirect_Ta;
1141	>
1142	>	if (b_is_Dipole || b_is_Quadrupole)
1143	>	*(idat.t2) += *(idat.sw) * indirect_Tb;
1144		}
911	–
912	–	idat.pot += epot;
913	–	idat.f1 += dVdr;
914	–
915	–	if (i_is_Dipole || i_is_Quadrupole)
916	–	idat.t1 -= cross(uz_i, duduz_i);
917	–	if (i_is_Quadrupole) {
918	–	idat.t1 -= cross(ux_i, dudux_i);
919	–	idat.t1 -= cross(uy_i, duduy_i);
920	–	}
921	–
922	–	if (j_is_Dipole || j_is_Quadrupole)
923	–	idat.t2 -= cross(uz_j, duduz_j);
924	–	if (j_is_Quadrupole) {
925	–	idat.t2 -= cross(uz_j, dudux_j);
926	–	idat.t2 -= cross(uz_j, duduy_j);
927	–	}
928	–
1145		return;
1146	<	}
1146	>	}
1147	>
1148	>	void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1149
932	–	void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
933	–
1150		if (!initialized_) initialize();
1151	+
1152	+	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1153
936	–	ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
937	–	ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
938	–
1154		// logicals
1155	<
1156	<	bool i_is_Charge = data1.is_Charge;
1157	<	bool i_is_Dipole = data1.is_Dipole;
1158	<
1159	<	bool j_is_Charge = data2.is_Charge;
945	<	bool j_is_Dipole = data2.is_Dipole;
946	<
947	<	RealType q_i, q_j;
1155	>	bool i_is_Charge = data.is_Charge;
1156	>	bool i_is_Dipole = data.is_Dipole;
1157	>	bool i_is_Fluctuating = data.is_Fluctuating;
1158	>	RealType C_a = data.fixedCharge;
1159	>	RealType self, preVal, DadDa;
1160
1161	<	// The skippedCharge computation is needed by the real-space cutoff methods
1162	<	// (i.e. shifted force and shifted potential)
1163	<
1164	<	if (i_is_Charge) {
1165	<	q_i = data1.charge;
1166	<	skdat.skippedCharge2 += q_i;
1161	>	if (i_is_Fluctuating) {
1162	>	C_a += *(sdat.flucQ);
1163	>	// dVdFQ is really a force, so this is negative the derivative
1164	>	*(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1165	>	(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1166	>	(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1167		}
1168
1169	<	if (j_is_Charge) {
1170	<	q_j = data2.charge;
959	<	skdat.skippedCharge1 += q_j;
960	<	}
961	<
962	<	// the rest of this function should only be necessary for reaction field.
963	<
964	<	if (summationMethod_ == esm_REACTION_FIELD) {
965	<	RealType riji, ri2, ri3;
966	<	RealType q_i, mu_i, ct_i;
967	<	RealType q_j, mu_j, ct_j;
968	<	RealType preVal, rfVal, vterm, dudr, pref, myPot;
969	<	Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
970	<
971	<	// some variables we'll need independent of electrostatic type:
1169	>	switch (summationMethod_) {
1170	>	case esm_REACTION_FIELD:
1171
973	–	riji = 1.0 / skdat.rij;
974	–	rhat = skdat.d  * riji;
975	–
976	–	if (i_is_Dipole) {
977	–	mu_i = data1.dipole_moment;
978	–	uz_i = skdat.eFrame1.getColumn(2);
979	–	ct_i = dot(uz_i, rhat);
980	–	duduz_i = V3Zero;
981	–	}
982	–
983	–	if (j_is_Dipole) {
984	–	mu_j = data2.dipole_moment;
985	–	uz_j = skdat.eFrame2.getColumn(2);
986	–	ct_j = dot(uz_j, rhat);
987	–	duduz_j = V3Zero;
988	–	}
989	–
1172		if (i_is_Charge) {
1173	<	if (j_is_Charge) {
1174	<	preVal = skdat.electroMult * pre11_ * q_i * q_j;
1175	<	rfVal = preRF_ * skdat.rij * skdat.rij;
1176	<	vterm = preVal * rfVal;
1177	<	myPot += skdat.sw * vterm;
996	<	dudr  = skdat.sw * preVal * 2.0 * rfVal * riji;
997	<	dVdr += dudr * rhat;
998	<	}
999	<
1000	<	if (j_is_Dipole) {
1001	<	ri2 = riji * riji;
1002	<	ri3 = ri2 * riji;
1003	<	pref = skdat.electroMult * pre12_ * q_i * mu_j;
1004	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1005	<	myPot += skdat.sw * vterm;
1006	<	dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1007	<	duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1008	<	}
1173	>	// Self potential [see Wang and Hermans, "Reaction Field
1174	>	// Molecular Dynamics Simulation with Friedman’s Image Charge
1175	>	// Method," J. Phys. Chem. 99, 12001-12007 (1995).]
1176	>	preVal = pre11_ * preRF_ * C_a * C_a;
1177	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal / cutoffRadius_;
1178		}
1179	+
1180		if (i_is_Dipole) {
1181	<	if (j_is_Charge) {
1182	<	ri2 = riji * riji;
1013	<	ri3 = ri2 * riji;
1014	<	pref = skdat.electroMult * pre12_ * q_j * mu_i;
1015	<	vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1016	<	myPot += skdat.sw * vterm;
1017	<	dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
1018	<	duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1019	<	}
1181	>	DadDa = data.dipole.lengthSquare();
1182	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa;
1183		}
1184
1185	<	// accumulate the forces and torques resulting from the self term
1023	<	skdat.pot += myPot;
1024	<	skdat.f1 += dVdr;
1185	>	break;
1186
1187	<	if (i_is_Dipole)
1188	<	skdat.t1 -= cross(uz_i, duduz_i);
1189	<	if (j_is_Dipole)
1190	<	skdat.t2 -= cross(uz_j, duduz_j);
1191	<	}
1031	<	}
1032	<
1033	<	void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1034	<	RealType mu1, preVal, chg1, self;
1035	<
1036	<	if (!initialized_) initialize();
1037	<
1038	<	ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1039	<
1040	<	// logicals
1041	<
1042	<	bool i_is_Charge = data.is_Charge;
1043	<	bool i_is_Dipole = data.is_Dipole;
1044	<
1045	<	if (summationMethod_ == esm_REACTION_FIELD) {
1046	<	if (i_is_Dipole) {
1047	<	mu1 = data.dipole_moment;
1048	<	preVal = pre22_ * preRF2_ * mu1 * mu1;
1049	<	scdat.pot -= 0.5 * preVal;
1050	<
1051	<	// The self-correction term adds into the reaction field vector
1052	<	Vector3d uz_i = scdat.eFrame.getColumn(2);
1053	<	Vector3d ei = preVal * uz_i;
1054	<
1055	<	// This looks very wrong.  A vector crossed with itself is zero.
1056	<	scdat.t -= cross(uz_i, ei);
1187	>	case esm_SHIFTED_FORCE:
1188	>	case esm_SHIFTED_POTENTIAL:
1189	>	if (i_is_Charge) {
1190	>	self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1191	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1192		}
1193	<	} else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1194	<	if (i_is_Charge) {
1195	<	chg1 = data.charge;
1061	<	if (screeningMethod_ == DAMPED) {
1062	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1063	<	} else {
1064	<	self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1065	<	}
1066	<	scdat.pot += self;
1067	<	}
1193	>	break;
1194	>	default:
1195	>	break;
1196		}
1197		}
1198	<
1199	<	RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1198	>
1199	>	RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1200		// This seems to work moderately well as a default.  There's no
1201		// inherent scale for 1/r interactions that we can standardize.
1202		// 12 angstroms seems to be a reasonably good guess for most

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines