ViewVC Help

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

root/OpenMD/trunk/src/nonbonded/Electrostatic.cpp

Comparing:
branches/development/src/nonbonded/Electrostatic.cpp (file contents), Revision 1734 by jmichalk, Tue Jun 5 17:49:36 2012 UTC vs.
trunk/src/nonbonded/Electrostatic.cpp (file contents), Revision 1895 by gezelter, Mon Jul 1 21:09:37 2013 UTC

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

Diff Legend

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