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 1586 by gezelter, Tue Jun 21 06:34:35 2011 UTC vs.
trunk/src/nonbonded/Electrostatic.cpp (file contents), Revision 1981 by gezelter, Mon Apr 14 18:32:51 2014 UTC

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

Diff Legend

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