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 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs.
trunk/src/nonbonded/Electrostatic.cpp (file contents), Revision 1925 by gezelter, Wed Aug 7 15:24:16 2013 UTC

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

Diff Legend

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