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 1767 by gezelter, Fri Jul 6 22:01:58 2012 UTC vs.
trunk/src/nonbonded/Electrostatic.cpp (file contents), Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

#	Line 35 | Line 35
35		*
36		* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37		* [2]  Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).
38	<	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
38	>	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39		* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40		* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43	+	#ifdef IS_MPI
44	+	#include <mpi.h>
45	+	#endif
46	+
47		#include <stdio.h>
48		#include <string.h>
49
50		#include <cmath>
51	+	#include <numeric>
52		#include "nonbonded/Electrostatic.hpp"
53		#include "utils/simError.h"
54		#include "types/NonBondedInteractionType.hpp"
#	Line 54 | Line 59
59		#include "nonbonded/SlaterIntegrals.hpp"
60		#include "utils/PhysicalConstants.hpp"
61		#include "math/erfc.hpp"
62	+	#include "math/SquareMatrix.hpp"
63	+	#include "primitives/Molecule.hpp"
64	+	#include "flucq/FluctuatingChargeForces.hpp"
65
66		namespace OpenMD {
67
68		Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
61	–	forceField_(NULL), info_(NULL),
69		haveCutoffRadius_(false),
70		haveDampingAlpha_(false),
71		haveDielectric_(false),
72	<	haveElectroSpline_(false)
73	<	{}
72	>	haveElectroSplines_(false),
73	>	info_(NULL), forceField_(NULL)
74	>
75	>	{
76	>	flucQ_ = new FluctuatingChargeForces(info_);
77	>	}
78
79	+	void Electrostatic::setForceField(ForceField *ff) {
80	+	forceField_ = ff;
81	+	flucQ_->setForceField(forceField_);
82	+	}
83	+
84	+	void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) {
85	+	simTypes_ = simtypes;
86	+	flucQ_->setSimulatedAtomTypes(simTypes_);
87	+	}
88	+
89		void Electrostatic::initialize() {
90
91		Globals* simParams_ = info_->getSimParams();
#	Line 74 | Line 95 | namespace OpenMD {
95		summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
96		summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
97		summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;
98	+	summationMap_["TAYLOR_SHIFTED"]     = esm_TAYLOR_SHIFTED;
99		summationMap_["REACTION_FIELD"]     = esm_REACTION_FIELD;
100		summationMap_["EWALD_FULL"]         = esm_EWALD_FULL;
101		summationMap_["EWALD_PME"]          = esm_EWALD_PME;
#	Line 88 | Line 110 | namespace OpenMD {
110		// Charge-Dipole, assuming charges are measured in electrons, and
111		// dipoles are measured in debyes
112		pre12_ = 69.13373;
113	<	// Dipole-Dipole, assuming dipoles are measured in debyes
113	>	// Dipole-Dipole, assuming dipoles are measured in Debye
114		pre22_ = 14.39325;
115		// Charge-Quadrupole, assuming charges are measured in electrons, and
116		// quadrupoles are measured in 10^-26 esu cm^2
117	<	// This unit is also known affectionately as an esu centi-barn.
117	>	// This unit is also known affectionately as an esu centibarn.
118		pre14_ = 69.13373;
119	<
119	>	// Dipole-Quadrupole, assuming dipoles are measured in debyes and
120	>	// quadrupoles in esu centibarns:
121	>	pre24_ = 14.39325;
122	>	// Quadrupole-Quadrupole, assuming esu centibarns:
123	>	pre44_ = 14.39325;
124	>
125		// conversions for the simulation box dipole moment
126		chargeToC_ = 1.60217733e-19;
127		angstromToM_ = 1.0e-10;
128		debyeToCm_ = 3.33564095198e-30;
129
130	<	// number of points for electrostatic splines
130	>	// Default number of points for electrostatic splines
131		np_ = 100;
132
133		// variables to handle different summation methods for long-range
#	Line 108 | Line 135 | namespace OpenMD {
135		summationMethod_ = esm_HARD;
136		screeningMethod_ = UNDAMPED;
137		dielectric_ = 1.0;
111	–	one_third_ = 1.0 / 3.0;
138
139		// check the summation method:
140		if (simParams_->haveElectrostaticSummationMethod()) {
#	Line 124 | Line 150 | namespace OpenMD {
150		"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
151		"\t(Input file specified %s .)\n"
152		"\telectrostaticSummationMethod must be one of: \"hard\",\n"
153	<	"\t\"shifted_potential\", \"shifted_force\", or \n"
154	<	"\t\"reaction_field\".\n", myMethod.c_str() );
153	>	"\t\"shifted_potential\", \"shifted_force\",\n"
154	>	"\t\"taylor_shifted\", or \"reaction_field\".\n",
155	>	myMethod.c_str() );
156		painCave.isFatal = 1;
157		simError();
158		}
#	Line 184 | Line 211 | namespace OpenMD {
211		simError();
212		}
213
214	<	if (screeningMethod_ == DAMPED) {
214	>	if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) {
215		if (!simParams_->haveDampingAlpha()) {
216		// first set a cutoff dependent alpha value
217		// we assume alpha depends linearly with rcut from 0 to 20.5 ang
218		dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
219	<	if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
193	<
219	>	if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
220		// throw warning
221		sprintf( painCave.errMsg,
222		"Electrostatic::initialize: dampingAlpha was not specified in the\n"
#	Line 206 | Line 232 | namespace OpenMD {
232		haveDampingAlpha_ = true;
233		}
234
235	<	// find all of the Electrostatic atom Types:
236	<	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
237	<	ForceField::AtomTypeContainer::MapTypeIterator i;
238	<	AtomType* at;
239	<
240	<	for (at = atomTypes->beginType(i); at != NULL;
241	<	at = atomTypes->nextType(i)) {
242	<
243	<	if (at->isElectrostatic())
218	<	addType(at);
219	<	}
220	<
221	<	cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
222	<	rcuti_ = 1.0 / cutoffRadius_;
223	<	rcuti2_ = rcuti_ * rcuti_;
224	<	rcuti3_ = rcuti2_ * rcuti_;
225	<	rcuti4_ = rcuti2_ * rcuti2_;
235	>
236	>	Etypes.clear();
237	>	Etids.clear();
238	>	FQtypes.clear();
239	>	FQtids.clear();
240	>	ElectrostaticMap.clear();
241	>	Jij.clear();
242	>	nElectro_ = 0;
243	>	nFlucq_ = 0;
244
245	<	if (screeningMethod_ == DAMPED) {
246	<
229	<	alpha2_ = dampingAlpha_ * dampingAlpha_;
230	<	alpha4_ = alpha2_ * alpha2_;
231	<	alpha6_ = alpha4_ * alpha2_;
232	<	alpha8_ = alpha4_ * alpha4_;
233	<
234	<	constEXP_ = exp(-alpha2_ * cutoffRadius2_);
235	<	invRootPi_ = 0.56418958354775628695;
236	<	alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
245	>	Etids.resize( forceField_->getNAtomType(), -1);
246	>	FQtids.resize( forceField_->getNAtomType(), -1);
247
248	<	c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
249	<	c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
250	<	c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
251	<	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_;
244	<	} else {
245	<	c1c_ = rcuti_;
246	<	c2c_ = c1c_ * rcuti_;
247	<	c3c_ = 3.0 * c2c_ * rcuti_;
248	<	c4c_ = 5.0 * c3c_ * rcuti2_;
249	<	c5c_ = 7.0 * c4c_ * rcuti2_;
250	<	c6c_ = 9.0 * c5c_ * rcuti2_;
248	>	set<AtomType*>::iterator at;
249	>	for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
250	>	if ((*at)->isElectrostatic()) nElectro_++;
251	>	if ((*at)->isFluctuatingCharge()) nFlucq_++;
252		}
253	<
253	>
254	>	Jij.resize(nFlucq_);
255	>
256	>	for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
257	>	if ((*at)->isElectrostatic()) addType(*at);
258	>	}
259	>
260		if (summationMethod_ == esm_REACTION_FIELD) {
261		preRF_ = (dielectric_ - 1.0) /
262	<	((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
256	<	preRF2_ = 2.0 * preRF_;
262	>	((2.0 * dielectric_ + 1.0) * pow(cutoffRadius_,3) );
263		}
264
265	+	RealType b0c, b1c, b2c, b3c, b4c, b5c;
266	+	RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5;
267	+	RealType a2, expTerm, invArootPi(0.0);
268	+
269	+	RealType r = cutoffRadius_;
270	+	RealType r2 = r * r;
271	+	RealType ric = 1.0 / r;
272	+	RealType ric2 = ric * ric;
273	+
274	+	if (screeningMethod_ == DAMPED) {
275	+	a2 = dampingAlpha_ * dampingAlpha_;
276	+	invArootPi = 1.0 / (dampingAlpha_ * sqrt(M_PI));
277	+	expTerm = exp(-a2 * r2);
278	+	// values of Smith's B_l functions at the cutoff radius:
279	+	b0c = erfc(dampingAlpha_ * r) / r;
280	+	b1c = (      b0c     + 2.0*a2     * expTerm * invArootPi) / r2;
281	+	b2c = (3.0 * b1c + pow(2.0*a2, 2) * expTerm * invArootPi) / r2;
282	+	b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
283	+	b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
284	+	b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
285	+	// Half the Smith self piece:
286	+	selfMult1_ = - a2 * invArootPi;
287	+	selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0;
288	+	selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0;
289	+	} else {
290	+	a2 = 0.0;
291	+	b0c = 1.0 / r;
292	+	b1c = (      b0c) / r2;
293	+	b2c = (3.0 * b1c) / r2;
294	+	b3c = (5.0 * b2c) / r2;
295	+	b4c = (7.0 * b3c) / r2;
296	+	b5c = (9.0 * b4c) / r2;
297	+	selfMult1_ = 0.0;
298	+	selfMult2_ = 0.0;
299	+	selfMult4_ = 0.0;
300	+	}
301	+
302	+	// higher derivatives of B_0 at the cutoff radius:
303	+	db0c_1 = -r * b1c;
304	+	db0c_2 =     -b1c + r2 * b2c;
305	+	db0c_3 =          3.0*r*b2c  - r2*r*b3c;
306	+	db0c_4 =          3.0*b2c  - 6.0*r2*b3c     + r2*r2*b4c;
307	+	db0c_5 =                    -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
308	+
309	+	if (summationMethod_ != esm_EWALD_FULL) {
310	+	selfMult1_ -= b0c;
311	+	selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) /  3.0;
312	+	selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0;
313	+	}
314	+
315	+	// working variables for the splines:
316	+	RealType ri, ri2;
317	+	RealType b0, b1, b2, b3, b4, b5;
318	+	RealType db0_1, db0_2, db0_3, db0_4, db0_5;
319	+	RealType f, fc, f0;
320	+	RealType g, gc, g0, g1, g2, g3, g4;
321	+	RealType h, hc, h1, h2, h3, h4;
322	+	RealType s, sc, s2, s3, s4;
323	+	RealType t, tc, t3, t4;
324	+	RealType u, uc, u4;
325	+
326	+	// working variables for Taylor expansion:
327	+	RealType rmRc, rmRc2, rmRc3, rmRc4;
328	+
329	+	// Approximate using splines using a maximum of 0.1 Angstroms
330	+	// between points.
331	+	int nptest = int((cutoffRadius_ + 2.0) / 0.1);
332	+	np_ = (np_ > nptest) ? np_ : nptest;
333	+
334		// Add a 2 angstrom safety window to deal with cutoffGroups that
335		// have charged atoms longer than the cutoffRadius away from each
336	<	// other.  Splining may not be the best choice here.  Direct calls
337	<	// to erfc might be preferrable.
336	>	// other.  Splining is almost certainly the best choice here.
337	>	// Direct calls to erfc would be preferrable if it is a very fast
338	>	// implementation.
339
340	<	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
341	<	RealType rval;
342	<	vector<RealType> rvals;
343	<	vector<RealType> yvals;
344	<	for (int i = 0; i < np_; i++) {
345	<	rval = RealType(i) * dx;
346	<	rvals.push_back(rval);
347	<	yvals.push_back(erfc(dampingAlpha_ * rval));
340	>	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_);
341	>
342	>	// Storage vectors for the computed functions
343	>	vector<RealType> rv;
344	>	vector<RealType> v01v;
345	>	vector<RealType> v11v;
346	>	vector<RealType> v21v, v22v;
347	>	vector<RealType> v31v, v32v;
348	>	vector<RealType> v41v, v42v, v43v;
349	>
350	>	for (int i = 1; i < np_ + 1; i++) {
351	>	r = RealType(i) * dx;
352	>	rv.push_back(r);
353	>
354	>	ri = 1.0 / r;
355	>	ri2 = ri * ri;
356	>
357	>	r2 = r * r;
358	>	expTerm = exp(-a2 * r2);
359	>
360	>	// Taylor expansion factors (no need for factorials this way):
361	>	rmRc = r - cutoffRadius_;
362	>	rmRc2 = rmRc  * rmRc / 2.0;
363	>	rmRc3 = rmRc2 * rmRc / 3.0;
364	>	rmRc4 = rmRc3 * rmRc / 4.0;
365	>
366	>	// values of Smith's B_l functions at r:
367	>	if (screeningMethod_ == DAMPED) {
368	>	b0 = erfc(dampingAlpha_ * r) * ri;
369	>	b1 = (      b0 +     2.0*a2     * expTerm * invArootPi) * ri2;
370	>	b2 = (3.0 * b1 + pow(2.0*a2, 2) * expTerm * invArootPi) * ri2;
371	>	b3 = (5.0 * b2 + pow(2.0*a2, 3) * expTerm * invArootPi) * ri2;
372	>	b4 = (7.0 * b3 + pow(2.0*a2, 4) * expTerm * invArootPi) * ri2;
373	>	b5 = (9.0 * b4 + pow(2.0*a2, 5) * expTerm * invArootPi) * ri2;
374	>	} else {
375	>	b0 = ri;
376	>	b1 = (      b0) * ri2;
377	>	b2 = (3.0 * b1) * ri2;
378	>	b3 = (5.0 * b2) * ri2;
379	>	b4 = (7.0 * b3) * ri2;
380	>	b5 = (9.0 * b4) * ri2;
381	>	}
382	>
383	>	// higher derivatives of B_0 at r:
384	>	db0_1 = -r * b1;
385	>	db0_2 =     -b1 + r2 * b2;
386	>	db0_3 =          3.0*r*b2   - r2*r*b3;
387	>	db0_4 =          3.0*b2   - 6.0*r2*b3     + r2*r2*b4;
388	>	db0_5 =                    -15.0*r*b3 + 10.0*r2*r*b4 - r2*r2*r*b5;
389	>
390	>	f = b0;
391	>	fc = b0c;
392	>	f0 = f - fc - rmRc*db0c_1;
393	>
394	>	g = db0_1;
395	>	gc = db0c_1;
396	>	g0 = g - gc;
397	>	g1 = g0 - rmRc *db0c_2;
398	>	g2 = g1 - rmRc2*db0c_3;
399	>	g3 = g2 - rmRc3*db0c_4;
400	>	g4 = g3 - rmRc4*db0c_5;
401	>
402	>	h = db0_2;
403	>	hc = db0c_2;
404	>	h1 = h - hc;
405	>	h2 = h1 - rmRc *db0c_3;
406	>	h3 = h2 - rmRc2*db0c_4;
407	>	h4 = h3 - rmRc3*db0c_5;
408	>
409	>	s = db0_3;
410	>	sc = db0c_3;
411	>	s2 = s - sc;
412	>	s3 = s2 - rmRc *db0c_4;
413	>	s4 = s3 - rmRc2*db0c_5;
414	>
415	>	t = db0_4;
416	>	tc = db0c_4;
417	>	t3 = t - tc;
418	>	t4 = t3 - rmRc *db0c_5;
419	>
420	>	u = db0_5;
421	>	uc = db0c_5;
422	>	u4 = u - uc;
423	>
424	>	// in what follows below, the various v functions are used for
425	>	// potentials and torques, while the w functions show up in the
426	>	// forces.
427	>
428	>	switch (summationMethod_) {
429	>	case esm_SHIFTED_FORCE:
430	>
431	>	v01 = f - fc - rmRc*gc;
432	>	v11 = g - gc - rmRc*hc;
433	>	v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
434	>	v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
435	>	v31 = (h-g*ri)*ri - (hc-gc*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
436	>	v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
437	>	- rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
438	>	v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
439	>	- rmRc*(sc - 3.0*(hc-gc*ric)*ric)*ric2;
440	>	v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric
441	>	- rmRc*(tc - (4.0*sc - 9.0*(hc - gc*ric)*ric)*ric)*ric;
442	>
443	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
444	>	- (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric)
445	>	- rmRc*(uc-3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
446	>
447	>	dv01 = g - gc;
448	>	dv11 = h - hc;
449	>	dv21 = (h - g*ri)*ri - (hc - gc*ric)*ric;
450	>	dv22 = (s - (h - g*ri)*ri) - (sc - (hc - gc*ric)*ric);
451	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri - (sc - 2.0*(hc-gc*ric)*ric)*ric;
452	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri)
453	>	- (tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
454	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2 - (sc - 3.0*(hc - gc*ric)*ric)*ric2;
455	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri
456	>	- (tc - (4.0*sc - 9.0*(hc-gc*ric)*ric)*ric)*ric;
457	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri)
458	>	- (uc - 3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
459	>
460	>	break;
461	>
462	>	case esm_TAYLOR_SHIFTED:
463	>
464	>	v01 = f0;
465	>	v11 = g1;
466	>	v21 = g2 * ri;
467	>	v22 = h2 - v21;
468	>	v31 = (h3 - g3 * ri) * ri;
469	>	v32 = s3 - 3.0*v31;
470	>	v41 = (h4 - g4 * ri) * ri2;
471	>	v42 = s4 * ri - 3.0*v41;
472	>	v43 = t4 - 6.0*v42 - 3.0*v41;
473	>
474	>	dv01 = g0;
475	>	dv11 = h1;
476	>	dv21 = (h2 - g2*ri)*ri;
477	>	dv22 = (s2 - (h2 - g2*ri)*ri);
478	>	dv31 = (s3 - 2.0*(h3-g3*ri)*ri)*ri;
479	>	dv32 = (t3 - 3.0*(s3-2.0*(h3-g3*ri)*ri)*ri);
480	>	dv41 = (s4 - 3.0*(h4 - g4*ri)*ri)*ri2;
481	>	dv42 = (t4 - (4.0*s4 - 9.0*(h4-g4*ri)*ri)*ri)*ri;
482	>	dv43 = (u4 - 3.0*(2.0*t4 - (7.0*s4 - 15.0*(h4 - g4*ri)*ri)*ri)*ri);
483	>
484	>	break;
485	>
486	>	case esm_SHIFTED_POTENTIAL:
487	>
488	>	v01 = f - fc;
489	>	v11 = g - gc;
490	>	v21 = g*ri - gc*ric;
491	>	v22 = h - g*ri - (hc - gc*ric);
492	>	v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
493	>	v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
494	>	v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
495	>	v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;
496	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
497	>	- (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
498	>
499	>	dv01 = g;
500	>	dv11 = h;
501	>	dv21 = (h - g*ri)*ri;
502	>	dv22 = (s - (h - g*ri)*ri);
503	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
504	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
505	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
506	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
507	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
508	>
509	>	break;
510	>
511	>	case esm_SWITCHING_FUNCTION:
512	>	case esm_HARD:
513	>	case esm_EWALD_FULL:
514	>
515	>	v01 = f;
516	>	v11 = g;
517	>	v21 = g*ri;
518	>	v22 = h - g*ri;
519	>	v31 = (h-g*ri)*ri;
520	>	v32 = (s - 3.0*(h-g*ri)*ri);
521	>	v41 = (h - g*ri)*ri2;
522	>	v42 = (s-3.0*(h-g*ri)*ri)*ri;
523	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri);
524	>
525	>	dv01 = g;
526	>	dv11 = h;
527	>	dv21 = (h - g*ri)*ri;
528	>	dv22 = (s - (h - g*ri)*ri);
529	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
530	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
531	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
532	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
533	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
534	>
535	>	break;
536	>
537	>	case esm_REACTION_FIELD:
538	>
539	>	// following DL_POLY's lead for shifting the image charge potential:
540	>	f = b0 + preRF_ * r2;
541	>	fc = b0c + preRF_ * cutoffRadius_ * cutoffRadius_;
542	>
543	>	g = db0_1 + preRF_ * 2.0 * r;
544	>	gc = db0c_1 + preRF_ * 2.0 * cutoffRadius_;
545	>
546	>	h = db0_2 + preRF_ * 2.0;
547	>	hc = db0c_2 + preRF_ * 2.0;
548	>
549	>	v01 = f - fc;
550	>	v11 = g - gc;
551	>	v21 = g*ri - gc*ric;
552	>	v22 = h - g*ri - (hc - gc*ric);
553	>	v31 = (h-g*ri)*ri - (hc-gc*ric)*ric;
554	>	v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
555	>	v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
556	>	v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;
557	>	v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
558	>	- (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
559	>
560	>	dv01 = g;
561	>	dv11 = h;
562	>	dv21 = (h - g*ri)*ri;
563	>	dv22 = (s - (h - g*ri)*ri);
564	>	dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
565	>	dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
566	>	dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
567	>	dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
568	>	dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
569	>
570	>	break;
571	>
572	>	case esm_EWALD_PME:
573	>	case esm_EWALD_SPME:
574	>	default :
575	>	map<string, ElectrostaticSummationMethod>::iterator i;
576	>	std::string meth;
577	>	for (i = summationMap_.begin(); i != summationMap_.end(); ++i) {
578	>	if ((*i).second == summationMethod_) meth = (*i).first;
579	>	}
580	>	sprintf( painCave.errMsg,
581	>	"Electrostatic::initialize: electrostaticSummationMethod %s \n"
582	>	"\thas not been implemented yet. Please select one of:\n"
583	>	"\t\"hard\", \"shifted_potential\", or \"shifted_force\"\n",
584	>	meth.c_str() );
585	>	painCave.isFatal = 1;
586	>	simError();
587	>	break;
588	>	}
589	>
590	>	// Add these computed values to the storage vectors for spline creation:
591	>	v01v.push_back(v01);
592	>	v11v.push_back(v11);
593	>	v21v.push_back(v21);
594	>	v22v.push_back(v22);
595	>	v31v.push_back(v31);
596	>	v32v.push_back(v32);
597	>	v41v.push_back(v41);
598	>	v42v.push_back(v42);
599	>	v43v.push_back(v43);
600		}
273	–	erfcSpline_ = new CubicSpline();
274	–	erfcSpline_->addPoints(rvals, yvals);
275	–	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	+	haveElectroSplines_ = true;
625	+
626		initialized_ = true;
627		}
628
629		void Electrostatic::addType(AtomType* atomType){
630	<
630	>
631		ElectrostaticAtomData electrostaticAtomData;
632		electrostaticAtomData.is_Charge = false;
633		electrostaticAtomData.is_Dipole = false;
285	–	electrostaticAtomData.is_SplitDipole = false;
634		electrostaticAtomData.is_Quadrupole = false;
635		electrostaticAtomData.is_Fluctuating = false;
636
#	Line 297 | Line 645 | namespace OpenMD {
645		if (ma.isMultipole()) {
646		if (ma.isDipole()) {
647		electrostaticAtomData.is_Dipole = true;
648	<	electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
648	>	electrostaticAtomData.dipole = ma.getDipole();
649		}
302	–	if (ma.isSplitDipole()) {
303	–	electrostaticAtomData.is_SplitDipole = true;
304	–	electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
305	–	}
650		if (ma.isQuadrupole()) {
307	–	// Quadrupoles in OpenMD are set as the diagonal elements
308	–	// of the diagonalized traceless quadrupole moment tensor.
309	–	// The column vectors of the unitary matrix that diagonalizes
310	–	// the quadrupole moment tensor become the eFrame (or the
311	–	// electrostatic version of the body-fixed frame.
651		electrostaticAtomData.is_Quadrupole = true;
652	<	electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
652	>	electrostaticAtomData.quadrupole = ma.getQuadrupole();
653		}
654		}
655
#	Line 324 | Line 663 | namespace OpenMD {
663		electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
664		}
665
666	<	pair<map<int,AtomType*>::iterator,bool> ret;
667	<	ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
668	<	atomType) );
666	>	int atid = atomType->getIdent();
667	>	int etid = Etypes.size();
668	>	int fqtid = FQtypes.size();
669	>
670	>	pair<set<int>::iterator,bool> ret;
671	>	ret = Etypes.insert( atid );
672		if (ret.second == false) {
673		sprintf( painCave.errMsg,
674		"Electrostatic already had a previous entry with ident %d\n",
675	<	atomType->getIdent() );
675	>	atid);
676		painCave.severity = OPENMD_INFO;
677		painCave.isFatal = 0;
678		simError();
679		}
680
681	<	ElectrostaticMap[atomType] = electrostaticAtomData;
681	>	Etids[ atid ] = etid;
682	>	ElectrostaticMap.push_back(electrostaticAtomData);
683
684	<	// Now, iterate over all known types and add to the mixing map:
685	<
686	<	map<AtomType*, ElectrostaticAtomData>::iterator it;
687	<	for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
688	<	AtomType* atype2 = (*it).first;
689	<	ElectrostaticAtomData eaData2 = (*it).second;
690	<	if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
691	<
684	>	if (electrostaticAtomData.is_Fluctuating) {
685	>	ret = FQtypes.insert( atid );
686	>	if (ret.second == false) {
687	>	sprintf( painCave.errMsg,
688	>	"Electrostatic already had a previous fluctuating charge entry with ident %d\n",
689	>	atid );
690	>	painCave.severity = OPENMD_INFO;
691	>	painCave.isFatal = 0;
692	>	simError();
693	>	}
694	>	FQtids[atid] = fqtid;
695	>	Jij[fqtid].resize(nFlucq_);
696	>
697	>	// Now, iterate over all known fluctuating and add to the
698	>	// coulomb integral map:
699	>
700	>	std::set<int>::iterator it;
701	>	for( it = FQtypes.begin(); it != FQtypes.end(); ++it) {
702	>	int etid2 = Etids[ (*it) ];
703	>	int fqtid2 = FQtids[ (*it) ];
704	>	ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ];
705		RealType a = electrostaticAtomData.slaterZeta;
706		RealType b = eaData2.slaterZeta;
707		int m = electrostaticAtomData.slaterN;
708		int n = eaData2.slaterN;
709	<
709	>
710		// Create the spline of the coulombic integral for s-type
711		// Slater orbitals.  Add a 2 angstrom safety window to deal
712		// with cutoffGroups that have charged atoms longer than the
713		// cutoffRadius away from each other.
714	<
714	>
715		RealType rval;
716		RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
717		vector<RealType> rvals;
718	<	vector<RealType> J1vals;
719	<	vector<RealType> J2vals;
720	<	// don't start at i = 0, as rval = 0 is undefined for the slater overlap integrals.
718	>	vector<RealType> Jvals;
719	>	// don't start at i = 0, as rval = 0 is undefined for the
720	>	// slater overlap integrals.
721		for (int i = 1; i < np_+1; i++) {
722		rval = RealType(i) * dr;
723		rvals.push_back(rval);
724	<	J1vals.push_back(sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
725	<	// may not be necessary if Slater coulomb integral is symmetric
726	<	J2vals.push_back(sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
724	>	Jvals.push_back(sSTOCoulInt( a, b, m, n, rval *
725	>	PhysicalConstants::angstromToBohr ) *
726	>	PhysicalConstants::hartreeToKcal );
727		}
372	–
373	–	CubicSpline* J1 = new CubicSpline();
374	–	J1->addPoints(rvals, J1vals);
375	–	CubicSpline* J2 = new CubicSpline();
376	–	J2->addPoints(rvals, J2vals);
728
729	<	pair<AtomType*, AtomType*> key1, key2;
730	<	key1 = make_pair(atomType, atype2);
731	<	key2 = make_pair(atype2, atomType);
732	<
733	<	Jij[key1] = J1;
734	<	Jij[key2] = J2;
735	<	}
385	<	}
386	<
729	>	CubicSpline* J = new CubicSpline();
730	>	J->addPoints(rvals, Jvals);
731	>	Jij[fqtid][fqtid2] = J;
732	>	Jij[fqtid2].resize( nFlucq_ );
733	>	Jij[fqtid2][fqtid] = J;
734	>	}
735	>	}
736		return;
737		}
738
739		void Electrostatic::setCutoffRadius( RealType rCut ) {
740		cutoffRadius_ = rCut;
392	–	rrf_ = cutoffRadius_;
741		haveCutoffRadius_ = true;
742		}
743
396	–	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
397	–	rt_ = rSwitch;
398	–	}
744		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
745		summationMethod_ = esm;
746		}
#	Line 413 | Line 758 | namespace OpenMD {
758
759		void Electrostatic::calcForce(InteractionData &idat) {
760
416	–	// utility variables.  Should clean these up and use the Vector3d and
417	–	// Mat3x3d to replace as many as we can in future versions:
418	–
419	–	RealType q_i, q_j, mu_i, mu_j, d_i, d_j;
420	–	RealType qxx_i, qyy_i, qzz_i;
421	–	RealType qxx_j, qyy_j, qzz_j;
422	–	RealType cx_i, cy_i, cz_i;
423	–	RealType cx_j, cy_j, cz_j;
424	–	RealType cx2, cy2, cz2;
425	–	RealType ct_i, ct_j, ct_ij, a1;
426	–	RealType riji, ri, ri2, ri3, ri4;
427	–	RealType pref, vterm, epot, dudr;
428	–	RealType vpair(0.0);
429	–	RealType scale, sc2;
430	–	RealType pot_term, preVal, rfVal;
431	–	RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
432	–	RealType preSw, preSwSc;
433	–	RealType c1, c2, c3, c4;
434	–	RealType erfcVal(1.0), derfcVal(0.0);
435	–	RealType BigR;
436	–	RealType two(2.0), three(3.0);
437	–
438	–	Vector3d Q_i, Q_j;
439	–	Vector3d ux_i, uy_i, uz_i;
440	–	Vector3d ux_j, uy_j, uz_j;
441	–	Vector3d dudux_i, duduy_i, duduz_i;
442	–	Vector3d dudux_j, duduy_j, duduz_j;
443	–	Vector3d rhatdot2, rhatc4;
444	–	Vector3d dVdr;
445	–
446	–	// variables for indirect (reaction field) interactions for excluded pairs:
447	–	RealType indirect_Pot(0.0);
448	–	RealType indirect_vpair(0.0);
449	–	Vector3d indirect_dVdr(V3Zero);
450	–	Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
451	–
452	–	RealType coulInt, vFluc1(0.0), vFluc2(0.0);
453	–	pair<RealType, RealType> res;
454	–
455	–	// splines for coulomb integrals
456	–	CubicSpline* J1;
457	–	CubicSpline* J2;
458	–
761		if (!initialized_) initialize();
762
763	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
764	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
763	>	data1 = ElectrostaticMap[Etids[idat.atid1]];
764	>	data2 = ElectrostaticMap[Etids[idat.atid2]];
765	>
766	>	U = 0.0;  // Potential
767	>	F.zero();  // Force
768	>	Ta.zero(); // Torque on site a
769	>	Tb.zero(); // Torque on site b
770	>	Ea.zero(); // Electric field at site a
771	>	Eb.zero(); // Electric field at site b
772	>	Pa = 0.0;  // Site potential at site a
773	>	Pb = 0.0;  // Site potential at site b
774	>	dUdCa = 0.0; // fluctuating charge force at site a
775	>	dUdCb = 0.0; // fluctuating charge force at site a
776
777	<	// some variables we'll need independent of electrostatic type:
777	>	// Indirect interactions mediated by the reaction field.
778	>	indirect_Pot = 0.0;   // Potential
779	>	indirect_F.zero();    // Force
780	>	indirect_Ta.zero();   // Torque on site a
781	>	indirect_Tb.zero();   // Torque on site b
782
783	<	riji = 1.0 /  *(idat.rij) ;
784	<	Vector3d rhat =  *(idat.d)   * riji;
783	>	// Excluded potential that is still computed for fluctuating charges
784	>	excluded_Pot= 0.0;
785
786	+	// some variables we'll need independent of electrostatic type:
787	+
788	+	ri = 1.0 /  *(idat.rij);
789	+	rhat =  *(idat.d)  * ri;
790	+
791		// logicals
792
793	<	bool i_is_Charge = data1.is_Charge;
794	<	bool i_is_Dipole = data1.is_Dipole;
795	<	bool i_is_SplitDipole = data1.is_SplitDipole;
796	<	bool i_is_Quadrupole = data1.is_Quadrupole;
475	<	bool i_is_Fluctuating = data1.is_Fluctuating;
793	>	a_is_Charge = data1.is_Charge;
794	>	a_is_Dipole = data1.is_Dipole;
795	>	a_is_Quadrupole = data1.is_Quadrupole;
796	>	a_is_Fluctuating = data1.is_Fluctuating;
797
798	<	bool j_is_Charge = data2.is_Charge;
799	<	bool j_is_Dipole = data2.is_Dipole;
800	<	bool j_is_SplitDipole = data2.is_SplitDipole;
801	<	bool j_is_Quadrupole = data2.is_Quadrupole;
802	<	bool j_is_Fluctuating = data2.is_Fluctuating;
798	>	b_is_Charge = data2.is_Charge;
799	>	b_is_Dipole = data2.is_Dipole;
800	>	b_is_Quadrupole = data2.is_Quadrupole;
801	>	b_is_Fluctuating = data2.is_Fluctuating;
802	>
803	>	// Obtain all of the required radial function values from the
804	>	// spline structures:
805
806	<	if (i_is_Charge) {
807	<	q_i = data1.fixedCharge;
806	>	// needed for fields (and forces):
807	>	if (a_is_Charge || b_is_Charge) {
808	>	v01s->getValueAndDerivativeAt( *(idat.rij), v01, dv01);
809	>	}
810	>	if (a_is_Dipole || b_is_Dipole) {
811	>	v11s->getValueAndDerivativeAt( *(idat.rij), v11, dv11);
812	>	v11or = ri * v11;
813	>	}
814	>	if (a_is_Quadrupole || b_is_Quadrupole ||  (a_is_Dipole && b_is_Dipole)) {
815	>	v21s->getValueAndDerivativeAt( *(idat.rij), v21, dv21);
816	>	v22s->getValueAndDerivativeAt( *(idat.rij), v22, dv22);
817	>	v22or = ri * v22;
818	>	}
819
820	<	if (i_is_Fluctuating) {
821	<	q_i += *(idat.flucQ1);
820	>	// needed for potentials (and forces and torques):
821	>	if ((a_is_Dipole && b_is_Quadrupole) ||
822	>	(b_is_Dipole && a_is_Quadrupole)) {
823	>	v31s->getValueAndDerivativeAt( *(idat.rij), v31, dv31);
824	>	v32s->getValueAndDerivativeAt( *(idat.rij), v32, dv32);
825	>	v31or = v31 * ri;
826	>	v32or = v32 * ri;
827	>	}
828	>	if (a_is_Quadrupole && b_is_Quadrupole) {
829	>	v41s->getValueAndDerivativeAt( *(idat.rij), v41, dv41);
830	>	v42s->getValueAndDerivativeAt( *(idat.rij), v42, dv42);
831	>	v43s->getValueAndDerivativeAt( *(idat.rij), v43, dv43);
832	>	v42or = v42 * ri;
833	>	v43or = v43 * ri;
834	>	}
835	>
836	>	// calculate the single-site contributions (fields, etc).
837	>
838	>	if (a_is_Charge) {
839	>	C_a = data1.fixedCharge;
840	>
841	>	if (a_is_Fluctuating) {
842	>	C_a += *(idat.flucQ1);
843		}
844
845		if (idat.excluded) {
846	<	*(idat.skippedCharge2) += q_i;
846	>	*(idat.skippedCharge2) += C_a;
847	>	} else {
848	>	// only do the field and site potentials if we're not excluded:
849	>	Eb -= C_a *  pre11_ * dv01 * rhat;
850	>	Pb += C_a *  pre11_ * v01;
851		}
852		}
494	–
495	–	if (i_is_Dipole) {
496	–	mu_i = data1.dipole_moment;
497	–	uz_i = idat.eFrame1->getColumn(2);
498	–
499	–	ct_i = dot(uz_i, rhat);
500	–
501	–	if (i_is_SplitDipole)
502	–	d_i = data1.split_dipole_distance;
503	–
504	–	duduz_i = V3Zero;
505	–	}
853
854	<	if (i_is_Quadrupole) {
855	<	Q_i = data1.quadrupole_moments;
856	<	qxx_i = Q_i.x();
857	<	qyy_i = Q_i.y();
858	<	qzz_i = Q_i.z();
859	<
860	<	ux_i = idat.eFrame1->getColumn(0);
861	<	uy_i = idat.eFrame1->getColumn(1);
515	<	uz_i = idat.eFrame1->getColumn(2);
854	>	if (a_is_Dipole) {
855	>	D_a = *(idat.dipole1);
856	>	rdDa = dot(rhat, D_a);
857	>	rxDa = cross(rhat, D_a);
858	>	if (!idat.excluded) {
859	>	Eb -=  pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
860	>	Pb +=  pre12_ * v11 * rdDa;
861	>	}
862
517	–	cx_i = dot(ux_i, rhat);
518	–	cy_i = dot(uy_i, rhat);
519	–	cz_i = dot(uz_i, rhat);
520	–
521	–	dudux_i = V3Zero;
522	–	duduy_i = V3Zero;
523	–	duduz_i = V3Zero;
863		}
864	<
865	<	if (j_is_Charge) {
866	<	q_j = data2.fixedCharge;
867	<
868	<	if (j_is_Fluctuating)
869	<	q_j += *(idat.flucQ2);
870	<
871	<	if (idat.excluded) {
872	<	*(idat.skippedCharge1) += q_j;
864	>
865	>	if (a_is_Quadrupole) {
866	>	Q_a = *(idat.quadrupole1);
867	>	trQa =  Q_a.trace();
868	>	Qar =   Q_a * rhat;
869	>	rQa = rhat * Q_a;
870	>	rdQar = dot(rhat, Qar);
871	>	rxQar = cross(rhat, Qar);
872	>	if (!idat.excluded) {
873	>	Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
874	>	+ rdQar * rhat * (dv22 - 2.0*v22or));
875	>	Pb += pre14_ * (v21 * trQa + v22 * rdQar);
876		}
877		}
878	<
879	<
880	<	if (j_is_Dipole) {
539	<	mu_j = data2.dipole_moment;
540	<	uz_j = idat.eFrame2->getColumn(2);
878	>
879	>	if (b_is_Charge) {
880	>	C_b = data2.fixedCharge;
881
882	<	ct_j = dot(uz_j, rhat);
883	<
544	<	if (j_is_SplitDipole)
545	<	d_j = data2.split_dipole_distance;
882	>	if (b_is_Fluctuating)
883	>	C_b += *(idat.flucQ2);
884
885	<	duduz_j = V3Zero;
885	>	if (idat.excluded) {
886	>	*(idat.skippedCharge1) += C_b;
887	>	} else {
888	>	// only do the field if we're not excluded:
889	>	Ea += C_b *  pre11_ * dv01 * rhat;
890	>	Pa += C_b *  pre11_ * v01;
891	>
892	>	}
893		}
894
895	<	if (j_is_Quadrupole) {
896	<	Q_j = data2.quadrupole_moments;
897	<	qxx_j = Q_j.x();
898	<	qyy_j = Q_j.y();
899	<	qzz_j = Q_j.z();
900	<
901	<	ux_j = idat.eFrame2->getColumn(0);
902	<	uy_j = idat.eFrame2->getColumn(1);
558	<	uz_j = idat.eFrame2->getColumn(2);
559	<
560	<	cx_j = dot(ux_j, rhat);
561	<	cy_j = dot(uy_j, rhat);
562	<	cz_j = dot(uz_j, rhat);
563	<
564	<	dudux_j = V3Zero;
565	<	duduy_j = V3Zero;
566	<	duduz_j = V3Zero;
895	>	if (b_is_Dipole) {
896	>	D_b = *(idat.dipole2);
897	>	rdDb = dot(rhat, D_b);
898	>	rxDb = cross(rhat, D_b);
899	>	if (!idat.excluded) {
900	>	Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
901	>	Pa += pre12_ * v11 * rdDb;
902	>	}
903		}
904
905	<	if (i_is_Fluctuating && j_is_Fluctuating) {
906	<	J1 = Jij[idat.atypes];
907	<	J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
905	>	if (b_is_Quadrupole) {
906	>	Q_b = *(idat.quadrupole2);
907	>	trQb =  Q_b.trace();
908	>	Qbr =   Q_b * rhat;
909	>	rQb = rhat * Q_b;
910	>	rdQbr = dot(rhat, Qbr);
911	>	rxQbr = cross(rhat, Qbr);
912	>	if (!idat.excluded) {
913	>	Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
914	>	+ rdQbr * rhat * (dv22 - 2.0*v22or));
915	>	Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
916	>	}
917		}
918	+
919
920	<	epot = 0.0;
921	<	dVdr = V3Zero;
922	<
923	<	if (i_is_Charge) {
920	>	if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
921	>	J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
922	>	}
923	>
924	>	if (a_is_Charge) {
925
926	<	if (j_is_Charge) {
927	<	if (screeningMethod_ == DAMPED) {
928	<	// assemble the damping variables
929	<	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
583	<	erfcVal = res.first;
584	<	derfcVal = res.second;
926	>	if (b_is_Charge) {
927	>	pref =  pre11_ * *(idat.electroMult);
928	>	U  += C_a * C_b * pref * v01;
929	>	F  += C_a * C_b * pref * dv01 * rhat;
930
931	<	//erfcVal = erfc(dampingAlpha_ * *(idat.rij));
932	<	//derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
931	>	// If this is an excluded pair, there are still indirect
932	>	// interactions via the reaction field we must worry about:
933
934	<	c1 = erfcVal * riji;
935	<	c2 = (-derfcVal + c1) * riji;
936	<	} else {
937	<	c1 = riji;
593	<	c2 = c1 * riji;
934	>	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
935	>	rfContrib = preRF_ * pref * C_a * C_b * *(idat.r2);
936	>	indirect_Pot += rfContrib;
937	>	indirect_F   += rfContrib * 2.0 * ri * rhat;
938		}
939
940	<	preVal =  *(idat.electroMult) * pre11_;
941	<
942	<	if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
599	<	vterm = preVal * (c1 - c1c_);
600	<	dudr  = - *(idat.sw)  * preVal * c2;
940	>	// Fluctuating charge forces are handled via Coulomb integrals
941	>	// for excluded pairs (i.e. those connected via bonds) and
942	>	// with the standard charge-charge interaction otherwise.
943
944	<	} else if (summationMethod_ == esm_SHIFTED_FORCE)  {
945	<	vterm = preVal * ( c1 - c1c_ + c2c_*( *(idat.rij)  - cutoffRadius_) );
946	<	dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
947	<
948	<	} else if (summationMethod_ == esm_REACTION_FIELD) {
607	<	rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
608	<
609	<	vterm = preVal * ( riji + rfVal );
610	<	dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
611	<
612	<	// if this is an excluded pair, there are still indirect
613	<	// interactions via the reaction field we must worry about:
614	<
615	<	if (idat.excluded) {
616	<	indirect_vpair += preVal * rfVal;
617	<	indirect_Pot += *(idat.sw) * preVal * rfVal;
618	<	indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
944	>	if (idat.excluded) {
945	>	if (a_is_Fluctuating || b_is_Fluctuating) {
946	>	coulInt = J->getValueAt( *(idat.rij) );
947	>	if (a_is_Fluctuating) dUdCa += C_b * coulInt;
948	>	if (b_is_Fluctuating) dUdCb += C_a * coulInt;
949		}
620	–
950		} else {
951	+	if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
952	+	if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
953	+	}
954	+	}
955
956	<	vterm = preVal * riji * erfcVal;
957	<	dudr  = -  *(idat.sw)  * preVal * c2;
958	<
959	<	}
960	<
628	<	vpair += vterm * q_i * q_j;
629	<	epot +=  *(idat.sw)  * vterm * q_i * q_j;
630	<	dVdr += dudr * rhat * q_i * q_j;
956	>	if (b_is_Dipole) {
957	>	pref =  pre12_ * *(idat.electroMult);
958	>	U  += C_a * pref * v11 * rdDb;
959	>	F  += C_a * pref * ((dv11 - v11or) * rdDb * rhat + v11or * D_b);
960	>	Tb += C_a * pref * v11 * rxDb;
961
962	<	if (i_is_Fluctuating) {
633	<	if (idat.excluded) {
634	<	// vFluc1 is the difference between the direct coulomb integral
635	<	// and the normal 1/r-like  interaction between point charges.
636	<	coulInt = J1->getValueAt( *(idat.rij) );
637	<	vFluc1 = coulInt - (*(idat.sw) * vterm);
638	<	} else {
639	<	vFluc1 = 0.0;
640	<	}
641	<	*(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j;
642	<	}
962	>	if (a_is_Fluctuating) dUdCa += pref * v11 * rdDb;
963
964	<	if (j_is_Fluctuating) {
965	<	if (idat.excluded) {
966	<	// vFluc2 is the difference between the direct coulomb integral
967	<	// and the normal 1/r-like  interaction between point charges.
968	<	coulInt = J2->getValueAt( *(idat.rij) );
969	<	vFluc2 = coulInt - (*(idat.sw) * vterm);
970	<	} else {
971	<	vFluc2 = 0.0;
972	<	}
653	<	*(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i;
964	>	// Even if we excluded this pair from direct interactions, we
965	>	// still have the reaction-field-mediated charge-dipole
966	>	// interaction:
967	>
968	>	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
969	>	rfContrib = C_a * pref * preRF_ * 2.0 * *(idat.rij);
970	>	indirect_Pot += rfContrib * rdDb;
971	>	indirect_F   += rfContrib * D_b / (*idat.rij);
972	>	indirect_Tb  += C_a * pref * preRF_ * rxDb;
973		}
655	–
656	–
974		}
975
976	<	if (j_is_Dipole) {
977	<	// pref is used by all the possible methods
978	<	pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
979	<	preSw =  *(idat.sw)  * pref;
976	>	if (b_is_Quadrupole) {
977	>	pref = pre14_ * *(idat.electroMult);
978	>	U  +=  C_a * pref * (v21 * trQb + v22 * rdQbr);
979	>	F  +=  C_a * pref * (trQb * dv21 * rhat + 2.0 * Qbr * v22or);
980	>	F  +=  C_a * pref * rdQbr * rhat * (dv22 - 2.0*v22or);
981	>	Tb +=  C_a * pref * 2.0 * rxQbr * v22;
982
983	<	if (summationMethod_ == esm_REACTION_FIELD) {
984	<	ri2 = riji * riji;
985	<	ri3 = ri2 * riji;
667	<
668	<	vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
669	<	vpair += vterm;
670	<	epot +=  *(idat.sw)  * vterm;
983	>	if (a_is_Fluctuating) dUdCa += pref * (v21 * trQb + v22 * rdQbr);
984	>	}
985	>	}
986
987	<	dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
673	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
987	>	if (a_is_Dipole) {
988
989	<	// Even if we excluded this pair from direct interactions,
990	<	// we still have the reaction-field-mediated charge-dipole
677	<	// interaction:
989	>	if (b_is_Charge) {
990	>	pref = pre12_ * *(idat.electroMult);
991
992	<	if (idat.excluded) {
993	<	indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
994	<	indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
682	<	indirect_dVdr += preSw * preRF2_ * uz_j;
683	<	indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
684	<	}
685	<
686	<	} else {
687	<	// determine the inverse r used if we have split dipoles
688	<	if (j_is_SplitDipole) {
689	<	BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
690	<	ri = 1.0 / BigR;
691	<	scale =  *(idat.rij)  * ri;
692	<	} else {
693	<	ri = riji;
694	<	scale = 1.0;
695	<	}
696	<
697	<	sc2 = scale * scale;
992	>	U  -= C_b * pref * v11 * rdDa;
993	>	F  -= C_b * pref * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
994	>	Ta -= C_b * pref * v11 * rxDa;
995
996	<	if (screeningMethod_ == DAMPED) {
700	<	// assemble the damping variables
701	<	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
702	<	//erfcVal = res.first;
703	<	//derfcVal = res.second;
704	<	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
705	<	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
706	<	c1 = erfcVal * ri;
707	<	c2 = (-derfcVal + c1) * ri;
708	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
709	<	} else {
710	<	c1 = ri;
711	<	c2 = c1 * ri;
712	<	c3 = 3.0 * c2 * ri;
713	<	}
714	<
715	<	c2ri = c2 * ri;
996	>	if (b_is_Fluctuating) dUdCb -= pref * v11 * rdDa;
997
998	<	// calculate the potential
999	<	pot_term =  scale * c2;
1000	<	vterm = -pref * ct_j * pot_term;
1001	<	vpair += vterm;
1002	<	epot +=  *(idat.sw)  * vterm;
1003	<
1004	<	// calculate derivatives for forces and torques
1005	<
725	<	dVdr += -preSw * (uz_j * c2ri - ct_j * rhat * sc2 * c3);
726	<	duduz_j += -preSw * pot_term * rhat;
727	<
998	>	// Even if we excluded this pair from direct interactions,
999	>	// we still have the reaction-field-mediated charge-dipole
1000	>	// interaction:
1001	>	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
1002	>	rfContrib = C_b * pref * preRF_ * 2.0 * *(idat.rij);
1003	>	indirect_Pot -= rfContrib * rdDa;
1004	>	indirect_F   -= rfContrib * D_a / (*idat.rij);
1005	>	indirect_Ta  -= C_b * pref * preRF_ * rxDa;
1006		}
729	–	if (i_is_Fluctuating) {
730	–	*(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
731	–	}
1007		}
1008
1009	<	if (j_is_Quadrupole) {
1010	<	// first precalculate some necessary variables
1011	<	cx2 = cx_j * cx_j;
1012	<	cy2 = cy_j * cy_j;
738	<	cz2 = cz_j * cz_j;
739	<	pref =   *(idat.electroMult) * pre14_ * q_i * one_third_;
740	<
741	<	if (screeningMethod_ == DAMPED) {
742	<	// assemble the damping variables
743	<	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
744	<	//erfcVal = res.first;
745	<	//derfcVal = res.second;
746	<	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
747	<	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
748	<	c1 = erfcVal * riji;
749	<	c2 = (-derfcVal + c1) * riji;
750	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
751	<	c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * riji * riji;
752	<	} else {
753	<	c1 = riji;
754	<	c2 = c1 * riji;
755	<	c3 = 3.0 * c2 * riji;
756	<	c4 = 5.0 * c3 * riji * riji;
757	<	}
1009	>	if (b_is_Dipole) {
1010	>	pref = pre22_ * *(idat.electroMult);
1011	>	DadDb = dot(D_a, D_b);
1012	>	DaxDb = cross(D_a, D_b);
1013
1014	<	// precompute variables for convenience
1015	<	preSw =  *(idat.sw)  * pref;
1016	<	c2ri = c2 * riji;
1017	<	c3ri = c3 * riji;
1018	<	c4rij = c4 *  *(idat.rij) ;
1019	<	rhatdot2 = two * rhat * c3;
1020	<	rhatc4 = rhat * c4rij;
1021	<
1022	<	// calculate the potential
1023	<	pot_term = ( qxx_j * (cx2*c3 - c2ri) +
1024	<	qyy_j * (cy2*c3 - c2ri) +
1025	<	qzz_j * (cz2*c3 - c2ri) );
1026	<	vterm = pref * pot_term;
772	<	vpair += vterm;
773	<	epot +=  *(idat.sw)  * vterm;
774	<
775	<	// calculate derivatives for the forces and torques
776	<
777	<	dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
778	<	qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
779	<	qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
780	<
781	<	dudux_j += preSw * qxx_j * cx_j * rhatdot2;
782	<	duduy_j += preSw * qyy_j * cy_j * rhatdot2;
783	<	duduz_j += preSw * qzz_j * cz_j * rhatdot2;
784	<	if (i_is_Fluctuating) {
785	<	*(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
1014	>	U  -= pref * (DadDb * v21 + rdDa * rdDb * v22);
1015	>	F  -= pref * (dv21 * DadDb * rhat + v22or * (rdDb * D_a + rdDa * D_b));
1016	>	F  -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
1017	>	Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
1018	>	Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
1019	>	// Even if we excluded this pair from direct interactions, we
1020	>	// still have the reaction-field-mediated dipole-dipole
1021	>	// interaction:
1022	>	if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
1023	>	rfContrib = -pref * preRF_ * 2.0;
1024	>	indirect_Pot += rfContrib * DadDb;
1025	>	indirect_Ta  += rfContrib * DaxDb;
1026	>	indirect_Tb  -= rfContrib * DaxDb;
1027		}
1028	+	}
1029
1030	+	if (b_is_Quadrupole) {
1031	+	pref = pre24_ * *(idat.electroMult);
1032	+	DadQb = D_a * Q_b;
1033	+	DadQbr = dot(D_a, Qbr);
1034	+	DaxQbr = cross(D_a, Qbr);
1035	+
1036	+	U  -= pref * ((trQb*rdDa + 2.0*DadQbr)*v31 + rdDa*rdQbr*v32);
1037	+	F  -= pref * (trQb*D_a + 2.0*DadQb) * v31or;
1038	+	F  -= pref * (trQb*rdDa + 2.0*DadQbr) * (dv31-v31or) * rhat;
1039	+	F  -= pref * (D_a*rdQbr + 2.0*rdDa*rQb) * v32or;
1040	+	F  -= pref * (rdDa * rdQbr * rhat * (dv32-3.0*v32or));
1041	+	Ta += pref * ((-trQb*rxDa + 2.0 * DaxQbr)*v31 - rxDa*rdQbr*v32);
1042	+	Tb += pref * ((2.0*cross(DadQb, rhat) - 2.0*DaxQbr)*v31
1043	+	- 2.0*rdDa*rxQbr*v32);
1044		}
1045		}
790	–
791	–	if (i_is_Dipole) {
1046
1047	<	if (j_is_Charge) {
1048	<	// variables used by all the methods
1049	<	pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1050	<	preSw =  *(idat.sw)  * pref;
1047	>	if (a_is_Quadrupole) {
1048	>	if (b_is_Charge) {
1049	>	pref = pre14_ * *(idat.electroMult);
1050	>	U  += C_b * pref * (v21 * trQa + v22 * rdQar);
1051	>	F  += C_b * pref * (trQa * rhat * dv21 + 2.0 * Qar * v22or);
1052	>	F  += C_b * pref * rdQar * rhat * (dv22 - 2.0*v22or);
1053	>	Ta += C_b * pref * 2.0 * rxQar * v22;
1054
1055	<	if (summationMethod_ == esm_REACTION_FIELD) {
1055	>	if (b_is_Fluctuating) dUdCb += pref * (v21 * trQa + v22 * rdQar);
1056	>	}
1057	>	if (b_is_Dipole) {
1058	>	pref = pre24_ * *(idat.electroMult);
1059	>	DbdQa = D_b * Q_a;
1060	>	DbdQar = dot(D_b, Qar);
1061	>	DbxQar = cross(D_b, Qar);
1062
1063	<	ri2 = riji * riji;
1064	<	ri3 = ri2 * riji;
1063	>	U  += pref * ((trQa*rdDb + 2.0*DbdQar)*v31 + rdDb*rdQar*v32);
1064	>	F  += pref * (trQa*D_b + 2.0*DbdQa) * v31or;
1065	>	F  += pref * (trQa*rdDb + 2.0*DbdQar) * (dv31-v31or) * rhat;
1066	>	F  += pref * (D_b*rdQar + 2.0*rdDb*rQa) * v32or;
1067	>	F  += pref * (rdDb * rdQar * rhat * (dv32-3.0*v32or));
1068	>	Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31
1069	>	+ 2.0*rdDb*rxQar*v32);
1070	>	Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1071	>	}
1072	>	if (b_is_Quadrupole) {
1073	>	pref = pre44_ * *(idat.electroMult);  // yes
1074	>	QaQb = Q_a * Q_b;
1075	>	trQaQb = QaQb.trace();
1076	>	rQaQb = rhat * QaQb;
1077	>	QaQbr = QaQb * rhat;
1078	>	QaxQb = mCross(Q_a, Q_b);
1079	>	rQaQbr = dot(rQa, Qbr);
1080	>	rQaxQbr = cross(rQa, Qbr);
1081	>
1082	>	U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1083	>	U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
1084	>	U  += pref * (rdQar * rdQbr) * v43;
1085
1086	<	vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1087	<	vpair += vterm;
1088	<	epot +=  *(idat.sw)  * vterm;
806	<
807	<	dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
808	<
809	<	duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1086	>	F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*dv41;
1087	>	F  += pref*rhat*(trQa*rdQbr + trQb*rdQar + 4.0*rQaQbr)*(dv42-2.0*v42or);
1088	>	F  += pref * rhat * (rdQar * rdQbr)*(dv43 - 4.0*v43or);
1089
1090	<	// Even if we excluded this pair from direct interactions,
1091	<	// we still have the reaction-field-mediated charge-dipole
1092	<	// interaction:
1090	>	F  += pref * 2.0 * (trQb*rQa + trQa*rQb) * v42or;
1091	>	F  += pref * 4.0 * (rQaQb + QaQbr) * v42or;
1092	>	F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb) * v43or;
1093
1094	<	if (idat.excluded) {
1095	<	indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
1096	<	indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
1097	<	indirect_dVdr += -preSw * preRF2_ * uz_i;
1098	<	indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
820	<	}
821	<
822	<	} else {
823	<
824	<	// determine inverse r if we are using split dipoles
825	<	if (i_is_SplitDipole) {
826	<	BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
827	<	ri = 1.0 / BigR;
828	<	scale =  *(idat.rij)  * ri;
829	<	} else {
830	<	ri = riji;
831	<	scale = 1.0;
832	<	}
833	<
834	<	sc2 = scale * scale;
835	<
836	<	if (screeningMethod_ == DAMPED) {
837	<	// assemble the damping variables
838	<	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
839	<	//erfcVal = res.first;
840	<	//derfcVal = res.second;
841	<	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
842	<	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
843	<	c1 = erfcVal * ri;
844	<	c2 = (-derfcVal + c1) * ri;
845	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
846	<	} else {
847	<	c1 = ri;
848	<	c2 = c1 * ri;
849	<	c3 = 3.0 * c2 * ri;
850	<	}
851	<
852	<	c2ri = c2 * ri;
853	<
854	<	// calculate the potential
855	<	pot_term = c2 * scale;
856	<	vterm = pref * ct_i * pot_term;
857	<	vpair += vterm;
858	<	epot +=  *(idat.sw)  * vterm;
1094	>	Ta += pref * (- 4.0 * QaxQb * v41);
1095	>	Ta += pref * (- 2.0 * trQb * cross(rQa, rhat)
1096	>	+ 4.0 * cross(rhat, QaQbr)
1097	>	- 4.0 * rQaxQbr) * v42;
1098	>	Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
1099
860	–	// calculate derivatives for the forces and torques
861	–	dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
862	–	duduz_i += preSw * pot_term * rhat;
863	–	}
1100
1101	<	if (j_is_Fluctuating) {
1102	<	*(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1103	<	}
1101	>	Tb += pref * (+ 4.0 * QaxQb * v41);
1102	>	Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1103	>	- 4.0 * cross(rQaQb, rhat)
1104	>	+ 4.0 * rQaxQbr) * v42;
1105	>	// Possible replacement for line 2 above:
1106	>	//             + 4.0 * cross(rhat, QbQar)
1107
1108	+	Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1109		}
1110	<
871	<	if (j_is_Dipole) {
872	<	// variables used by all methods
873	<	ct_ij = dot(uz_i, uz_j);
874	<
875	<	pref =  *(idat.electroMult) * pre22_ * mu_i * mu_j;
876	<	preSw =  *(idat.sw)  * pref;
877	<
878	<	if (summationMethod_ == esm_REACTION_FIELD) {
879	<	ri2 = riji * riji;
880	<	ri3 = ri2 * riji;
881	<	ri4 = ri2 * ri2;
1110	>	}
1111
1112	<	vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
1113	<	preRF2_ * ct_ij );
1114	<	vpair += vterm;
886	<	epot +=  *(idat.sw)  * vterm;
887	<
888	<	a1 = 5.0 * ct_i * ct_j - ct_ij;
889	<
890	<	dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
891	<
892	<	duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
893	<	duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
894	<
895	<	if (idat.excluded) {
896	<	indirect_vpair +=  - pref * preRF2_ * ct_ij;
897	<	indirect_Pot +=    - preSw * preRF2_ * ct_ij;
898	<	indirect_duduz_i += -preSw * preRF2_ * uz_j;
899	<	indirect_duduz_j += -preSw * preRF2_ * uz_i;
900	<	}
901	<
902	<	} else {
903	<
904	<	if (i_is_SplitDipole) {
905	<	if (j_is_SplitDipole) {
906	<	BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
907	<	} else {
908	<	BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
909	<	}
910	<	ri = 1.0 / BigR;
911	<	scale =  *(idat.rij)  * ri;
912	<	} else {
913	<	if (j_is_SplitDipole) {
914	<	BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
915	<	ri = 1.0 / BigR;
916	<	scale =  *(idat.rij)  * ri;
917	<	} else {
918	<	ri = riji;
919	<	scale = 1.0;
920	<	}
921	<	}
922	<	if (screeningMethod_ == DAMPED) {
923	<	// assemble damping variables
924	<	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
925	<	//erfcVal = res.first;
926	<	//derfcVal = res.second;
927	<	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
928	<	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
929	<	c1 = erfcVal * ri;
930	<	c2 = (-derfcVal + c1) * ri;
931	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
932	<	c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * ri * ri;
933	<	} else {
934	<	c1 = ri;
935	<	c2 = c1 * ri;
936	<	c3 = 3.0 * c2 * ri;
937	<	c4 = 5.0 * c3 * ri * ri;
938	<	}
939	<
940	<	// precompute variables for convenience
941	<	sc2 = scale * scale;
942	<	cti3 = ct_i * sc2 * c3;
943	<	ctj3 = ct_j * sc2 * c3;
944	<	ctidotj = ct_i * ct_j * sc2;
945	<	preSwSc = preSw * scale;
946	<	c2ri = c2 * ri;
947	<	c3ri = c3 * ri;
948	<	c4rij = c4 *  *(idat.rij) ;
949	<
950	<	// calculate the potential
951	<	pot_term = (ct_ij * c2ri - ctidotj * c3);
952	<	vterm = pref * pot_term;
953	<	vpair += vterm;
954	<	epot +=  *(idat.sw)  * vterm;
955	<
956	<	// calculate derivatives for the forces and torques
957	<	dVdr += preSwSc * ( ctidotj * rhat * c4rij  -
958	<	(ct_i*uz_j + ct_j*uz_i + ct_ij*rhat) * c3ri);
959	<
960	<	duduz_i += preSw * (uz_j * c2ri - ctj3 * rhat);
961	<	duduz_j += preSw * (uz_i * c2ri - cti3 * rhat);
962	<	}
963	<	}
1112	>	if (idat.doElectricField) {
1113	>	*(idat.eField1) += Ea * *(idat.electroMult);
1114	>	*(idat.eField2) += Eb * *(idat.electroMult);
1115		}
1116
1117	<	if (i_is_Quadrupole) {
1118	<	if (j_is_Charge) {
1119	<	// precompute some necessary variables
969	<	cx2 = cx_i * cx_i;
970	<	cy2 = cy_i * cy_i;
971	<	cz2 = cz_i * cz_i;
972	<
973	<	pref =  *(idat.electroMult) * pre14_ * q_j * one_third_;
974	<
975	<	if (screeningMethod_ == DAMPED) {
976	<	// assemble the damping variables
977	<	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
978	<	//erfcVal = res.first;
979	<	//derfcVal = res.second;
980	<	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
981	<	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
982	<	c1 = erfcVal * riji;
983	<	c2 = (-derfcVal + c1) * riji;
984	<	c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
985	<	c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * riji * riji;
986	<	} else {
987	<	c1 = riji;
988	<	c2 = c1 * riji;
989	<	c3 = 3.0 * c2 * riji;
990	<	c4 = 5.0 * c3 * riji * riji;
991	<	}
992	<
993	<	// precompute some variables for convenience
994	<	preSw =  *(idat.sw)  * pref;
995	<	c2ri = c2 * riji;
996	<	c3ri = c3 * riji;
997	<	c4rij = c4 *  *(idat.rij) ;
998	<	rhatdot2 = two * rhat * c3;
999	<	rhatc4 = rhat * c4rij;
1000	<
1001	<	// calculate the potential
1002	<	pot_term = ( qxx_i * (cx2 * c3 - c2ri) +
1003	<	qyy_i * (cy2 * c3 - c2ri) +
1004	<	qzz_i * (cz2 * c3 - c2ri) );
1005	<
1006	<	vterm = pref * pot_term;
1007	<	vpair += vterm;
1008	<	epot +=  *(idat.sw)  * vterm;
1009	<
1010	<	// calculate the derivatives for the forces and torques
1011	<
1012	<	dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
1013	<	qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
1014	<	qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
1015	<
1016	<	dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1017	<	duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1018	<	duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1019	<
1020	<	if (j_is_Fluctuating) {
1021	<	*(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1022	<	}
1023	<
1024	<	}
1117	>	if (idat.doSitePotential) {
1118	>	*(idat.sPot1) += Pa * *(idat.electroMult);
1119	>	*(idat.sPot2) += Pb * *(idat.electroMult);
1120		}
1121
1122	+	if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
1123	+	if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw);
1124
1125		if (!idat.excluded) {
1029	–	*(idat.vpair) += vpair;
1030	–	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
1031	–	*(idat.f1) += dVdr;
1126
1127	<	if (i_is_Dipole || i_is_Quadrupole)
1128	<	*(idat.t1) -= cross(uz_i, duduz_i);
1129	<	if (i_is_Quadrupole) {
1036	<	*(idat.t1) -= cross(ux_i, dudux_i);
1037	<	*(idat.t1) -= cross(uy_i, duduy_i);
1038	<	}
1127	>	*(idat.vpair) += U;
1128	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += U * *(idat.sw);
1129	>	*(idat.f1) += F * *(idat.sw);
1130
1131	<	if (j_is_Dipole || j_is_Quadrupole)
1132	<	*(idat.t2) -= cross(uz_j, duduz_j);
1042	<	if (j_is_Quadrupole) {
1043	<	*(idat.t2) -= cross(uz_j, dudux_j);
1044	<	*(idat.t2) -= cross(uz_j, duduy_j);
1045	<	}
1131	>	if (a_is_Dipole || a_is_Quadrupole)
1132	>	*(idat.t1) += Ta * *(idat.sw);
1133
1134	+	if (b_is_Dipole || b_is_Quadrupole)
1135	+	*(idat.t2) += Tb * *(idat.sw);
1136	+
1137		} else {
1138
1139		// only accumulate the forces and torques resulting from the
1140		// indirect reaction field terms.
1141
1142	<	*(idat.vpair) += indirect_vpair;
1142	>	*(idat.vpair) += indirect_Pot;
1143	>	(*(idat.excludedPot))[ELECTROSTATIC_FAMILY] +=  excluded_Pot;
1144	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += *(idat.sw) * indirect_Pot;
1145	>	*(idat.f1) += *(idat.sw) * indirect_F;
1146
1147	<	(*(idat.excludedPot))[ELECTROSTATIC_FAMILY] +=   (*(idat.sw) * vterm +
1148	<	vFluc1 ) * q_i * q_j;
1149	<	(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1150	<	*(idat.f1) += indirect_dVdr;
1151	<
1059	<	if (i_is_Dipole)
1060	<	*(idat.t1) -= cross(uz_i, indirect_duduz_i);
1061	<	if (j_is_Dipole)
1062	<	*(idat.t2) -= cross(uz_j, indirect_duduz_j);
1147	>	if (a_is_Dipole || a_is_Quadrupole)
1148	>	*(idat.t1) += *(idat.sw) * indirect_Ta;
1149	>
1150	>	if (b_is_Dipole || b_is_Quadrupole)
1151	>	*(idat.t2) += *(idat.sw) * indirect_Tb;
1152		}
1064	–
1153		return;
1154	<	}
1154	>	}
1155
1156		void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1157	<	RealType mu1, preVal, self;
1157	>
1158		if (!initialized_) initialize();
1159
1160	<	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1161	<
1160	>	ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]];
1161	>
1162		// logicals
1163		bool i_is_Charge = data.is_Charge;
1164		bool i_is_Dipole = data.is_Dipole;
1165	+	bool i_is_Quadrupole = data.is_Quadrupole;
1166		bool i_is_Fluctuating = data.is_Fluctuating;
1167	<	RealType chg1 = data.fixedCharge;
1168	<
1080	<	if (i_is_Fluctuating) {
1081	<	chg1 += *(sdat.flucQ);
1082	<	// dVdFQ is really a force, so this is negative the derivative
1083	<	*(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1084	<	(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1085	<	(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1086	<	}
1167	>	RealType C_a = data.fixedCharge;
1168	>	RealType self(0.0), preVal, DdD(0.0), trQ, trQQ;
1169
1170	<	if (summationMethod_ == esm_REACTION_FIELD) {
1171	<	if (i_is_Dipole) {
1172	<	mu1 = data.dipole_moment;
1091	<	preVal = pre22_ * preRF2_ * mu1 * mu1;
1092	<	(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1170	>	if (i_is_Dipole) {
1171	>	DdD = data.dipole.lengthSquare();
1172	>	}
1173
1174	<	// The self-correction term adds into the reaction field vector
1175	<	Vector3d uz_i = sdat.eFrame->getColumn(2);
1096	<	Vector3d ei = preVal * uz_i;
1174	>	if (i_is_Fluctuating) {
1175	>	C_a += *(sdat.flucQ);
1176
1177	<	// This looks very wrong.  A vector crossed with itself is zero.
1178	<	*(sdat.t) -= cross(uz_i, ei);
1179	<	}
1180	<	} else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1181	<	if (i_is_Charge) {
1182	<	if (screeningMethod_ == DAMPED) {
1183	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1184	<	} else {
1185	<	self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1186	<	}
1187	<	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1177	>	flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ),
1178	>	(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY],
1179	>	*(sdat.flucQfrc) );
1180	>
1181	>	}
1182	>
1183	>	switch (summationMethod_) {
1184	>	case esm_REACTION_FIELD:
1185	>
1186	>	if (i_is_Charge) {
1187	>	// Self potential [see Wang and Hermans, "Reaction Field
1188	>	// Molecular Dynamics Simulation with Friedman’s Image Charge
1189	>	// Method," J. Phys. Chem. 99, 12001-12007 (1995).]
1190	>	preVal = pre11_ * preRF_ * C_a * C_a;
1191	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal / cutoffRadius_;
1192		}
1193	+
1194	+	if (i_is_Dipole) {
1195	+	(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD;
1196	+	}
1197	+
1198	+	break;
1199	+
1200	+	case esm_SHIFTED_FORCE:
1201	+	case esm_SHIFTED_POTENTIAL:
1202	+	case esm_TAYLOR_SHIFTED:
1203	+	case esm_EWALD_FULL:
1204	+	if (i_is_Charge)
1205	+	self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge));
1206	+	if (i_is_Dipole)
1207	+	self += selfMult2_ * pre22_ * DdD;
1208	+	if (i_is_Quadrupole) {
1209	+	trQ = data.quadrupole.trace();
1210	+	trQQ = (data.quadrupole * data.quadrupole).trace();
1211	+	self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ);
1212	+	if (i_is_Charge)
1213	+	self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ;
1214	+	}
1215	+	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1216	+	break;
1217	+	default:
1218	+	break;
1219		}
1220		}
1221	<
1221	>
1222		RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1223		// This seems to work moderately well as a default.  There's no
1224		// inherent scale for 1/r interactions that we can standardize.
#	Line 1117 | Line 1226 | namespace OpenMD {
1226		// cases.
1227		return 12.0;
1228		}
1229	+
1230	+
1231	+	void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1232	+
1233	+	RealType kPot = 0.0;
1234	+	RealType kVir = 0.0;
1235	+
1236	+	const RealType mPoleConverter = 0.20819434; // converts from the
1237	+	// internal units of
1238	+	// Debye (for dipoles)
1239	+	// or Debye-angstroms
1240	+	// (for quadrupoles) to
1241	+	// electron angstroms or
1242	+	// electron-angstroms^2
1243	+
1244	+	const RealType eConverter = 332.0637778; // convert the
1245	+	// Charge-Charge
1246	+	// electrostatic
1247	+	// interactions into kcal /
1248	+	// mol assuming distances
1249	+	// are measured in
1250	+	// angstroms.
1251	+
1252	+	Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
1253	+	Vector3d box = hmat.diagonals();
1254	+	RealType boxMax = box.max();
1255	+
1256	+	//int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) );
1257	+	int kMax = 7;
1258	+	int kSqMax = kMax*kMax + 2;
1259	+
1260	+	int kLimit = kMax+1;
1261	+	int kLim2 = 2*kMax+1;
1262	+	int kSqLim = kSqMax;
1263	+
1264	+	vector<RealType> AK(kSqLim+1, 0.0);
1265	+	RealType xcl = 2.0 * M_PI / box.x();
1266	+	RealType ycl = 2.0 * M_PI / box.y();
1267	+	RealType zcl = 2.0 * M_PI / box.z();
1268	+	RealType rcl = 2.0 * M_PI / boxMax;
1269	+	RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z());
1270	+
1271	+	if(dampingAlpha_ < 1.0e-12) return;
1272	+
1273	+	RealType ralph = -0.25/pow(dampingAlpha_,2);
1274	+
1275	+	// Calculate and store exponential factors
1276	+
1277	+	vector<vector<RealType> > elc;
1278	+	vector<vector<RealType> > emc;
1279	+	vector<vector<RealType> > enc;
1280	+	vector<vector<RealType> > els;
1281	+	vector<vector<RealType> > ems;
1282	+	vector<vector<RealType> > ens;
1283	+
1284	+	int nMax = info_->getNAtoms();
1285	+
1286	+	elc.resize(kLimit+1);
1287	+	emc.resize(kLimit+1);
1288	+	enc.resize(kLimit+1);
1289	+	els.resize(kLimit+1);
1290	+	ems.resize(kLimit+1);
1291	+	ens.resize(kLimit+1);
1292	+
1293	+	for (int j = 0; j < kLimit+1; j++) {
1294	+	elc[j].resize(nMax);
1295	+	emc[j].resize(nMax);
1296	+	enc[j].resize(nMax);
1297	+	els[j].resize(nMax);
1298	+	ems[j].resize(nMax);
1299	+	ens[j].resize(nMax);
1300	+	}
1301	+
1302	+	Vector3d t( 2.0 * M_PI );
1303	+	t.Vdiv(t, box);
1304	+
1305	+	SimInfo::MoleculeIterator mi;
1306	+	Molecule::AtomIterator ai;
1307	+	int i;
1308	+	Vector3d r;
1309	+	Vector3d tt;
1310	+
1311	+	for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1312	+	mol = info_->nextMolecule(mi)) {
1313	+	for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1314	+	atom = mol->nextAtom(ai)) {
1315	+
1316	+	i = atom->getLocalIndex();
1317	+	r = atom->getPos();
1318	+	info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
1319	+
1320	+	tt.Vmul(t, r);
1321	+
1322	+	elc[1][i] = 1.0;
1323	+	emc[1][i] = 1.0;
1324	+	enc[1][i] = 1.0;
1325	+	els[1][i] = 0.0;
1326	+	ems[1][i] = 0.0;
1327	+	ens[1][i] = 0.0;
1328	+
1329	+	elc[2][i] = cos(tt.x());
1330	+	emc[2][i] = cos(tt.y());
1331	+	enc[2][i] = cos(tt.z());
1332	+	els[2][i] = sin(tt.x());
1333	+	ems[2][i] = sin(tt.y());
1334	+	ens[2][i] = sin(tt.z());
1335	+
1336	+	for(int l = 3; l <= kLimit; l++) {
1337	+	elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i];
1338	+	emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i];
1339	+	enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i];
1340	+	els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i];
1341	+	ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i];
1342	+	ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i];
1343	+	}
1344	+	}
1345	+	}
1346	+
1347	+	// Calculate and store AK coefficients:
1348	+
1349	+	RealType eksq = 1.0;
1350	+	RealType expf = 0.0;
1351	+	if (ralph < 0.0) expf = exp(ralph*rcl*rcl);
1352	+	for (i = 1; i <= kSqLim; i++) {
1353	+	RealType rksq = float(i)*rcl*rcl;
1354	+	eksq = expf*eksq;
1355	+	AK[i] = eConverter * eksq/rksq;
1356	+	}
1357	+
1358	+	/*
1359	+	* Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz)
1360	+	* the values of ll, mm and nn are selected so that the symmetry of
1361	+	* reciprocal lattice is taken into account i.e. the following
1362	+	* rules apply.
1363	+	*
1364	+	* ll ranges over the values 0 to kMax only.
1365	+	*
1366	+	* mm ranges over 0 to kMax when ll=0 and over
1367	+	*            -kMax to kMax otherwise.
1368	+	* nn ranges over 1 to kMax when ll=mm=0 and over
1369	+	*            -kMax to kMax otherwise.
1370	+	*
1371	+	* Hence the result of the summation must be doubled at the end.
1372	+	*/
1373	+
1374	+	std::vector<RealType> clm(nMax, 0.0);
1375	+	std::vector<RealType> slm(nMax, 0.0);
1376	+	std::vector<RealType> ckr(nMax, 0.0);
1377	+	std::vector<RealType> skr(nMax, 0.0);
1378	+	std::vector<RealType> ckc(nMax, 0.0);
1379	+	std::vector<RealType> cks(nMax, 0.0);
1380	+	std::vector<RealType> dkc(nMax, 0.0);
1381	+	std::vector<RealType> dks(nMax, 0.0);
1382	+	std::vector<RealType> qkc(nMax, 0.0);
1383	+	std::vector<RealType> qks(nMax, 0.0);
1384	+	std::vector<Vector3d> dxk(nMax, V3Zero);
1385	+	std::vector<Vector3d> qxk(nMax, V3Zero);
1386	+	RealType rl, rm, rn;
1387	+	Vector3d kVec;
1388	+	Vector3d Qk;
1389	+	Mat3x3d k2;
1390	+	RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1391	+	int atid;
1392	+	ElectrostaticAtomData data;
1393	+	RealType C, dk, qk;
1394	+	Vector3d D;
1395	+	Mat3x3d  Q;
1396	+
1397	+	int mMin = kLimit;
1398	+	int nMin = kLimit + 1;
1399	+	for (int l = 1; l <= kLimit; l++) {
1400	+	int ll = l - 1;
1401	+	rl = xcl * float(ll);
1402	+	for (int mmm = mMin; mmm <= kLim2; mmm++) {
1403	+	int mm = mmm - kLimit;
1404	+	int m = abs(mm) + 1;
1405	+	rm = ycl * float(mm);
1406	+	// Set temporary products of exponential terms
1407	+	for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1408	+	mol = info_->nextMolecule(mi)) {
1409	+	for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1410	+	atom = mol->nextAtom(ai)) {
1411	+
1412	+	i = atom->getLocalIndex();
1413	+	if(mm < 0) {
1414	+	clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i];
1415	+	slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i];
1416	+	} else {
1417	+	clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i];
1418	+	slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i];
1419	+	}
1420	+	}
1421	+	}
1422	+	for (int nnn = nMin; nnn <= kLim2; nnn++) {
1423	+	int nn = nnn - kLimit;
1424	+	int n = abs(nn) + 1;
1425	+	rn = zcl * float(nn);
1426	+	// Test on magnitude of k vector:
1427	+	int kk=ll*ll + mm*mm + nn*nn;
1428	+	if(kk <= kSqLim) {
1429	+	kVec = Vector3d(rl, rm, rn);
1430	+	k2 = outProduct(kVec, kVec);
1431	+	// Calculate exp(ikr) terms
1432	+	for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1433	+	mol = info_->nextMolecule(mi)) {
1434	+	for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1435	+	atom = mol->nextAtom(ai)) {
1436	+	i = atom->getLocalIndex();
1437	+
1438	+	if (nn < 0) {
1439	+	ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i];
1440	+	skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i];
1441	+
1442	+	} else {
1443	+	ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i];
1444	+	skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i];
1445	+	}
1446	+	}
1447	+	}
1448	+
1449	+	// Calculate scalar and vector products for each site:
1450	+
1451	+	for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1452	+	mol = info_->nextMolecule(mi)) {
1453	+	for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1454	+	atom = mol->nextAtom(ai)) {
1455	+	i = atom->getLocalIndex();
1456	+	int atid = atom->getAtomType()->getIdent();
1457	+	data = ElectrostaticMap[Etids[atid]];
1458	+
1459	+	if (data.is_Charge) {
1460	+	C = data.fixedCharge;
1461	+	if (atom->isFluctuatingCharge()) C += atom->getFlucQPos();
1462	+	ckc[i] = C * ckr[i];
1463	+	cks[i] = C * skr[i];
1464	+	}
1465	+
1466	+	if (data.is_Dipole) {
1467	+	D = atom->getDipole() * mPoleConverter;
1468	+	dk = dot(D, kVec);
1469	+	dxk[i] = cross(D, kVec);
1470	+	dkc[i] = dk * ckr[i];
1471	+	dks[i] = dk * skr[i];
1472	+	}
1473	+	if (data.is_Quadrupole) {
1474	+	Q = atom->getQuadrupole() * mPoleConverter;
1475	+	Qk = Q * kVec;
1476	+	qk = dot(kVec, Qk);
1477	+	qxk[i] = -cross(kVec, Qk);
1478	+	qkc[i] = qk * ckr[i];
1479	+	qks[i] = qk * skr[i];
1480	+	}
1481	+	}
1482	+	}
1483	+
1484	+	// calculate vector sums
1485	+
1486	+	ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0);
1487	+	ckss = std::accumulate(cks.begin(),cks.end(),0.0);
1488	+	dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0);
1489	+	dkss = std::accumulate(dks.begin(),dks.end(),0.0);
1490	+	qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0);
1491	+	qkss = std::accumulate(qks.begin(),qks.end(),0.0);
1492	+
1493	+	#ifdef IS_MPI
1494	+	MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE,
1495	+	MPI_SUM, MPI_COMM_WORLD);
1496	+	MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE,
1497	+	MPI_SUM, MPI_COMM_WORLD);
1498	+	MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE,
1499	+	MPI_SUM, MPI_COMM_WORLD);
1500	+	MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE,
1501	+	MPI_SUM, MPI_COMM_WORLD);
1502	+	MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE,
1503	+	MPI_SUM, MPI_COMM_WORLD);
1504	+	MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE,
1505	+	MPI_SUM, MPI_COMM_WORLD);
1506	+	#endif
1507	+
1508	+	// Accumulate potential energy and virial contribution:
1509	+
1510	+	kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss)
1511	+	+ (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs));
1512	+
1513	+	kVir += 2.0 * rvol  * AK[kk]*(ckcs*ckcs+ckss*ckss
1514	+	+4.0*(ckss*dkcs-ckcs*dkss)
1515	+	+3.0*(dkcs*dkcs+dkss*dkss)
1516	+	-6.0*(ckss*qkss+ckcs*qkcs)
1517	+	+8.0*(dkss*qkcs-dkcs*qkss)
1518	+	+5.0*(qkss*qkss+qkcs*qkcs));
1519	+
1520	+	// Calculate force and torque for each site:
1521	+
1522	+	for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1523	+	mol = info_->nextMolecule(mi)) {
1524	+	for(Atom* atom = mol->beginAtom(ai); atom != NULL;
1525	+	atom = mol->nextAtom(ai)) {
1526	+
1527	+	i = atom->getLocalIndex();
1528	+	atid = atom->getAtomType()->getIdent();
1529	+	data = ElectrostaticMap[Etids[atid]];
1530	+
1531	+	RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs)
1532	+	- (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss));
1533	+	RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs)
1534	+	-ckr[i]*(ckss+dkcs-qkss));
1535	+	RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)
1536	+	+skr[i]*(ckss+dkcs-qkss));
1537	+
1538	+	atom->addFrc( 4.0 * rvol * qfrc * kVec );
1539	+
1540	+	if (atom->isFluctuatingCharge()) {
1541	+	atom->addFlucQFrc( - 2.0 * rvol * qtrq2 );
1542	+	}
1543	+
1544	+	if (data.is_Dipole) {
1545	+	atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] );
1546	+	}
1547	+	if (data.is_Quadrupole) {
1548	+	atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] );
1549	+	}
1550	+	}
1551	+	}
1552	+	}
1553	+	}
1554	+	nMin = 1;
1555	+	}
1556	+	mMin = 1;
1557	+	}
1558	+	pot += kPot;
1559	+	}
1560	+
1561	+	void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded,
1562	+	RealType &spot1, RealType &spot2) {
1563	+
1564	+	if (!initialized_) {
1565	+	cerr << "initializing\n";
1566	+	initialize();
1567	+	cerr << "done\n";
1568	+	}
1569	+
1570	+	const RealType mPoleConverter = 0.20819434;
1571	+
1572	+	AtomType* atype1 = a1->getAtomType();
1573	+	AtomType* atype2 = a2->getAtomType();
1574	+	int atid1 = atype1->getIdent();
1575	+	int atid2 = atype2->getIdent();
1576	+	data1 = ElectrostaticMap[Etids[atid1]];
1577	+	data2 = ElectrostaticMap[Etids[atid2]];
1578	+
1579	+	Pa = 0.0;  // Site potential at site a
1580	+	Pb = 0.0;  // Site potential at site b
1581	+
1582	+	Vector3d d = a2->getPos() - a1->getPos();
1583	+	info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d);
1584	+	RealType rij = d.length();
1585	+	// some variables we'll need independent of electrostatic type:
1586	+
1587	+	RealType ri = 1.0 /  rij;
1588	+	rhat = d  * ri;
1589	+
1590	+
1591	+	if ((rij >= cutoffRadius_) || excluded) {
1592	+	spot1 = 0.0;
1593	+	spot2 = 0.0;
1594	+	return;
1595	+	}
1596	+
1597	+	// logicals
1598	+
1599	+	a_is_Charge = data1.is_Charge;
1600	+	a_is_Dipole = data1.is_Dipole;
1601	+	a_is_Quadrupole = data1.is_Quadrupole;
1602	+	a_is_Fluctuating = data1.is_Fluctuating;
1603	+
1604	+	b_is_Charge = data2.is_Charge;
1605	+	b_is_Dipole = data2.is_Dipole;
1606	+	b_is_Quadrupole = data2.is_Quadrupole;
1607	+	b_is_Fluctuating = data2.is_Fluctuating;
1608	+
1609	+	// Obtain all of the required radial function values from the
1610	+	// spline structures:
1611	+
1612	+
1613	+	if (a_is_Charge || b_is_Charge) {
1614	+	v01 = v01s->getValueAt(rij);
1615	+	}
1616	+	if (a_is_Dipole || b_is_Dipole) {
1617	+	v11 = v11s->getValueAt(rij);
1618	+	v11or = ri * v11;
1619	+	}
1620	+	if (a_is_Quadrupole || b_is_Quadrupole) {
1621	+	v21 = v21s->getValueAt(rij);
1622	+	v22 = v22s->getValueAt(rij);
1623	+	v22or = ri * v22;
1624	+	}
1625	+
1626	+	if (a_is_Charge) {
1627	+	C_a = data1.fixedCharge;
1628	+
1629	+	if (a_is_Fluctuating) {
1630	+	C_a += a1->getFlucQPos();
1631	+	}
1632	+
1633	+	Pb += C_a *  pre11_ * v01;
1634	+	}
1635	+
1636	+	if (a_is_Dipole) {
1637	+	D_a = a1->getDipole() * mPoleConverter;
1638	+	rdDa = dot(rhat, D_a);
1639	+	Pb +=  pre12_ * v11 * rdDa;
1640	+	}
1641	+
1642	+	if (a_is_Quadrupole) {
1643	+	Q_a = a1->getQuadrupole() * mPoleConverter;
1644	+	trQa =  Q_a.trace();
1645	+	Qar =   Q_a * rhat;
1646	+	rdQar = dot(rhat, Qar);
1647	+	Pb += pre14_ * (v21 * trQa + v22 * rdQar);
1648	+	}
1649	+
1650	+	if (b_is_Charge) {
1651	+	C_b = data2.fixedCharge;
1652	+
1653	+	if (b_is_Fluctuating)
1654	+	C_b += a2->getFlucQPos();
1655	+
1656	+	Pa += C_b *  pre11_ * v01;
1657	+	}
1658	+
1659	+	if (b_is_Dipole) {
1660	+	D_b = a2->getDipole() * mPoleConverter;
1661	+	rdDb = dot(rhat, D_b);
1662	+	Pa += pre12_ * v11 * rdDb;
1663	+	}
1664	+
1665	+	if (b_is_Quadrupole) {
1666	+	Q_a = a2->getQuadrupole() * mPoleConverter;
1667	+	trQb =  Q_b.trace();
1668	+	Qbr =   Q_b * rhat;
1669	+	rdQbr = dot(rhat, Qbr);
1670	+	Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
1671	+	}
1672	+
1673	+	spot1 = Pa;
1674	+	spot2 = Pb;
1675	+	}
1676		}

Diff Legend

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