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 1616 by gezelter, Tue Aug 30 15:45:35 2011 UTC vs.
Revision 1825 by gezelter, Wed Jan 9 19:27:52 2013 UTC

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

Diff Legend

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