ViewVC Help

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

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

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

Diff Legend

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