[ViewVC] Diff of: OpenMD/branches/development/src/nonbonded/Electrostatic.cpp

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs.
Revision 1767 by gezelter, Fri Jul 6 22:01:58 2012 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).
37	>	* [2] Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).
38		* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	<	* [4] Vardeman & Gezelter, in progress (2009).
39	>	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	>	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		#include <stdio.h>
#	Line 46 \| Line 47
47		#include "nonbonded/Electrostatic.hpp"
48		#include "utils/simError.h"
49		#include "types/NonBondedInteractionType.hpp"
50	<	#include "types/DirectionalAtomType.hpp"
50	>	#include "types/FixedChargeAdapter.hpp"
51	>	#include "types/FluctuatingChargeAdapter.hpp"
52	>	#include "types/MultipoleAdapter.hpp"
53	>	#include "io/Globals.hpp"
54	>	#include "nonbonded/SlaterIntegrals.hpp"
55	>	#include "utils/PhysicalConstants.hpp"
56	>	#include "math/erfc.hpp"
57
51	–
58		namespace OpenMD {
59
60		Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
61	<	forceField_(NULL) {}
61	>	forceField_(NULL), info_(NULL),
62	>	haveCutoffRadius_(false),
63	>	haveDampingAlpha_(false),
64	>	haveDielectric_(false),
65	>	haveElectroSpline_(false)
66	>	{}
67
68		void Electrostatic::initialize() {
69	+
70	+	Globals* simParams_ = info_->getSimParams();
71	+
72	+	summationMap_["HARD"] = esm_HARD;
73	+	summationMap_["NONE"] = esm_HARD;
74	+	summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
75	+	summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
76	+	summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
77	+	summationMap_["REACTION_FIELD"] = esm_REACTION_FIELD;
78	+	summationMap_["EWALD_FULL"] = esm_EWALD_FULL;
79	+	summationMap_["EWALD_PME"] = esm_EWALD_PME;
80	+	summationMap_["EWALD_SPME"] = esm_EWALD_SPME;
81	+	screeningMap_["DAMPED"] = DAMPED;
82	+	screeningMap_["UNDAMPED"] = UNDAMPED;
83	+
84		// these prefactors convert the multipole interactions into kcal / mol
85		// all were computed assuming distances are measured in angstroms
86		// Charge-Charge, assuming charges are measured in electrons
#	Line 79 \| Line 105 \| namespace OpenMD {
105
106		// variables to handle different summation methods for long-range
107		// electrostatics:
108	<	summationMethod_ = NONE;
108	>	summationMethod_ = esm_HARD;
109		screeningMethod_ = UNDAMPED;
110		dielectric_ = 1.0;
111		one_third_ = 1.0 / 3.0;
86	–	haveDefaultCutoff_ = false;
87	–	haveDampingAlpha_ = false;
88	–	haveDielectric_ = false;
89	–	haveElectroSpline_ = false;
112
113	+	// check the summation method:
114	+	if (simParams_->haveElectrostaticSummationMethod()) {
115	+	string myMethod = simParams_->getElectrostaticSummationMethod();
116	+	toUpper(myMethod);
117	+	map<string, ElectrostaticSummationMethod>::iterator i;
118	+	i = summationMap_.find(myMethod);
119	+	if ( i != summationMap_.end() ) {
120	+	summationMethod_ = (*i).second;
121	+	} else {
122	+	// throw error
123	+	sprintf( painCave.errMsg,
124	+	"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
125	+	"\t(Input file specified %s .)\n"
126	+	"\telectrostaticSummationMethod must be one of: \"hard\",\n"
127	+	"\t\"shifted_potential\", \"shifted_force\", or \n"
128	+	"\t\"reaction_field\".\n", myMethod.c_str() );
129	+	painCave.isFatal = 1;
130	+	simError();
131	+	}
132	+	} else {
133	+	// set ElectrostaticSummationMethod to the cutoffMethod:
134	+	if (simParams_->haveCutoffMethod()){
135	+	string myMethod = simParams_->getCutoffMethod();
136	+	toUpper(myMethod);
137	+	map<string, ElectrostaticSummationMethod>::iterator i;
138	+	i = summationMap_.find(myMethod);
139	+	if ( i != summationMap_.end() ) {
140	+	summationMethod_ = (*i).second;
141	+	}
142	+	}
143	+	}
144	+
145	+	if (summationMethod_ == esm_REACTION_FIELD) {
146	+	if (!simParams_->haveDielectric()) {
147	+	// throw warning
148	+	sprintf( painCave.errMsg,
149	+	"SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
150	+	"\tA default value of %f will be used for the dielectric.\n", dielectric_);
151	+	painCave.isFatal = 0;
152	+	painCave.severity = OPENMD_INFO;
153	+	simError();
154	+	} else {
155	+	dielectric_ = simParams_->getDielectric();
156	+	}
157	+	haveDielectric_ = true;
158	+	}
159	+
160	+	if (simParams_->haveElectrostaticScreeningMethod()) {
161	+	string myScreen = simParams_->getElectrostaticScreeningMethod();
162	+	toUpper(myScreen);
163	+	map<string, ElectrostaticScreeningMethod>::iterator i;
164	+	i = screeningMap_.find(myScreen);
165	+	if ( i != screeningMap_.end()) {
166	+	screeningMethod_ = (*i).second;
167	+	} else {
168	+	sprintf( painCave.errMsg,
169	+	"SimInfo error: Unknown electrostaticScreeningMethod.\n"
170	+	"\t(Input file specified %s .)\n"
171	+	"\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
172	+	"or \"damped\".\n", myScreen.c_str() );
173	+	painCave.isFatal = 1;
174	+	simError();
175	+	}
176	+	}
177	+
178	+	// check to make sure a cutoff value has been set:
179	+	if (!haveCutoffRadius_) {
180	+	sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
181	+	"Cutoff value!\n");
182	+	painCave.severity = OPENMD_ERROR;
183	+	painCave.isFatal = 1;
184	+	simError();
185	+	}
186	+
187	+	if (screeningMethod_ == DAMPED) {
188	+	if (!simParams_->haveDampingAlpha()) {
189	+	// first set a cutoff dependent alpha value
190	+	// we assume alpha depends linearly with rcut from 0 to 20.5 ang
191	+	dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
192	+	if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
193	+
194	+	// throw warning
195	+	sprintf( painCave.errMsg,
196	+	"Electrostatic::initialize: dampingAlpha was not specified in the\n"
197	+	"\tinput file. A default value of %f (1/ang) will be used for the\n"
198	+	"\tcutoff of %f (ang).\n",
199	+	dampingAlpha_, cutoffRadius_);
200	+	painCave.severity = OPENMD_INFO;
201	+	painCave.isFatal = 0;
202	+	simError();
203	+	} else {
204	+	dampingAlpha_ = simParams_->getDampingAlpha();
205	+	}
206	+	haveDampingAlpha_ = true;
207	+	}
208	+
209		// find all of the Electrostatic atom Types:
210		ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
211		ForceField::AtomTypeContainer::MapTypeIterator i;
212		AtomType* at;
213	<
213	>
214		for (at = atomTypes->beginType(i); at != NULL;
215		at = atomTypes->nextType(i)) {
216
#	Line 100 \| Line 218 \| namespace OpenMD {
218		addType(at);
219		}
220
221	<	// check to make sure a cutoff value has been set:
222	<	if (!haveDefaultCutoff_) {
105	<	sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
106	<	"Cutoff value!\n");
107	<	painCave.severity = OPENMD_ERROR;
108	<	painCave.isFatal = 1;
109	<	simError();
110	<	}
111	<
112	<	defaultCutoff2_ = defaultCutoff_ * defaultCutoff_;
113	<	rcuti_ = 1.0 / defaultCutoff_;
221	>	cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
222	>	rcuti_ = 1.0 / cutoffRadius_;
223		rcuti2_ = rcuti_ * rcuti_;
224		rcuti3_ = rcuti2_ * rcuti_;
225		rcuti4_ = rcuti2_ * rcuti2_;
226
227		if (screeningMethod_ == DAMPED) {
228	<	if (!haveDampingAlpha_) {
120	<	sprintf( painCave.errMsg, "Electrostatic::initialize has no "
121	<	"DampingAlpha value!\n");
122	<	painCave.severity = OPENMD_ERROR;
123	<	painCave.isFatal = 1;
124	<	simError();
125	<	}
126	<
228	>
229		alpha2_ = dampingAlpha_ * dampingAlpha_;
230		alpha4_ = alpha2_ * alpha2_;
231		alpha6_ = alpha4_ * alpha2_;
232		alpha8_ = alpha4_ * alpha4_;
233
234	<	constEXP_ = exp(-alpha2_ * defaultCutoff2_);
234	>	constEXP_ = exp(-alpha2_ * cutoffRadius2_);
235		invRootPi_ = 0.56418958354775628695;
236		alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
237
238	<	c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * rcuti_;
238	>	c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
239		c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
240		c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
241		c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_;
#	Line 148 \| Line 250 \| namespace OpenMD {
250		c6c_ = 9.0 * c5c_ * rcuti2_;
251		}
252
253	<	if (summationMethod_ == REACTION_FIELD) {
254	<	if (haveDielectric_) {
255	<	preRF_ = (dielectric_ - 1.0) /
256	<	((2.0 * dielectric_ + 1.0) * defaultCutoff2_ * defaultCutoff_);
155	<	preRF2_ = 2.0 * preRF_;
156	<	} else {
157	<	sprintf( painCave.errMsg, "Electrostatic::initialize has no Dielectric"
158	<	" value!\n");
159	<	painCave.severity = OPENMD_ERROR;
160	<	painCave.isFatal = 1;
161	<	simError();
162	<	}
253	>	if (summationMethod_ == esm_REACTION_FIELD) {
254	>	preRF_ = (dielectric_ - 1.0) /
255	>	((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
256	>	preRF2_ = 2.0 * preRF_;
257		}
258	<
259	<	RealType dx = defaultCutoff_ / RealType(np_ - 1);
258	>
259	>	// Add a 2 angstrom safety window to deal with cutoffGroups that
260	>	// have charged atoms longer than the cutoffRadius away from each
261	>	// other. Splining may not be the best choice here. Direct calls
262	>	// to erfc might be preferrable.
263	>
264	>	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
265		RealType rval;
266		vector<RealType> rvals;
267		vector<RealType> yvals;
#	Line 185 \| Line 284 \| namespace OpenMD {
284		electrostaticAtomData.is_Dipole = false;
285		electrostaticAtomData.is_SplitDipole = false;
286		electrostaticAtomData.is_Quadrupole = false;
287	+	electrostaticAtomData.is_Fluctuating = false;
288
289	<	if (atomType->isCharge()) {
190	<	GenericData* data = atomType->getPropertyByName("Charge");
289	>	FixedChargeAdapter fca = FixedChargeAdapter(atomType);
290
291	<	if (data == NULL) {
193	<	sprintf( painCave.errMsg, "Electrostatic::addType could not find "
194	<	"Charge\n"
195	<	"\tparameters for atomType %s.\n",
196	<	atomType->getName().c_str());
197	<	painCave.severity = OPENMD_ERROR;
198	<	painCave.isFatal = 1;
199	<	simError();
200	<	}
201	<
202	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
203	<	if (doubleData == NULL) {
204	<	sprintf( painCave.errMsg,
205	<	"Electrostatic::addType could not convert GenericData to "
206	<	"Charge for\n"
207	<	"\tatom type %s\n", atomType->getName().c_str());
208	<	painCave.severity = OPENMD_ERROR;
209	<	painCave.isFatal = 1;
210	<	simError();
211	<	}
291	>	if (fca.isFixedCharge()) {
292		electrostaticAtomData.is_Charge = true;
293	<	electrostaticAtomData.charge = doubleData->getData();
293	>	electrostaticAtomData.fixedCharge = fca.getCharge();
294		}
295
296	<	if (atomType->isDirectional()) {
297	<	DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
298	<
219	<	if (daType->isDipole()) {
220	<	GenericData* data = daType->getPropertyByName("Dipole");
221	<
222	<	if (data == NULL) {
223	<	sprintf( painCave.errMsg,
224	<	"Electrostatic::addType could not find Dipole\n"
225	<	"\tparameters for atomType %s.\n",
226	<	daType->getName().c_str());
227	<	painCave.severity = OPENMD_ERROR;
228	<	painCave.isFatal = 1;
229	<	simError();
230	<	}
231	<
232	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
233	<	if (doubleData == NULL) {
234	<	sprintf( painCave.errMsg,
235	<	"Electrostatic::addType could not convert GenericData to "
236	<	"Dipole Moment\n"
237	<	"\tfor atom type %s\n", daType->getName().c_str());
238	<	painCave.severity = OPENMD_ERROR;
239	<	painCave.isFatal = 1;
240	<	simError();
241	<	}
296	>	MultipoleAdapter ma = MultipoleAdapter(atomType);
297	>	if (ma.isMultipole()) {
298	>	if (ma.isDipole()) {
299		electrostaticAtomData.is_Dipole = true;
300	<	electrostaticAtomData.dipole_moment = doubleData->getData();
301	<	}
302	<
246	<	if (daType->isSplitDipole()) {
247	<	GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
248	<
249	<	if (data == NULL) {
250	<	sprintf(painCave.errMsg,
251	<	"Electrostatic::addType could not find SplitDipoleDistance\n"
252	<	"\tparameter for atomType %s.\n",
253	<	daType->getName().c_str());
254	<	painCave.severity = OPENMD_ERROR;
255	<	painCave.isFatal = 1;
256	<	simError();
257	<	}
258	<
259	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
260	<	if (doubleData == NULL) {
261	<	sprintf( painCave.errMsg,
262	<	"Electrostatic::addType could not convert GenericData to "
263	<	"SplitDipoleDistance for\n"
264	<	"\tatom type %s\n", daType->getName().c_str());
265	<	painCave.severity = OPENMD_ERROR;
266	<	painCave.isFatal = 1;
267	<	simError();
268	<	}
300	>	electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
301	>	}
302	>	if (ma.isSplitDipole()) {
303		electrostaticAtomData.is_SplitDipole = true;
304	<	electrostaticAtomData.split_dipole_distance = doubleData->getData();
304	>	electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
305		}
306	<
307	<	if (daType->isQuadrupole()) {
308	<	GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
309	<
310	<	if (data == NULL) {
311	<	sprintf( painCave.errMsg,
278	<	"Electrostatic::addType could not find QuadrupoleMoments\n"
279	<	"\tparameter for atomType %s.\n",
280	<	daType->getName().c_str());
281	<	painCave.severity = OPENMD_ERROR;
282	<	painCave.isFatal = 1;
283	<	simError();
284	<	}
285	<
286	<	Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
287	<	if (v3dData == NULL) {
288	<	sprintf( painCave.errMsg,
289	<	"Electrostatic::addType could not convert GenericData to "
290	<	"Quadrupole Moments for\n"
291	<	"\tatom type %s\n", daType->getName().c_str());
292	<	painCave.severity = OPENMD_ERROR;
293	<	painCave.isFatal = 1;
294	<	simError();
295	<	}
306	>	if (ma.isQuadrupole()) {
307	>	// Quadrupoles in OpenMD are set as the diagonal elements
308	>	// of the diagonalized traceless quadrupole moment tensor.
309	>	// The column vectors of the unitary matrix that diagonalizes
310	>	// the quadrupole moment tensor become the eFrame (or the
311	>	// electrostatic version of the body-fixed frame.
312		electrostaticAtomData.is_Quadrupole = true;
313	<	electrostaticAtomData.quadrupole_moments = v3dData->getData();
313	>	electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
314		}
315		}
316
317	<	AtomTypeProperties atp = atomType->getATP();
317	>	FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
318
319	+	if (fqa.isFluctuatingCharge()) {
320	+	electrostaticAtomData.is_Fluctuating = true;
321	+	electrostaticAtomData.electronegativity = fqa.getElectronegativity();
322	+	electrostaticAtomData.hardness = fqa.getHardness();
323	+	electrostaticAtomData.slaterN = fqa.getSlaterN();
324	+	electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
325	+	}
326	+
327		pair<map<int,AtomType*>::iterator,bool> ret;
328	<	ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
328	>	ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
329	>	atomType) );
330		if (ret.second == false) {
331		sprintf( painCave.errMsg,
332		"Electrostatic already had a previous entry with ident %d\n",
333	<	atp.ident);
333	>	atomType->getIdent() );
334		painCave.severity = OPENMD_INFO;
335		painCave.isFatal = 0;
336		simError();
337		}
338
339	<	ElectrostaticMap[atomType] = electrostaticAtomData;
339	>	ElectrostaticMap[atomType] = electrostaticAtomData;
340	>
341	>	// Now, iterate over all known types and add to the mixing map:
342	>
343	>	map<AtomType*, ElectrostaticAtomData>::iterator it;
344	>	for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
345	>	AtomType* atype2 = (*it).first;
346	>	ElectrostaticAtomData eaData2 = (*it).second;
347	>	if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
348	>
349	>	RealType a = electrostaticAtomData.slaterZeta;
350	>	RealType b = eaData2.slaterZeta;
351	>	int m = electrostaticAtomData.slaterN;
352	>	int n = eaData2.slaterN;
353	>
354	>	// Create the spline of the coulombic integral for s-type
355	>	// Slater orbitals. Add a 2 angstrom safety window to deal
356	>	// with cutoffGroups that have charged atoms longer than the
357	>	// cutoffRadius away from each other.
358	>
359	>	RealType rval;
360	>	RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
361	>	vector<RealType> rvals;
362	>	vector<RealType> J1vals;
363	>	vector<RealType> J2vals;
364	>	// don't start at i = 0, as rval = 0 is undefined for the slater overlap integrals.
365	>	for (int i = 1; i < np_+1; i++) {
366	>	rval = RealType(i) * dr;
367	>	rvals.push_back(rval);
368	>	J1vals.push_back(sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
369	>	// may not be necessary if Slater coulomb integral is symmetric
370	>	J2vals.push_back(sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
371	>	}
372	>
373	>	CubicSpline* J1 = new CubicSpline();
374	>	J1->addPoints(rvals, J1vals);
375	>	CubicSpline* J2 = new CubicSpline();
376	>	J2->addPoints(rvals, J2vals);
377	>
378	>	pair<AtomType, AtomType> key1, key2;
379	>	key1 = make_pair(atomType, atype2);
380	>	key2 = make_pair(atype2, atomType);
381	>
382	>	Jij[key1] = J1;
383	>	Jij[key2] = J2;
384	>	}
385	>	}
386	>
387		return;
388		}
389
390	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
391	<	RealType theRSW ) {
392	<	defaultCutoff_ = theECR;
393	<	rrf_ = defaultCutoff_;
322	<	rt_ = theRSW;
323	<	haveDefaultCutoff_ = true;
390	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
391	>	cutoffRadius_ = rCut;
392	>	rrf_ = cutoffRadius_;
393	>	haveCutoffRadius_ = true;
394		}
395	+
396	+	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
397	+	rt_ = rSwitch;
398	+	}
399		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
400		summationMethod_ = esm;
401		}
#	Line 337 \| Line 411 \| namespace OpenMD {
411		haveDielectric_ = true;
412		}
413
414	<	void Electrostatic::calcForce(InteractionData idat) {
414	>	void Electrostatic::calcForce(InteractionData &idat) {
415
416		// utility variables. Should clean these up and use the Vector3d and
417		// Mat3x3d to replace as many as we can in future versions:
#	Line 351 \| Line 425 \| namespace OpenMD {
425		RealType ct_i, ct_j, ct_ij, a1;
426		RealType riji, ri, ri2, ri3, ri4;
427		RealType pref, vterm, epot, dudr;
428	+	RealType vpair(0.0);
429		RealType scale, sc2;
430		RealType pot_term, preVal, rfVal;
431		RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
432		RealType preSw, preSwSc;
433		RealType c1, c2, c3, c4;
434	<	RealType erfcVal, derfcVal;
434	>	RealType erfcVal(1.0), derfcVal(0.0);
435		RealType BigR;
436	+	RealType two(2.0), three(3.0);
437
438		Vector3d Q_i, Q_j;
439		Vector3d ux_i, uy_i, uz_i;
#	Line 367 \| Line 443 \| namespace OpenMD {
443		Vector3d rhatdot2, rhatc4;
444		Vector3d dVdr;
445
446	+	// variables for indirect (reaction field) interactions for excluded pairs:
447	+	RealType indirect_Pot(0.0);
448	+	RealType indirect_vpair(0.0);
449	+	Vector3d indirect_dVdr(V3Zero);
450	+	Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
451	+
452	+	RealType coulInt, vFluc1(0.0), vFluc2(0.0);
453		pair<RealType, RealType> res;
454
455	+	// splines for coulomb integrals
456	+	CubicSpline* J1;
457	+	CubicSpline* J2;
458	+
459		if (!initialized_) initialize();
460
461	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
462	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
461	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
462	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
463
464		// some variables we'll need independent of electrostatic type:
465
466	<	riji = 1.0 / idat.rij;
467	<	Vector3d rhat = idat.d * riji;
466	>	riji = 1.0 / *(idat.rij) ;
467	>	Vector3d rhat = (idat.d) riji;
468
469		// logicals
470
#	Line 385 \| Line 472 \| namespace OpenMD {
472		bool i_is_Dipole = data1.is_Dipole;
473		bool i_is_SplitDipole = data1.is_SplitDipole;
474		bool i_is_Quadrupole = data1.is_Quadrupole;
475	+	bool i_is_Fluctuating = data1.is_Fluctuating;
476
477		bool j_is_Charge = data2.is_Charge;
478		bool j_is_Dipole = data2.is_Dipole;
479		bool j_is_SplitDipole = data2.is_SplitDipole;
480		bool j_is_Quadrupole = data2.is_Quadrupole;
481	+	bool j_is_Fluctuating = data2.is_Fluctuating;
482
483	<	if (i_is_Charge)
484	<	q_i = data1.charge;
483	>	if (i_is_Charge) {
484	>	q_i = data1.fixedCharge;
485
486	+	if (i_is_Fluctuating) {
487	+	q_i += *(idat.flucQ1);
488	+	}
489	+
490	+	if (idat.excluded) {
491	+	*(idat.skippedCharge2) += q_i;
492	+	}
493	+	}
494	+
495		if (i_is_Dipole) {
496		mu_i = data1.dipole_moment;
497	<	uz_i = idat.eFrame1.getColumn(2);
497	>	uz_i = idat.eFrame1->getColumn(2);
498
499		ct_i = dot(uz_i, rhat);
500
#	Line 412 \| Line 510 \| namespace OpenMD {
510		qyy_i = Q_i.y();
511		qzz_i = Q_i.z();
512
513	<	ux_i = idat.eFrame1.getColumn(0);
514	<	uy_i = idat.eFrame1.getColumn(1);
515	<	uz_i = idat.eFrame1.getColumn(2);
513	>	ux_i = idat.eFrame1->getColumn(0);
514	>	uy_i = idat.eFrame1->getColumn(1);
515	>	uz_i = idat.eFrame1->getColumn(2);
516
517		cx_i = dot(ux_i, rhat);
518		cy_i = dot(uy_i, rhat);
#	Line 425 \| Line 523 \| namespace OpenMD {
523		duduz_i = V3Zero;
524		}
525
526	<	if (j_is_Charge)
527	<	q_j = data2.charge;
526	>	if (j_is_Charge) {
527	>	q_j = data2.fixedCharge;
528
529	+	if (j_is_Fluctuating)
530	+	q_j += *(idat.flucQ2);
531	+
532	+	if (idat.excluded) {
533	+	*(idat.skippedCharge1) += q_j;
534	+	}
535	+	}
536	+
537	+
538		if (j_is_Dipole) {
539		mu_j = data2.dipole_moment;
540	<	uz_j = idat.eFrame2.getColumn(2);
540	>	uz_j = idat.eFrame2->getColumn(2);
541
542		ct_j = dot(uz_j, rhat);
543
#	Line 446 \| Line 553 \| namespace OpenMD {
553		qyy_j = Q_j.y();
554		qzz_j = Q_j.z();
555
556	<	ux_j = idat.eFrame2.getColumn(0);
557	<	uy_j = idat.eFrame2.getColumn(1);
558	<	uz_j = idat.eFrame2.getColumn(2);
556	>	ux_j = idat.eFrame2->getColumn(0);
557	>	uy_j = idat.eFrame2->getColumn(1);
558	>	uz_j = idat.eFrame2->getColumn(2);
559
560		cx_j = dot(ux_j, rhat);
561		cy_j = dot(uy_j, rhat);
#	Line 459 \| Line 566 \| namespace OpenMD {
566		duduz_j = V3Zero;
567		}
568
569	+	if (i_is_Fluctuating && j_is_Fluctuating) {
570	+	J1 = Jij[idat.atypes];
571	+	J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
572	+	}
573	+
574		epot = 0.0;
575		dVdr = V3Zero;
576
#	Line 467 \| Line 579 \| namespace OpenMD {
579		if (j_is_Charge) {
580		if (screeningMethod_ == DAMPED) {
581		// assemble the damping variables
582	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
582	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
583		erfcVal = res.first;
584		derfcVal = res.second;
585	+
586	+	//erfcVal = erfc(dampingAlpha_ * *(idat.rij));
587	+	//derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
588	+
589		c1 = erfcVal * riji;
590		c2 = (-derfcVal + c1) * riji;
591		} else {
#	Line 477 \| Line 593 \| namespace OpenMD {
593		c2 = c1 * riji;
594		}
595
596	<	preVal = idat.electroMult * pre11_ * q_i * q_j;
596	>	preVal = (idat.electroMult) pre11_;
597
598	<	if (summationMethod_ == SHIFTED_POTENTIAL) {
598	>	if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
599		vterm = preVal * (c1 - c1c_);
600	<	dudr = -idat.sw * preVal * c2;
600	>	dudr = - (idat.sw) preVal * c2;
601
602	<	} else if (summationMethod_ == SHIFTED_FORCE) {
603	<	vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) );
604	<	dudr = idat.sw * preVal * (c2c_ - c2);
602	>	} else if (summationMethod_ == esm_SHIFTED_FORCE) {
603	>	vterm = preVal * ( c1 - c1c_ + c2c_( (idat.rij) - cutoffRadius_) );
604	>	dudr = (idat.sw) preVal * (c2c_ - c2);
605
606	<	} else if (summationMethod_ == REACTION_FIELD) {
607	<	rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
606	>	} else if (summationMethod_ == esm_REACTION_FIELD) {
607	>	rfVal = preRF_ * (idat.rij) *(idat.rij);
608	>
609		vterm = preVal * ( riji + rfVal );
610	<	dudr = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
610	>	dudr = (idat.sw) preVal * ( 2.0 * rfVal - riji ) * riji;
611	>
612	>	// if this is an excluded pair, there are still indirect
613	>	// interactions via the reaction field we must worry about:
614
615	+	if (idat.excluded) {
616	+	indirect_vpair += preVal * rfVal;
617	+	indirect_Pot += (idat.sw) preVal * rfVal;
618	+	indirect_dVdr += (idat.sw) preVal * two * rfVal * riji * rhat;
619	+	}
620	+
621		} else {
496	–	vterm = preVal * riji * erfcVal;
622
623	<	dudr = - idat.sw * preVal * c2;
623	>	vterm = preVal * riji * erfcVal;
624	>	dudr = - (idat.sw) preVal * c2;
625	>
626	>	}
627	>
628	>	vpair += vterm * q_i * q_j;
629	>	epot += (idat.sw) vterm * q_i * q_j;
630	>	dVdr += dudr * rhat * q_i * q_j;
631
632	+	if (i_is_Fluctuating) {
633	+	if (idat.excluded) {
634	+	// vFluc1 is the difference between the direct coulomb integral
635	+	// and the normal 1/r-like interaction between point charges.
636	+	coulInt = J1->getValueAt( *(idat.rij) );
637	+	vFluc1 = coulInt - ((idat.sw) vterm);
638	+	} else {
639	+	vFluc1 = 0.0;
640	+	}
641	+	(idat.dVdFQ1) += ( (idat.sw) * vterm + vFluc1 ) * q_j;
642		}
501	–
502	–	idat.vpair += vterm;
503	–	epot += idat.sw * vterm;
643
644	<	dVdr += dudr * rhat;
644	>	if (j_is_Fluctuating) {
645	>	if (idat.excluded) {
646	>	// vFluc2 is the difference between the direct coulomb integral
647	>	// and the normal 1/r-like interaction between point charges.
648	>	coulInt = J2->getValueAt( *(idat.rij) );
649	>	vFluc2 = coulInt - ((idat.sw) vterm);
650	>	} else {
651	>	vFluc2 = 0.0;
652	>	}
653	>	(idat.dVdFQ2) += ( (idat.sw) * vterm + vFluc2 ) * q_i;
654	>	}
655	>
656	>
657		}
658
659		if (j_is_Dipole) {
660		// pref is used by all the possible methods
661	<	pref = idat.electroMult * pre12_ * q_i * mu_j;
662	<	preSw = idat.sw * pref;
661	>	pref = (idat.electroMult) pre12_ * q_i * mu_j;
662	>	preSw = (idat.sw) pref;
663
664	<	if (summationMethod_ == REACTION_FIELD) {
664	>	if (summationMethod_ == esm_REACTION_FIELD) {
665		ri2 = riji * riji;
666		ri3 = ri2 * riji;
667
668	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
669	<	idat.vpair += vterm;
670	<	epot += idat.sw * vterm;
668	>	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
669	>	vpair += vterm;
670	>	epot += (idat.sw) vterm;
671
672	<	dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
673	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);
672	>	dVdr += -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
673	>	duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
674
675	+	// Even if we excluded this pair from direct interactions,
676	+	// we still have the reaction-field-mediated charge-dipole
677	+	// interaction:
678	+
679	+	if (idat.excluded) {
680	+	indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
681	+	indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
682	+	indirect_dVdr += preSw * preRF2_ * uz_j;
683	+	indirect_duduz_j += preSw * rhat * preRF2_ * *(idat.rij);
684	+	}
685	+
686		} else {
687		// determine the inverse r used if we have split dipoles
688		if (j_is_SplitDipole) {
689	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
689	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
690		ri = 1.0 / BigR;
691	<	scale = idat.rij * ri;
691	>	scale = (idat.rij) ri;
692		} else {
693		ri = riji;
694		scale = 1.0;
#	Line 536 \| Line 698 \| namespace OpenMD {
698
699		if (screeningMethod_ == DAMPED) {
700		// assemble the damping variables
701	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
702	<	erfcVal = res.first;
703	<	derfcVal = res.second;
701	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
702	>	//erfcVal = res.first;
703	>	//derfcVal = res.second;
704	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
705	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
706		c1 = erfcVal * ri;
707		c2 = (-derfcVal + c1) * ri;
708		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 553 \| Line 717 \| namespace OpenMD {
717		// calculate the potential
718		pot_term = scale * c2;
719		vterm = -pref * ct_j * pot_term;
720	<	idat.vpair += vterm;
721	<	epot += idat.sw * vterm;
720	>	vpair += vterm;
721	>	epot += (idat.sw) vterm;
722
723		// calculate derivatives for forces and torques
724
#	Line 562 \| Line 726 \| namespace OpenMD {
726		duduz_j += -preSw * pot_term * rhat;
727
728		}
729	+	if (i_is_Fluctuating) {
730	+	(idat.dVdFQ1) += ( (idat.sw) * vterm ) / q_i;
731	+	}
732		}
733
734		if (j_is_Quadrupole) {
#	Line 569 \| Line 736 \| namespace OpenMD {
736		cx2 = cx_j * cx_j;
737		cy2 = cy_j * cy_j;
738		cz2 = cz_j * cz_j;
739	<	pref = idat.electroMult * pre14_ * q_i * one_third_;
739	>	pref = (idat.electroMult) pre14_ * q_i * one_third_;
740
741		if (screeningMethod_ == DAMPED) {
742		// assemble the damping variables
743	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
744	<	erfcVal = res.first;
745	<	derfcVal = res.second;
743	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
744	>	//erfcVal = res.first;
745	>	//derfcVal = res.second;
746	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
747	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
748		c1 = erfcVal * riji;
749		c2 = (-derfcVal + c1) * riji;
750		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 588 \| Line 757 \| namespace OpenMD {
757		}
758
759		// precompute variables for convenience
760	<	preSw = idat.sw * pref;
760	>	preSw = (idat.sw) pref;
761		c2ri = c2 * riji;
762		c3ri = c3 * riji;
763	<	c4rij = c4 * idat.rij;
764	<	rhatdot2 = 2.0 * rhat * c3;
763	>	c4rij = c4 * *(idat.rij) ;
764	>	rhatdot2 = two * rhat * c3;
765		rhatc4 = rhat * c4rij;
766
767		// calculate the potential
#	Line 600 \| Line 769 \| namespace OpenMD {
769		qyy_j * (cy2*c3 - c2ri) +
770		qzz_j * (cz2*c3 - c2ri) );
771		vterm = pref * pot_term;
772	<	idat.vpair += vterm;
773	<	epot += idat.sw * vterm;
772	>	vpair += vterm;
773	>	epot += (idat.sw) vterm;
774
775		// calculate derivatives for the forces and torques
776
777	<	dVdr += -preSw * ( qxx_j* (cx2rhatc4 - (2.0cx_jux_j + rhat)c3ri) +
778	<	qyy_j* (cy2rhatc4 - (2.0cy_juy_j + rhat)c3ri) +
779	<	qzz_j* (cz2rhatc4 - (2.0cz_juz_j + rhat)c3ri));
777	>	dVdr += -preSw * ( qxx_j* (cx2rhatc4 - (twocx_jux_j + rhat)c3ri) +
778	>	qyy_j* (cy2rhatc4 - (twocy_juy_j + rhat)c3ri) +
779	>	qzz_j* (cz2rhatc4 - (twocz_juz_j + rhat)c3ri));
780
781		dudux_j += preSw * qxx_j * cx_j * rhatdot2;
782		duduy_j += preSw * qyy_j * cy_j * rhatdot2;
783		duduz_j += preSw * qzz_j * cz_j * rhatdot2;
784	+	if (i_is_Fluctuating) {
785	+	(idat.dVdFQ1) += ( (idat.sw) * vterm ) / q_i;
786	+	}
787	+
788		}
789		}
790
#	Line 619 \| Line 792 \| namespace OpenMD {
792
793		if (j_is_Charge) {
794		// variables used by all the methods
795	<	pref = idat.electroMult * pre12_ * q_j * mu_i;
796	<	preSw = idat.sw * pref;
795	>	pref = (idat.electroMult) pre12_ * q_j * mu_i;
796	>	preSw = (idat.sw) pref;
797
798	<	if (summationMethod_ == REACTION_FIELD) {
798	>	if (summationMethod_ == esm_REACTION_FIELD) {
799
800		ri2 = riji * riji;
801		ri3 = ri2 * riji;
802
803	<	vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
804	<	idat.vpair += vterm;
805	<	epot += idat.sw * vterm;
803	>	vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
804	>	vpair += vterm;
805	>	epot += (idat.sw) vterm;
806
807	<	dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
807	>	dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
808
809	<	duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
809	>	duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
810	>
811	>	// Even if we excluded this pair from direct interactions,
812	>	// we still have the reaction-field-mediated charge-dipole
813	>	// interaction:
814	>
815	>	if (idat.excluded) {
816	>	indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
817	>	indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
818	>	indirect_dVdr += -preSw * preRF2_ * uz_i;
819	>	indirect_duduz_i += -preSw * rhat * preRF2_ * *(idat.rij);
820	>	}
821
822		} else {
823
824		// determine inverse r if we are using split dipoles
825		if (i_is_SplitDipole) {
826	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
826	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
827		ri = 1.0 / BigR;
828	<	scale = idat.rij * ri;
828	>	scale = (idat.rij) ri;
829		} else {
830		ri = riji;
831		scale = 1.0;
#	Line 651 \| Line 835 \| namespace OpenMD {
835
836		if (screeningMethod_ == DAMPED) {
837		// assemble the damping variables
838	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
839	<	erfcVal = res.first;
840	<	derfcVal = res.second;
838	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
839	>	//erfcVal = res.first;
840	>	//derfcVal = res.second;
841	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
842	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
843		c1 = erfcVal * ri;
844		c2 = (-derfcVal + c1) * ri;
845		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 668 \| Line 854 \| namespace OpenMD {
854		// calculate the potential
855		pot_term = c2 * scale;
856		vterm = pref * ct_i * pot_term;
857	<	idat.vpair += vterm;
858	<	epot += idat.sw * vterm;
857	>	vpair += vterm;
858	>	epot += (idat.sw) vterm;
859
860		// calculate derivatives for the forces and torques
861		dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
862		duduz_i += preSw * pot_term * rhat;
863		}
864	+
865	+	if (j_is_Fluctuating) {
866	+	(idat.dVdFQ2) += ( (idat.sw) * vterm ) / q_j;
867	+	}
868	+
869		}
870
871		if (j_is_Dipole) {
872		// variables used by all methods
873		ct_ij = dot(uz_i, uz_j);
874
875	<	pref = idat.electroMult * pre22_ * mu_i * mu_j;
876	<	preSw = idat.sw * pref;
875	>	pref = (idat.electroMult) pre22_ * mu_i * mu_j;
876	>	preSw = (idat.sw) pref;
877
878	<	if (summationMethod_ == REACTION_FIELD) {
878	>	if (summationMethod_ == esm_REACTION_FIELD) {
879		ri2 = riji * riji;
880		ri3 = ri2 * riji;
881		ri4 = ri2 * ri2;
882
883		vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
884		preRF2_ * ct_ij );
885	<	idat.vpair += vterm;
886	<	epot += idat.sw * vterm;
885	>	vpair += vterm;
886	>	epot += (idat.sw) vterm;
887
888		a1 = 5.0 * ct_i * ct_j - ct_ij;
889
890	<	dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
890	>	dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
891
892	<	duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
893	<	duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
892	>	duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
893	>	duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
894
895	+	if (idat.excluded) {
896	+	indirect_vpair += - pref * preRF2_ * ct_ij;
897	+	indirect_Pot += - preSw * preRF2_ * ct_ij;
898	+	indirect_duduz_i += -preSw * preRF2_ * uz_j;
899	+	indirect_duduz_j += -preSw * preRF2_ * uz_i;
900	+	}
901	+
902		} else {
903
904		if (i_is_SplitDipole) {
905		if (j_is_SplitDipole) {
906	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
906	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i + 0.25 * d_j * d_j);
907		} else {
908	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
908	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
909		}
910		ri = 1.0 / BigR;
911	<	scale = idat.rij * ri;
911	>	scale = (idat.rij) ri;
912		} else {
913		if (j_is_SplitDipole) {
914	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
914	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
915		ri = 1.0 / BigR;
916	<	scale = idat.rij * ri;
916	>	scale = (idat.rij) ri;
917		} else {
918		ri = riji;
919		scale = 1.0;
#	Line 723 \| Line 921 \| namespace OpenMD {
921		}
922		if (screeningMethod_ == DAMPED) {
923		// assemble damping variables
924	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
925	<	erfcVal = res.first;
926	<	derfcVal = res.second;
924	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
925	>	//erfcVal = res.first;
926	>	//derfcVal = res.second;
927	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
928	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
929		c1 = erfcVal * ri;
930		c2 = (-derfcVal + c1) * ri;
931		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 745 \| Line 945 \| namespace OpenMD {
945		preSwSc = preSw * scale;
946		c2ri = c2 * ri;
947		c3ri = c3 * ri;
948	<	c4rij = c4 * idat.rij;
948	>	c4rij = c4 * *(idat.rij) ;
949
950		// calculate the potential
951		pot_term = (ct_ij * c2ri - ctidotj * c3);
952		vterm = pref * pot_term;
953	<	idat.vpair += vterm;
954	<	epot += idat.sw * vterm;
953	>	vpair += vterm;
954	>	epot += (idat.sw) vterm;
955
956		// calculate derivatives for the forces and torques
957		dVdr += preSwSc * ( ctidotj * rhat * c4rij -
#	Line 770 \| Line 970 \| namespace OpenMD {
970		cy2 = cy_i * cy_i;
971		cz2 = cz_i * cz_i;
972
973	<	pref = idat.electroMult * pre14_ * q_j * one_third_;
973	>	pref = (idat.electroMult) pre14_ * q_j * one_third_;
974
975		if (screeningMethod_ == DAMPED) {
976		// assemble the damping variables
977	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
978	<	erfcVal = res.first;
979	<	derfcVal = res.second;
977	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
978	>	//erfcVal = res.first;
979	>	//derfcVal = res.second;
980	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
981	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
982		c1 = erfcVal * riji;
983		c2 = (-derfcVal + c1) * riji;
984		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 789 \| Line 991 \| namespace OpenMD {
991		}
992
993		// precompute some variables for convenience
994	<	preSw = idat.sw * pref;
994	>	preSw = (idat.sw) pref;
995		c2ri = c2 * riji;
996		c3ri = c3 * riji;
997	<	c4rij = c4 * idat.rij;
998	<	rhatdot2 = 2.0 * rhat * c3;
997	>	c4rij = c4 * *(idat.rij) ;
998	>	rhatdot2 = two * rhat * c3;
999		rhatc4 = rhat * c4rij;
1000
1001		// calculate the potential
#	Line 802 \| Line 1004 \| namespace OpenMD {
1004		qzz_i * (cz2 * c3 - c2ri) );
1005
1006		vterm = pref * pot_term;
1007	<	idat.vpair += vterm;
1008	<	epot += idat.sw * vterm;
1007	>	vpair += vterm;
1008	>	epot += (idat.sw) vterm;
1009
1010		// calculate the derivatives for the forces and torques
1011
1012	<	dVdr += -preSw * (qxx_i* (cx2rhatc4 - (2.0cx_iux_i + rhat)c3ri) +
1013	<	qyy_i* (cy2rhatc4 - (2.0cy_iuy_i + rhat)c3ri) +
1014	<	qzz_i* (cz2rhatc4 - (2.0cz_iuz_i + rhat)c3ri));
1012	>	dVdr += -preSw * (qxx_i* (cx2rhatc4 - (twocx_iux_i + rhat)c3ri) +
1013	>	qyy_i* (cy2rhatc4 - (twocy_iuy_i + rhat)c3ri) +
1014	>	qzz_i* (cz2rhatc4 - (twocz_iuz_i + rhat)c3ri));
1015
1016		dudux_i += preSw * qxx_i * cx_i * rhatdot2;
1017		duduy_i += preSw * qyy_i * cy_i * rhatdot2;
1018		duduz_i += preSw * qzz_i * cz_i * rhatdot2;
817	–	}
818	–	}
1019
1020	<	idat.pot += epot;
1021	<	idat.f1 += dVdr;
1020	>	if (j_is_Fluctuating) {
1021	>	(idat.dVdFQ2) += ( (idat.sw) * vterm ) / q_j;
1022	>	}
1023
1024	<	if (i_is_Dipole \|\| i_is_Quadrupole)
824	<	idat.t1 -= cross(uz_i, duduz_i);
825	<	if (i_is_Quadrupole) {
826	<	idat.t1 -= cross(ux_i, dudux_i);
827	<	idat.t1 -= cross(uy_i, duduy_i);
1024	>	}
1025		}
1026
830	–	if (j_is_Dipole \|\| j_is_Quadrupole)
831	–	idat.t2 -= cross(uz_j, duduz_j);
832	–	if (j_is_Quadrupole) {
833	–	idat.t2 -= cross(uz_j, dudux_j);
834	–	idat.t2 -= cross(uz_j, duduy_j);
835	–	}
1027
1028	<	return;
1029	<	}
1028	>	if (!idat.excluded) {
1029	>	*(idat.vpair) += vpair;
1030	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
1031	>	*(idat.f1) += dVdr;
1032	>
1033	>	if (i_is_Dipole \|\| i_is_Quadrupole)
1034	>	*(idat.t1) -= cross(uz_i, duduz_i);
1035	>	if (i_is_Quadrupole) {
1036	>	*(idat.t1) -= cross(ux_i, dudux_i);
1037	>	*(idat.t1) -= cross(uy_i, duduy_i);
1038	>	}
1039	>
1040	>	if (j_is_Dipole \|\| j_is_Quadrupole)
1041	>	*(idat.t2) -= cross(uz_j, duduz_j);
1042	>	if (j_is_Quadrupole) {
1043	>	*(idat.t2) -= cross(uz_j, dudux_j);
1044	>	*(idat.t2) -= cross(uz_j, duduy_j);
1045	>	}
1046
1047	<	void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
1047	>	} else {
1048
1049	<	if (!initialized_) initialize();
1050	<
844	<	ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
845	<	ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
846	<
847	<	// logicals
1049	>	// only accumulate the forces and torques resulting from the
1050	>	// indirect reaction field terms.
1051
1052	<	bool i_is_Charge = data1.is_Charge;
850	<	bool i_is_Dipole = data1.is_Dipole;
851	<
852	<	bool j_is_Charge = data2.is_Charge;
853	<	bool j_is_Dipole = data2.is_Dipole;
854	<
855	<	RealType q_i, q_j;
856	<
857	<	// The skippedCharge computation is needed by the real-space cutoff methods
858	<	// (i.e. shifted force and shifted potential)
859	<
860	<	if (i_is_Charge) {
861	<	q_i = data1.charge;
862	<	skdat.skippedCharge2 += q_i;
863	<	}
864	<
865	<	if (j_is_Charge) {
866	<	q_j = data2.charge;
867	<	skdat.skippedCharge1 += q_j;
868	<	}
869	<
870	<	// the rest of this function should only be necessary for reaction field.
871	<
872	<	if (summationMethod_ == REACTION_FIELD) {
873	<	RealType riji, ri2, ri3;
874	<	RealType q_i, mu_i, ct_i;
875	<	RealType q_j, mu_j, ct_j;
876	<	RealType preVal, rfVal, vterm, dudr, pref, myPot;
877	<	Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
878	<
879	<	// some variables we'll need independent of electrostatic type:
1052	>	*(idat.vpair) += indirect_vpair;
1053
1054	<	riji = 1.0 / skdat.rij;
1055	<	rhat = skdat.d * riji;
1056	<
1057	<	if (i_is_Dipole) {
885	<	mu_i = data1.dipole_moment;
886	<	uz_i = skdat.eFrame1.getColumn(2);
887	<	ct_i = dot(uz_i, rhat);
888	<	duduz_i = V3Zero;
889	<	}
890	<
891	<	if (j_is_Dipole) {
892	<	mu_j = data2.dipole_moment;
893	<	uz_j = skdat.eFrame2.getColumn(2);
894	<	ct_j = dot(uz_j, rhat);
895	<	duduz_j = V3Zero;
896	<	}
897	<
898	<	if (i_is_Charge) {
899	<	if (j_is_Charge) {
900	<	preVal = skdat.electroMult * pre11_ * q_i * q_j;
901	<	rfVal = preRF_ * skdat.rij * skdat.rij;
902	<	vterm = preVal * rfVal;
903	<	myPot += skdat.sw * vterm;
904	<	dudr = skdat.sw * preVal * 2.0 * rfVal * riji;
905	<	dVdr += dudr * rhat;
906	<	}
907	<
908	<	if (j_is_Dipole) {
909	<	ri2 = riji * riji;
910	<	ri3 = ri2 * riji;
911	<	pref = skdat.electroMult * pre12_ * q_i * mu_j;
912	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
913	<	myPot += skdat.sw * vterm;
914	<	dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
915	<	duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
916	<	}
917	<	}
918	<	if (i_is_Dipole) {
919	<	if (j_is_Charge) {
920	<	ri2 = riji * riji;
921	<	ri3 = ri2 * riji;
922	<	pref = skdat.electroMult * pre12_ * q_j * mu_i;
923	<	vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
924	<	myPot += skdat.sw * vterm;
925	<	dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
926	<	duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
927	<	}
928	<	}
1054	>	((idat.excludedPot))[ELECTROSTATIC_FAMILY] += ((idat.sw) * vterm +
1055	>	vFluc1 ) * q_i * q_j;
1056	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1057	>	*(idat.f1) += indirect_dVdr;
1058
930	–	// accumulate the forces and torques resulting from the self term
931	–	skdat.pot += myPot;
932	–	skdat.f1 += dVdr;
933	–
1059		if (i_is_Dipole)
1060	<	skdat.t1 -= cross(uz_i, duduz_i);
1060	>	*(idat.t1) -= cross(uz_i, indirect_duduz_i);
1061		if (j_is_Dipole)
1062	<	skdat.t2 -= cross(uz_j, duduz_j);
1062	>	*(idat.t2) -= cross(uz_j, indirect_duduz_j);
1063		}
1064	<	}
1064	>
1065	>	return;
1066	>	}
1067
1068	<	void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1069	<	RealType mu1, preVal, chg1, self;
943	<
1068	>	void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1069	>	RealType mu1, preVal, self;
1070		if (!initialized_) initialize();
1071	<
1072	<	ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1071	>
1072	>	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1073
1074		// logicals
949	–
1075		bool i_is_Charge = data.is_Charge;
1076		bool i_is_Dipole = data.is_Dipole;
1077	+	bool i_is_Fluctuating = data.is_Fluctuating;
1078	+	RealType chg1 = data.fixedCharge;
1079	+
1080	+	if (i_is_Fluctuating) {
1081	+	chg1 += *(sdat.flucQ);
1082	+	// dVdFQ is really a force, so this is negative the derivative
1083	+	(sdat.dVdFQ) -= (sdat.flucQ) * data.hardness + data.electronegativity;
1084	+	((sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (sdat.flucQ) *
1085	+	((sdat.flucQ) data.hardness * 0.5 + data.electronegativity);
1086	+	}
1087
1088	<	if (summationMethod_ == REACTION_FIELD) {
1088	>	if (summationMethod_ == esm_REACTION_FIELD) {
1089		if (i_is_Dipole) {
1090		mu1 = data.dipole_moment;
1091		preVal = pre22_ * preRF2_ * mu1 * mu1;
1092	<	scdat.pot -= 0.5 * preVal;
1092	>	((sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 preVal;
1093
1094		// The self-correction term adds into the reaction field vector
1095	<	Vector3d uz_i = scdat.eFrame.getColumn(2);
1095	>	Vector3d uz_i = sdat.eFrame->getColumn(2);
1096		Vector3d ei = preVal * uz_i;
1097
1098		// This looks very wrong. A vector crossed with itself is zero.
1099	<	scdat.t -= cross(uz_i, ei);
1099	>	*(sdat.t) -= cross(uz_i, ei);
1100		}
1101	<	} else if (summationMethod_ == SHIFTED_FORCE \|\| summationMethod_ == SHIFTED_POTENTIAL) {
1101	>	} else if (summationMethod_ == esm_SHIFTED_FORCE \|\| summationMethod_ == esm_SHIFTED_POTENTIAL) {
1102		if (i_is_Charge) {
968	–	chg1 = data.charge;
1103		if (screeningMethod_ == DAMPED) {
1104	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1104	>	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1105		} else {
1106	<	self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1106	>	self = - 0.5 * rcuti_ * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1107		}
1108	<	scdat.pot += self;
1108	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1109		}
1110		}
1111		}
1112	+
1113	+	RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
1114	+	// This seems to work moderately well as a default. There's no
1115	+	// inherent scale for 1/r interactions that we can standardize.
1116	+	// 12 angstroms seems to be a reasonably good guess for most
1117	+	// cases.
1118	+	return 12.0;
1119	+	}
1120		}

Diff Legend

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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents): Revision 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs. Revision 1767 by gezelter, Fri Jul 6 22:01:58 2012 UTC

Diff Legend

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs.
Revision 1767 by gezelter, Fri Jul 6 22:01:58 2012 UTC