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

Diff Legend

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