ViewVC Help

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

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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1750 by gezelter, Thu Jun 7 12:53:46 2012 UTC vs.
Revision 1808 by gezelter, Mon Oct 22 20:42:10 2012 UTC

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

Diff Legend

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