ViewVC Help

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

root/OpenMD/trunk/src/nonbonded/Electrostatic.cpp

Comparing:
branches/development/src/nonbonded/Electrostatic.cpp (file contents), Revision 1668 by gezelter, Fri Jan 6 19:03:05 2012 UTC vs.
trunk/src/nonbonded/Electrostatic.cpp (file contents), Revision 1915 by gezelter, Mon Jul 29 17:55:17 2013 UTC

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

Diff Legend

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