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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs.
Revision 1718 by gezelter, Thu May 24 01:29:59 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/MultipoleAdapter.hpp"
52		#include "io/Globals.hpp"
53	+	#include "nonbonded/SlaterIntegrals.hpp"
54	+	#include "utils/PhysicalConstants.hpp"
55
56	+
57		namespace OpenMD {
58
59		Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
60	<	forceField_(NULL) {}
60	>	forceField_(NULL), info_(NULL),
61	>	haveCutoffRadius_(false),
62	>	haveDampingAlpha_(false),
63	>	haveDielectric_(false),
64	>	haveElectroSpline_(false)
65	>	{}
66
67		void Electrostatic::initialize() {
68	<
69	<	Globals* simParams_;
68	>
69	>	Globals* simParams_ = info_->getSimParams();
70
71		summationMap_["HARD"] = esm_HARD;
72	+	summationMap_["NONE"] = esm_HARD;
73		summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
74		summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
75		summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
#	Line 97 \| Line 108 \| namespace OpenMD {
108		screeningMethod_ = UNDAMPED;
109		dielectric_ = 1.0;
110		one_third_ = 1.0 / 3.0;
100	–	haveCutoffRadius_ = false;
101	–	haveDampingAlpha_ = false;
102	–	haveDielectric_ = false;
103	–	haveElectroSpline_ = false;
111
112		// check the summation method:
113		if (simParams_->haveElectrostaticSummationMethod()) {
#	Line 115 \| Line 122 \| namespace OpenMD {
122		sprintf( painCave.errMsg,
123		"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
124		"\t(Input file specified %s .)\n"
125	<	"\telectrostaticSummationMethod must be one of: \"none\",\n"
125	>	"\telectrostaticSummationMethod must be one of: \"hard\",\n"
126		"\t\"shifted_potential\", \"shifted_force\", or \n"
127		"\t\"reaction_field\".\n", myMethod.c_str() );
128		painCave.isFatal = 1;
#	Line 248 \| Line 255 \| namespace OpenMD {
255		preRF2_ = 2.0 * preRF_;
256		}
257
258	<	RealType dx = cutoffRadius_ / RealType(np_ - 1);
258	>	// Add a 2 angstrom safety window to deal with cutoffGroups that
259	>	// have charged atoms longer than the cutoffRadius away from each
260	>	// other. Splining may not be the best choice here. Direct calls
261	>	// to erfc might be preferrable.
262	>
263	>	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
264		RealType rval;
265		vector<RealType> rvals;
266		vector<RealType> yvals;
#	Line 272 \| Line 284 \| namespace OpenMD {
284		electrostaticAtomData.is_SplitDipole = false;
285		electrostaticAtomData.is_Quadrupole = false;
286
287	<	if (atomType->isCharge()) {
276	<	GenericData* data = atomType->getPropertyByName("Charge");
287	>	FixedChargeAdapter fca = FixedChargeAdapter(atomType);
288
289	<	if (data == NULL) {
279	<	sprintf( painCave.errMsg, "Electrostatic::addType could not find "
280	<	"Charge\n"
281	<	"\tparameters for atomType %s.\n",
282	<	atomType->getName().c_str());
283	<	painCave.severity = OPENMD_ERROR;
284	<	painCave.isFatal = 1;
285	<	simError();
286	<	}
287	<
288	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
289	<	if (doubleData == NULL) {
290	<	sprintf( painCave.errMsg,
291	<	"Electrostatic::addType could not convert GenericData to "
292	<	"Charge for\n"
293	<	"\tatom type %s\n", atomType->getName().c_str());
294	<	painCave.severity = OPENMD_ERROR;
295	<	painCave.isFatal = 1;
296	<	simError();
297	<	}
289	>	if (fca.isFixedCharge()) {
290		electrostaticAtomData.is_Charge = true;
291	<	electrostaticAtomData.charge = doubleData->getData();
291	>	electrostaticAtomData.charge = fca.getCharge();
292		}
293
294	<	if (atomType->isDirectional()) {
295	<	DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
296	<
305	<	if (daType->isDipole()) {
306	<	GenericData* data = daType->getPropertyByName("Dipole");
307	<
308	<	if (data == NULL) {
309	<	sprintf( painCave.errMsg,
310	<	"Electrostatic::addType could not find Dipole\n"
311	<	"\tparameters for atomType %s.\n",
312	<	daType->getName().c_str());
313	<	painCave.severity = OPENMD_ERROR;
314	<	painCave.isFatal = 1;
315	<	simError();
316	<	}
317	<
318	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
319	<	if (doubleData == NULL) {
320	<	sprintf( painCave.errMsg,
321	<	"Electrostatic::addType could not convert GenericData to "
322	<	"Dipole Moment\n"
323	<	"\tfor atom type %s\n", daType->getName().c_str());
324	<	painCave.severity = OPENMD_ERROR;
325	<	painCave.isFatal = 1;
326	<	simError();
327	<	}
294	>	MultipoleAdapter ma = MultipoleAdapter(atomType);
295	>	if (ma.isMultipole()) {
296	>	if (ma.isDipole()) {
297		electrostaticAtomData.is_Dipole = true;
298	<	electrostaticAtomData.dipole_moment = doubleData->getData();
298	>	electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
299		}
300	<
332	<	if (daType->isSplitDipole()) {
333	<	GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
334	<
335	<	if (data == NULL) {
336	<	sprintf(painCave.errMsg,
337	<	"Electrostatic::addType could not find SplitDipoleDistance\n"
338	<	"\tparameter for atomType %s.\n",
339	<	daType->getName().c_str());
340	<	painCave.severity = OPENMD_ERROR;
341	<	painCave.isFatal = 1;
342	<	simError();
343	<	}
344	<
345	<	DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
346	<	if (doubleData == NULL) {
347	<	sprintf( painCave.errMsg,
348	<	"Electrostatic::addType could not convert GenericData to "
349	<	"SplitDipoleDistance for\n"
350	<	"\tatom type %s\n", daType->getName().c_str());
351	<	painCave.severity = OPENMD_ERROR;
352	<	painCave.isFatal = 1;
353	<	simError();
354	<	}
300	>	if (ma.isSplitDipole()) {
301		electrostaticAtomData.is_SplitDipole = true;
302	<	electrostaticAtomData.split_dipole_distance = doubleData->getData();
302	>	electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
303		}
304	<
359	<	if (daType->isQuadrupole()) {
360	<	GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
361	<
362	<	if (data == NULL) {
363	<	sprintf( painCave.errMsg,
364	<	"Electrostatic::addType could not find QuadrupoleMoments\n"
365	<	"\tparameter for atomType %s.\n",
366	<	daType->getName().c_str());
367	<	painCave.severity = OPENMD_ERROR;
368	<	painCave.isFatal = 1;
369	<	simError();
370	<	}
371	<
304	>	if (ma.isQuadrupole()) {
305		// Quadrupoles in OpenMD are set as the diagonal elements
306		// of the diagonalized traceless quadrupole moment tensor.
307		// The column vectors of the unitary matrix that diagonalizes
308		// the quadrupole moment tensor become the eFrame (or the
309		// electrostatic version of the body-fixed frame.
377	–
378	–	Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
379	–	if (v3dData == NULL) {
380	–	sprintf( painCave.errMsg,
381	–	"Electrostatic::addType could not convert GenericData to "
382	–	"Quadrupole Moments for\n"
383	–	"\tatom type %s\n", daType->getName().c_str());
384	–	painCave.severity = OPENMD_ERROR;
385	–	painCave.isFatal = 1;
386	–	simError();
387	–	}
310		electrostaticAtomData.is_Quadrupole = true;
311	<	electrostaticAtomData.quadrupole_moments = v3dData->getData();
311	>	electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
312		}
313		}
314
315	<	AtomTypeProperties atp = atomType->getATP();
315	>	FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
316
317	+	if (fqa.isFluctuatingCharge()) {
318	+	electrostaticAtomData.is_FluctuatingCharge = true;
319	+	electrostaticAtomData.electronegativity = fca.getElectronegativity();
320	+	electrostaticAtomData.hardness = fca.getHardness();
321	+	electrostaticAtomData.slaterN = fca.getSlaterN();
322	+	electrostaticAtomData.slaterZeta = fca.getSlaterZeta();
323	+	}
324	+
325		pair<map<int,AtomType*>::iterator,bool> ret;
326	<	ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
326	>	ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
327	>	atomType) );
328		if (ret.second == false) {
329		sprintf( painCave.errMsg,
330		"Electrostatic already had a previous entry with ident %d\n",
331	<	atp.ident);
331	>	atomType->getIdent() );
332		painCave.severity = OPENMD_INFO;
333		painCave.isFatal = 0;
334		simError();
335		}
336
337	<	ElectrostaticMap[atomType] = electrostaticAtomData;
337	>	ElectrostaticMap[atomType] = electrostaticAtomData;
338	>
339	>	// Now, iterate over all known types and add to the mixing map:
340	>
341	>	map<AtomType*, ElectrostaticAtomData>::iterator it;
342	>	for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
343	>	AtomType* atype2 = (*it).first;
344	>
345	>	if ((*it).is_FluctuatingCharge && electrostaticAtomData.is_FluctuatingCharge) {
346	>
347	>	RealType a = electrostaticAtomData.slaterZeta;
348	>	RealType b = (*it).slaterZeta;
349	>	int m = electrostaticAtomData.slaterN;
350	>	int n = (*it).slaterN;
351	>
352	>	// Create the spline of the coulombic integral for s-type
353	>	// Slater orbitals. Add a 2 angstrom safety window to deal
354	>	// with cutoffGroups that have charged atoms longer than the
355	>	// cutoffRadius away from each other.
356	>
357	>	RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
358	>	vector<RealType> rvals;
359	>	vector<RealType> J1vals;
360	>	vector<RealType> J2vals;
361	>	for (int i = 0; i < np_; i++) {
362	>	rval = RealType(i) * dr;
363	>	rvals.push_back(rval);
364	>	J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
365	>	J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
366	>	}
367	>
368	>	CubicSpline J1 = new CubicSpline();
369	>	J1->addPoints(rvals, J1vals);
370	>	CubicSpline J2 = new CubicSpline();
371	>	J2->addPoints(rvals, J2vals);
372	>
373	>	pair<AtomType, AtomType> key1, key2;
374	>	key1 = make_pair(atomType, atype2);
375	>	key2 = make_pair(atype2, atomType);
376	>
377	>	Jij[key1] = J1;
378	>	Jij[key2] = J2;
379	>	}
380	>	}
381	>
382		return;
383		}
384
385	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
386	<	RealType theRSW ) {
412	<	cutoffRadius_ = theECR;
385	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
386	>	cutoffRadius_ = rCut;
387		rrf_ = cutoffRadius_;
414	–	rt_ = theRSW;
388		haveCutoffRadius_ = true;
389		}
390	+
391	+	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
392	+	rt_ = rSwitch;
393	+	}
394		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
395		summationMethod_ = esm;
396		}
#	Line 443 \| Line 420 \| namespace OpenMD {
420		RealType ct_i, ct_j, ct_ij, a1;
421		RealType riji, ri, ri2, ri3, ri4;
422		RealType pref, vterm, epot, dudr;
423	+	RealType vpair(0.0);
424		RealType scale, sc2;
425		RealType pot_term, preVal, rfVal;
426		RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
427		RealType preSw, preSwSc;
428		RealType c1, c2, c3, c4;
429	<	RealType erfcVal, derfcVal;
429	>	RealType erfcVal(1.0), derfcVal(0.0);
430		RealType BigR;
431	+	RealType two(2.0), three(3.0);
432
433		Vector3d Q_i, Q_j;
434		Vector3d ux_i, uy_i, uz_i;
#	Line 459 \| Line 438 \| namespace OpenMD {
438		Vector3d rhatdot2, rhatc4;
439		Vector3d dVdr;
440
441	+	// variables for indirect (reaction field) interactions for excluded pairs:
442	+	RealType indirect_Pot(0.0);
443	+	RealType indirect_vpair(0.0);
444	+	Vector3d indirect_dVdr(V3Zero);
445	+	Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
446	+
447		pair<RealType, RealType> res;
448
449		if (!initialized_) initialize();
#	Line 483 \| Line 468 \| namespace OpenMD {
468		bool j_is_SplitDipole = data2.is_SplitDipole;
469		bool j_is_Quadrupole = data2.is_Quadrupole;
470
471	<	if (i_is_Charge)
471	>	if (i_is_Charge) {
472		q_i = data1.charge;
473	+	if (idat.excluded) {
474	+	*(idat.skippedCharge2) += q_i;
475	+	}
476	+	}
477
478		if (i_is_Dipole) {
479		mu_i = data1.dipole_moment;
#	Line 517 \| Line 506 \| namespace OpenMD {
506		duduz_i = V3Zero;
507		}
508
509	<	if (j_is_Charge)
509	>	if (j_is_Charge) {
510		q_j = data2.charge;
511	+	if (idat.excluded) {
512	+	*(idat.skippedCharge1) += q_j;
513	+	}
514	+	}
515
516	+
517		if (j_is_Dipole) {
518		mu_j = data2.dipole_moment;
519		uz_j = idat.eFrame2->getColumn(2);
#	Line 559 \| Line 553 \| namespace OpenMD {
553		if (j_is_Charge) {
554		if (screeningMethod_ == DAMPED) {
555		// assemble the damping variables
556	<	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
557	<	erfcVal = res.first;
558	<	derfcVal = res.second;
556	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
557	>	//erfcVal = res.first;
558	>	//derfcVal = res.second;
559	>
560	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
561	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
562	>
563		c1 = erfcVal * riji;
564		c2 = (-derfcVal + c1) * riji;
565		} else {
#	Line 580 \| Line 578 \| namespace OpenMD {
578		dudr = (idat.sw) preVal * (c2c_ - c2);
579
580		} else if (summationMethod_ == esm_REACTION_FIELD) {
581	<	rfVal = (idat.electroMult) preRF_ * (idat.rij) *(idat.rij) ;
581	>	rfVal = preRF_ * (idat.rij) *(idat.rij);
582	>
583		vterm = preVal * ( riji + rfVal );
584		dudr = (idat.sw) preVal * ( 2.0 * rfVal - riji ) * riji;
585	+
586	+	// if this is an excluded pair, there are still indirect
587	+	// interactions via the reaction field we must worry about:
588
589	+	if (idat.excluded) {
590	+	indirect_vpair += preVal * rfVal;
591	+	indirect_Pot += (idat.sw) preVal * rfVal;
592	+	indirect_dVdr += (idat.sw) preVal * two * rfVal * riji * rhat;
593	+	}
594	+
595		} else {
588	–	vterm = preVal * riji * erfcVal;
596
597	+	vterm = preVal * riji * erfcVal;
598		dudr = - (idat.sw) preVal * c2;
599
600		}
593	–
594	–	*(idat.vpair) += vterm;
595	–	epot += (idat.sw) vterm;
601
602	<	dVdr += dudr * rhat;
602	>	vpair += vterm;
603	>	epot += (idat.sw) vterm;
604	>	dVdr += dudr * rhat;
605		}
606
607		if (j_is_Dipole) {
#	Line 607 \| Line 614 \| namespace OpenMD {
614		ri3 = ri2 * riji;
615
616		vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
617	<	*(idat.vpair) += vterm;
617	>	vpair += vterm;
618		epot += (idat.sw) vterm;
619
620	<	dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
620	>	dVdr += -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
621		duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
622
623	+	// Even if we excluded this pair from direct interactions,
624	+	// we still have the reaction-field-mediated charge-dipole
625	+	// interaction:
626	+
627	+	if (idat.excluded) {
628	+	indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
629	+	indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
630	+	indirect_dVdr += preSw * preRF2_ * uz_j;
631	+	indirect_duduz_j += preSw * rhat * preRF2_ * *(idat.rij);
632	+	}
633	+
634		} else {
635		// determine the inverse r used if we have split dipoles
636		if (j_is_SplitDipole) {
#	Line 628 \| Line 646 \| namespace OpenMD {
646
647		if (screeningMethod_ == DAMPED) {
648		// assemble the damping variables
649	<	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
650	<	erfcVal = res.first;
651	<	derfcVal = res.second;
649	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
650	>	//erfcVal = res.first;
651	>	//derfcVal = res.second;
652	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
653	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
654		c1 = erfcVal * ri;
655		c2 = (-derfcVal + c1) * ri;
656		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 645 \| Line 665 \| namespace OpenMD {
665		// calculate the potential
666		pot_term = scale * c2;
667		vterm = -pref * ct_j * pot_term;
668	<	*(idat.vpair) += vterm;
668	>	vpair += vterm;
669		epot += (idat.sw) vterm;
670
671		// calculate derivatives for forces and torques
#	Line 665 \| Line 685 \| namespace OpenMD {
685
686		if (screeningMethod_ == DAMPED) {
687		// assemble the damping variables
688	<	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
689	<	erfcVal = res.first;
690	<	derfcVal = res.second;
688	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
689	>	//erfcVal = res.first;
690	>	//derfcVal = res.second;
691	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
692	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
693		c1 = erfcVal * riji;
694		c2 = (-derfcVal + c1) * riji;
695		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 684 \| Line 706 \| namespace OpenMD {
706		c2ri = c2 * riji;
707		c3ri = c3 * riji;
708		c4rij = c4 * *(idat.rij) ;
709	<	rhatdot2 = 2.0 * rhat * c3;
709	>	rhatdot2 = two * rhat * c3;
710		rhatc4 = rhat * c4rij;
711
712		// calculate the potential
#	Line 692 \| Line 714 \| namespace OpenMD {
714		qyy_j * (cy2*c3 - c2ri) +
715		qzz_j * (cz2*c3 - c2ri) );
716		vterm = pref * pot_term;
717	<	*(idat.vpair) += vterm;
717	>	vpair += vterm;
718		epot += (idat.sw) vterm;
719
720		// calculate derivatives for the forces and torques
721
722	<	dVdr += -preSw * ( qxx_j* (cx2rhatc4 - (2.0cx_jux_j + rhat)c3ri) +
723	<	qyy_j* (cy2rhatc4 - (2.0cy_juy_j + rhat)c3ri) +
724	<	qzz_j* (cz2rhatc4 - (2.0cz_juz_j + rhat)c3ri));
722	>	dVdr += -preSw * ( qxx_j* (cx2rhatc4 - (twocx_jux_j + rhat)c3ri) +
723	>	qyy_j* (cy2rhatc4 - (twocy_juy_j + rhat)c3ri) +
724	>	qzz_j* (cz2rhatc4 - (twocz_juz_j + rhat)c3ri));
725
726		dudux_j += preSw * qxx_j * cx_j * rhatdot2;
727		duduy_j += preSw * qyy_j * cy_j * rhatdot2;
#	Line 720 \| Line 742 \| namespace OpenMD {
742		ri3 = ri2 * riji;
743
744		vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
745	<	*(idat.vpair) += vterm;
745	>	vpair += vterm;
746		epot += (idat.sw) vterm;
747
748	<	dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
748	>	dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
749
750		duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
751	+
752	+	// Even if we excluded this pair from direct interactions,
753	+	// we still have the reaction-field-mediated charge-dipole
754	+	// interaction:
755	+
756	+	if (idat.excluded) {
757	+	indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
758	+	indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
759	+	indirect_dVdr += -preSw * preRF2_ * uz_i;
760	+	indirect_duduz_i += -preSw * rhat * preRF2_ * *(idat.rij);
761	+	}
762
763		} else {
764
#	Line 743 \| Line 776 \| namespace OpenMD {
776
777		if (screeningMethod_ == DAMPED) {
778		// assemble the damping variables
779	<	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
780	<	erfcVal = res.first;
781	<	derfcVal = res.second;
779	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
780	>	//erfcVal = res.first;
781	>	//derfcVal = res.second;
782	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
783	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
784		c1 = erfcVal * ri;
785		c2 = (-derfcVal + c1) * ri;
786		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 760 \| Line 795 \| namespace OpenMD {
795		// calculate the potential
796		pot_term = c2 * scale;
797		vterm = pref * ct_i * pot_term;
798	<	*(idat.vpair) += vterm;
798	>	vpair += vterm;
799		epot += (idat.sw) vterm;
800
801		// calculate derivatives for the forces and torques
#	Line 783 \| Line 818 \| namespace OpenMD {
818
819		vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
820		preRF2_ * ct_ij );
821	<	*(idat.vpair) += vterm;
821	>	vpair += vterm;
822		epot += (idat.sw) vterm;
823
824		a1 = 5.0 * ct_i * ct_j - ct_ij;
825
826	<	dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
826	>	dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
827
828	<	duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
829	<	duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
828	>	duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
829	>	duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
830
831	+	if (idat.excluded) {
832	+	indirect_vpair += - pref * preRF2_ * ct_ij;
833	+	indirect_Pot += - preSw * preRF2_ * ct_ij;
834	+	indirect_duduz_i += -preSw * preRF2_ * uz_j;
835	+	indirect_duduz_j += -preSw * preRF2_ * uz_i;
836	+	}
837	+
838		} else {
839
840		if (i_is_SplitDipole) {
#	Line 815 \| Line 857 \| namespace OpenMD {
857		}
858		if (screeningMethod_ == DAMPED) {
859		// assemble damping variables
860	<	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
861	<	erfcVal = res.first;
862	<	derfcVal = res.second;
860	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
861	>	//erfcVal = res.first;
862	>	//derfcVal = res.second;
863	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
864	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
865		c1 = erfcVal * ri;
866		c2 = (-derfcVal + c1) * ri;
867		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 842 \| Line 886 \| namespace OpenMD {
886		// calculate the potential
887		pot_term = (ct_ij * c2ri - ctidotj * c3);
888		vterm = pref * pot_term;
889	<	*(idat.vpair) += vterm;
889	>	vpair += vterm;
890		epot += (idat.sw) vterm;
891
892		// calculate derivatives for the forces and torques
#	Line 866 \| Line 910 \| namespace OpenMD {
910
911		if (screeningMethod_ == DAMPED) {
912		// assemble the damping variables
913	<	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
914	<	erfcVal = res.first;
915	<	derfcVal = res.second;
913	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
914	>	//erfcVal = res.first;
915	>	//derfcVal = res.second;
916	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
917	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
918		c1 = erfcVal * riji;
919		c2 = (-derfcVal + c1) * riji;
920		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 885 \| Line 931 \| namespace OpenMD {
931		c2ri = c2 * riji;
932		c3ri = c3 * riji;
933		c4rij = c4 * *(idat.rij) ;
934	<	rhatdot2 = 2.0 * rhat * c3;
934	>	rhatdot2 = two * rhat * c3;
935		rhatc4 = rhat * c4rij;
936
937		// calculate the potential
#	Line 894 \| Line 940 \| namespace OpenMD {
940		qzz_i * (cz2 * c3 - c2ri) );
941
942		vterm = pref * pot_term;
943	<	*(idat.vpair) += vterm;
943	>	vpair += vterm;
944		epot += (idat.sw) vterm;
945
946		// calculate the derivatives for the forces and torques
947
948	<	dVdr += -preSw * (qxx_i* (cx2rhatc4 - (2.0cx_iux_i + rhat)c3ri) +
949	<	qyy_i* (cy2rhatc4 - (2.0cy_iuy_i + rhat)c3ri) +
950	<	qzz_i* (cz2rhatc4 - (2.0cz_iuz_i + rhat)c3ri));
948	>	dVdr += -preSw * (qxx_i* (cx2rhatc4 - (twocx_iux_i + rhat)c3ri) +
949	>	qyy_i* (cy2rhatc4 - (twocy_iuy_i + rhat)c3ri) +
950	>	qzz_i* (cz2rhatc4 - (twocz_iuz_i + rhat)c3ri));
951
952		dudux_i += preSw * qxx_i * cx_i * rhatdot2;
953		duduy_i += preSw * qyy_i * cy_i * rhatdot2;
#	Line 909 \| Line 955 \| namespace OpenMD {
955		}
956		}
957
912	–	idat.pot[ELECTROSTATIC_FAMILY] += epot;
913	–	*(idat.f1) += dVdr;
958
959	<	if (i_is_Dipole \|\| i_is_Quadrupole)
960	<	*(idat.t1) -= cross(uz_i, duduz_i);
961	<	if (i_is_Quadrupole) {
962	<	*(idat.t1) -= cross(ux_i, dudux_i);
963	<	*(idat.t1) -= cross(uy_i, duduy_i);
964	<	}
965	<
966	<	if (j_is_Dipole \|\| j_is_Quadrupole)
967	<	*(idat.t2) -= cross(uz_j, duduz_j);
968	<	if (j_is_Quadrupole) {
969	<	*(idat.t2) -= cross(uz_j, dudux_j);
970	<	*(idat.t2) -= cross(uz_j, duduy_j);
971	<	}
959	>	if (!idat.excluded) {
960	>	*(idat.vpair) += vpair;
961	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
962	>	*(idat.f1) += dVdr;
963	>
964	>	if (i_is_Dipole \|\| i_is_Quadrupole)
965	>	*(idat.t1) -= cross(uz_i, duduz_i);
966	>	if (i_is_Quadrupole) {
967	>	*(idat.t1) -= cross(ux_i, dudux_i);
968	>	*(idat.t1) -= cross(uy_i, duduy_i);
969	>	}
970	>
971	>	if (j_is_Dipole \|\| j_is_Quadrupole)
972	>	*(idat.t2) -= cross(uz_j, duduz_j);
973	>	if (j_is_Quadrupole) {
974	>	*(idat.t2) -= cross(uz_j, dudux_j);
975	>	*(idat.t2) -= cross(uz_j, duduy_j);
976	>	}
977
978	<	return;
930	<	}
978	>	} else {
979
980	<	void Electrostatic::calcSkipCorrection(InteractionData &idat) {
980	>	// only accumulate the forces and torques resulting from the
981	>	// indirect reaction field terms.
982
983	<	if (!initialized_) initialize();
984	<
985	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
937	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
938	<
939	<	// logicals
940	<
941	<	bool i_is_Charge = data1.is_Charge;
942	<	bool i_is_Dipole = data1.is_Dipole;
943	<
944	<	bool j_is_Charge = data2.is_Charge;
945	<	bool j_is_Dipole = data2.is_Dipole;
946	<
947	<	RealType q_i, q_j;
948	<
949	<	// The skippedCharge computation is needed by the real-space cutoff methods
950	<	// (i.e. shifted force and shifted potential)
951	<
952	<	if (i_is_Charge) {
953	<	q_i = data1.charge;
954	<	*(idat.skippedCharge2) += q_i;
955	<	}
956	<
957	<	if (j_is_Charge) {
958	<	q_j = data2.charge;
959	<	*(idat.skippedCharge1) += q_j;
960	<	}
961	<
962	<	// the rest of this function should only be necessary for reaction field.
963	<
964	<	if (summationMethod_ == esm_REACTION_FIELD) {
965	<	RealType riji, ri2, ri3;
966	<	RealType mu_i, ct_i;
967	<	RealType mu_j, ct_j;
968	<	RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
969	<	Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
970	<
971	<	// some variables we'll need independent of electrostatic type:
983	>	*(idat.vpair) += indirect_vpair;
984	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
985	>	*(idat.f1) += indirect_dVdr;
986
973	–	riji = 1.0 / *(idat.rij) ;
974	–	rhat = (idat.d) riji;
975	–
976	–	if (i_is_Dipole) {
977	–	mu_i = data1.dipole_moment;
978	–	uz_i = idat.eFrame1->getColumn(2);
979	–	ct_i = dot(uz_i, rhat);
980	–	duduz_i = V3Zero;
981	–	}
982	–
983	–	if (j_is_Dipole) {
984	–	mu_j = data2.dipole_moment;
985	–	uz_j = idat.eFrame2->getColumn(2);
986	–	ct_j = dot(uz_j, rhat);
987	–	duduz_j = V3Zero;
988	–	}
989	–
990	–	if (i_is_Charge) {
991	–	if (j_is_Charge) {
992	–	preVal = (idat.electroMult) pre11_ * q_i * q_j;
993	–	rfVal = preRF_ * (idat.rij) *(idat.rij) ;
994	–	vterm = preVal * rfVal;
995	–	myPot += (idat.sw) vterm;
996	–	dudr = (idat.sw) preVal * 2.0 * rfVal * riji;
997	–	dVdr += dudr * rhat;
998	–	}
999	–
1000	–	if (j_is_Dipole) {
1001	–	ri2 = riji * riji;
1002	–	ri3 = ri2 * riji;
1003	–	pref = (idat.electroMult) pre12_ * q_i * mu_j;
1004	–	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
1005	–	myPot += (idat.sw) vterm;
1006	–	dVdr += - (idat.sw) pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1007	–	duduz_j += - (idat.sw) pref * rhat * (ri2 - preRF2_ * *(idat.rij) );
1008	–	}
1009	–	}
1010	–	if (i_is_Dipole) {
1011	–	if (j_is_Charge) {
1012	–	ri2 = riji * riji;
1013	–	ri3 = ri2 * riji;
1014	–	pref = (idat.electroMult) pre12_ * q_j * mu_i;
1015	–	vterm = - pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
1016	–	myPot += (idat.sw) vterm;
1017	–	dVdr += (idat.sw) pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
1018	–	duduz_i += (idat.sw) pref * rhat * (ri2 - preRF2_ * *(idat.rij));
1019	–	}
1020	–	}
1021	–
1022	–	// accumulate the forces and torques resulting from the self term
1023	–	idat.pot[ELECTROSTATIC_FAMILY] += myPot;
1024	–	*(idat.f1) += dVdr;
1025	–
987		if (i_is_Dipole)
988	<	*(idat.t1) -= cross(uz_i, duduz_i);
988	>	*(idat.t1) -= cross(uz_i, indirect_duduz_i);
989		if (j_is_Dipole)
990	<	*(idat.t2) -= cross(uz_j, duduz_j);
990	>	*(idat.t2) -= cross(uz_j, indirect_duduz_j);
991		}
992	<	}
992	>
993	>
994	>	return;
995	>	}
996
997		void Electrostatic::calcSelfCorrection(SelfData &sdat) {
998		RealType mu1, preVal, chg1, self;
999
1000		if (!initialized_) initialize();
1001	<
1001	>
1002		ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1003
1004		// logicals
1041	–
1005		bool i_is_Charge = data.is_Charge;
1006		bool i_is_Dipole = data.is_Dipole;
1007
#	Line 1046 \| Line 1009 \| namespace OpenMD {
1009		if (i_is_Dipole) {
1010		mu1 = data.dipole_moment;
1011		preVal = pre22_ * preRF2_ * mu1 * mu1;
1012	<	sdat.pot[2] -= 0.5 * preVal;
1012	>	((sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 preVal;
1013
1014		// The self-correction term adds into the reaction field vector
1015		Vector3d uz_i = sdat.eFrame->getColumn(2);
#	Line 1063 \| Line 1026 \| namespace OpenMD {
1026		} else {
1027		self = - 0.5 * rcuti_ * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1028		}
1029	<	sdat.pot[ELECTROSTATIC_FAMILY] += self;
1029	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1030		}
1031		}
1032		}

Diff Legend

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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents): Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs. Revision 1718 by gezelter, Thu May 24 01:29:59 2012 UTC

Diff Legend

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs.
Revision 1718 by gezelter, Thu May 24 01:29:59 2012 UTC