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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1535 by gezelter, Fri Dec 31 18:31:56 2010 UTC vs.
Revision 1616 by gezelter, Tue Aug 30 15:45:35 2011 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).
40		*/
#	Line 52 \| Line 52 \| namespace OpenMD {
52		namespace OpenMD {
53
54		Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55	<	forceField_(NULL) {}
55	>	forceField_(NULL), info_(NULL),
56	>	haveCutoffRadius_(false),
57	>	haveDampingAlpha_(false),
58	>	haveDielectric_(false),
59	>	haveElectroSpline_(false)
60	>	{}
61
62		void Electrostatic::initialize() {
63	+
64	+	Globals* simParams_ = info_->getSimParams();
65
59	–	Globals* simParams_;
60	–
66		summationMap_["HARD"] = esm_HARD;
67	+	summationMap_["NONE"] = esm_HARD;
68		summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
69		summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
70		summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
#	Line 97 \| Line 103 \| namespace OpenMD {
103		screeningMethod_ = UNDAMPED;
104		dielectric_ = 1.0;
105		one_third_ = 1.0 / 3.0;
100	–	haveCutoffRadius_ = false;
101	–	haveDampingAlpha_ = false;
102	–	haveDielectric_ = false;
103	–	haveElectroSpline_ = false;
106
107		// check the summation method:
108		if (simParams_->haveElectrostaticSummationMethod()) {
#	Line 113 \| Line 115 \| namespace OpenMD {
115		} else {
116		// throw error
117		sprintf( painCave.errMsg,
118	<	"SimInfo error: Unknown electrostaticSummationMethod.\n"
118	>	"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
119		"\t(Input file specified %s .)\n"
120	<	"\telectrostaticSummationMethod must be one of: \"none\",\n"
120	>	"\telectrostaticSummationMethod must be one of: \"hard\",\n"
121		"\t\"shifted_potential\", \"shifted_force\", or \n"
122		"\t\"reaction_field\".\n", myMethod.c_str() );
123		painCave.isFatal = 1;
#	Line 248 \| Line 250 \| namespace OpenMD {
250		preRF2_ = 2.0 * preRF_;
251		}
252
253	<	RealType dx = cutoffRadius_ / RealType(np_ - 1);
253	>	// Add a 2 angstrom safety window to deal with cutoffGroups that
254	>	// have charged atoms longer than the cutoffRadius away from each
255	>	// other. Splining may not be the best choice here. Direct calls
256	>	// to erfc might be preferrable.
257	>
258	>	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
259		RealType rval;
260		vector<RealType> rvals;
261		vector<RealType> yvals;
#	Line 407 \| Line 414 \| namespace OpenMD {
414		return;
415		}
416
417	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
418	<	RealType theRSW ) {
412	<	cutoffRadius_ = theECR;
417	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
418	>	cutoffRadius_ = rCut;
419		rrf_ = cutoffRadius_;
414	–	rt_ = theRSW;
420		haveCutoffRadius_ = true;
421		}
422	+
423	+	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
424	+	rt_ = rSwitch;
425	+	}
426		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
427		summationMethod_ = esm;
428		}
#	Line 429 \| Line 438 \| namespace OpenMD {
438		haveDielectric_ = true;
439		}
440
441	<	void Electrostatic::calcForce(InteractionData idat) {
441	>	void Electrostatic::calcForce(InteractionData &idat) {
442
443		// utility variables. Should clean these up and use the Vector3d and
444		// Mat3x3d to replace as many as we can in future versions:
#	Line 443 \| Line 452 \| namespace OpenMD {
452		RealType ct_i, ct_j, ct_ij, a1;
453		RealType riji, ri, ri2, ri3, ri4;
454		RealType pref, vterm, epot, dudr;
455	+	RealType vpair(0.0);
456		RealType scale, sc2;
457		RealType pot_term, preVal, rfVal;
458		RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
459		RealType preSw, preSwSc;
460		RealType c1, c2, c3, c4;
461	<	RealType erfcVal, derfcVal;
461	>	RealType erfcVal(1.0), derfcVal(0.0);
462		RealType BigR;
463
464		Vector3d Q_i, Q_j;
#	Line 459 \| Line 469 \| namespace OpenMD {
469		Vector3d rhatdot2, rhatc4;
470		Vector3d dVdr;
471
472	+	// variables for indirect (reaction field) interactions for excluded pairs:
473	+	RealType indirect_Pot(0.0);
474	+	RealType indirect_vpair(0.0);
475	+	Vector3d indirect_dVdr(V3Zero);
476	+	Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
477	+
478		pair<RealType, RealType> res;
479
480		if (!initialized_) initialize();
481
482	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
483	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
482	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
483	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
484
485		// some variables we'll need independent of electrostatic type:
486
487	<	riji = 1.0 / idat.rij;
488	<	Vector3d rhat = idat.d * riji;
487	>	riji = 1.0 / *(idat.rij) ;
488	>	Vector3d rhat = (idat.d) riji;
489
490		// logicals
491
#	Line 483 \| Line 499 \| namespace OpenMD {
499		bool j_is_SplitDipole = data2.is_SplitDipole;
500		bool j_is_Quadrupole = data2.is_Quadrupole;
501
502	<	if (i_is_Charge)
502	>	if (i_is_Charge) {
503		q_i = data1.charge;
504	+	if (idat.excluded) {
505	+	*(idat.skippedCharge2) += q_i;
506	+	}
507	+	}
508
509		if (i_is_Dipole) {
510		mu_i = data1.dipole_moment;
511	<	uz_i = idat.eFrame1.getColumn(2);
511	>	uz_i = idat.eFrame1->getColumn(2);
512
513		ct_i = dot(uz_i, rhat);
514
#	Line 504 \| Line 524 \| namespace OpenMD {
524		qyy_i = Q_i.y();
525		qzz_i = Q_i.z();
526
527	<	ux_i = idat.eFrame1.getColumn(0);
528	<	uy_i = idat.eFrame1.getColumn(1);
529	<	uz_i = idat.eFrame1.getColumn(2);
527	>	ux_i = idat.eFrame1->getColumn(0);
528	>	uy_i = idat.eFrame1->getColumn(1);
529	>	uz_i = idat.eFrame1->getColumn(2);
530
531		cx_i = dot(ux_i, rhat);
532		cy_i = dot(uy_i, rhat);
#	Line 517 \| Line 537 \| namespace OpenMD {
537		duduz_i = V3Zero;
538		}
539
540	<	if (j_is_Charge)
540	>	if (j_is_Charge) {
541		q_j = data2.charge;
542	+	if (idat.excluded) {
543	+	*(idat.skippedCharge1) += q_j;
544	+	}
545	+	}
546
547	+
548		if (j_is_Dipole) {
549		mu_j = data2.dipole_moment;
550	<	uz_j = idat.eFrame2.getColumn(2);
550	>	uz_j = idat.eFrame2->getColumn(2);
551
552		ct_j = dot(uz_j, rhat);
553
#	Line 538 \| Line 563 \| namespace OpenMD {
563		qyy_j = Q_j.y();
564		qzz_j = Q_j.z();
565
566	<	ux_j = idat.eFrame2.getColumn(0);
567	<	uy_j = idat.eFrame2.getColumn(1);
568	<	uz_j = idat.eFrame2.getColumn(2);
566	>	ux_j = idat.eFrame2->getColumn(0);
567	>	uy_j = idat.eFrame2->getColumn(1);
568	>	uz_j = idat.eFrame2->getColumn(2);
569
570		cx_j = dot(ux_j, rhat);
571		cy_j = dot(uy_j, rhat);
#	Line 559 \| Line 584 \| namespace OpenMD {
584		if (j_is_Charge) {
585		if (screeningMethod_ == DAMPED) {
586		// assemble the damping variables
587	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
588	<	erfcVal = res.first;
589	<	derfcVal = res.second;
587	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
588	>	//erfcVal = res.first;
589	>	//derfcVal = res.second;
590	>
591	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
592	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
593	>
594		c1 = erfcVal * riji;
595		c2 = (-derfcVal + c1) * riji;
596		} else {
#	Line 569 \| Line 598 \| namespace OpenMD {
598		c2 = c1 * riji;
599		}
600
601	<	preVal = idat.electroMult * pre11_ * q_i * q_j;
601	>	preVal = (idat.electroMult) pre11_ * q_i * q_j;
602
603		if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
604		vterm = preVal * (c1 - c1c_);
605	<	dudr = -idat.sw * preVal * c2;
605	>	dudr = - (idat.sw) preVal * c2;
606
607		} else if (summationMethod_ == esm_SHIFTED_FORCE) {
608	<	vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
609	<	dudr = idat.sw * preVal * (c2c_ - c2);
608	>	vterm = preVal * ( c1 - c1c_ + c2c_( (idat.rij) - cutoffRadius_) );
609	>	dudr = (idat.sw) preVal * (c2c_ - c2);
610
611		} else if (summationMethod_ == esm_REACTION_FIELD) {
612	<	rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
612	>	rfVal = preRF_ * (idat.rij) *(idat.rij);
613	>
614		vterm = preVal * ( riji + rfVal );
615	<	dudr = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
615	>	dudr = (idat.sw) preVal * ( 2.0 * rfVal - riji ) * riji;
616	>
617	>	// if this is an excluded pair, there are still indirect
618	>	// interactions via the reaction field we must worry about:
619
620	+	if (idat.excluded) {
621	+	indirect_vpair += preVal * rfVal;
622	+	indirect_Pot += (idat.sw) preVal * rfVal;
623	+	indirect_dVdr += (idat.sw) preVal * 2.0 * rfVal * riji * rhat;
624	+	}
625	+
626		} else {
588	–	vterm = preVal * riji * erfcVal;
627
628	<	dudr = - idat.sw * preVal * c2;
628	>	vterm = preVal * riji * erfcVal;
629	>	dudr = - (idat.sw) preVal * c2;
630
631		}
593	–
594	–	idat.vpair += vterm;
595	–	epot += idat.sw * vterm;
632
633	<	dVdr += dudr * rhat;
633	>	vpair += vterm;
634	>	epot += (idat.sw) vterm;
635	>	dVdr += dudr * rhat;
636		}
637
638		if (j_is_Dipole) {
639		// pref is used by all the possible methods
640	<	pref = idat.electroMult * pre12_ * q_i * mu_j;
641	<	preSw = idat.sw * pref;
640	>	pref = (idat.electroMult) pre12_ * q_i * mu_j;
641	>	preSw = (idat.sw) pref;
642
643		if (summationMethod_ == esm_REACTION_FIELD) {
644		ri2 = riji * riji;
645		ri3 = ri2 * riji;
646
647	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
648	<	idat.vpair += vterm;
649	<	epot += idat.sw * vterm;
647	>	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
648	>	vpair += vterm;
649	>	epot += (idat.sw) vterm;
650
651		dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
652	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);
652	>	duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
653
654	+	// Even if we excluded this pair from direct interactions,
655	+	// we still have the reaction-field-mediated charge-dipole
656	+	// interaction:
657	+
658	+	if (idat.excluded) {
659	+	indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
660	+	indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
661	+	indirect_dVdr += preSw * preRF2_ * uz_j;
662	+	indirect_duduz_j += preSw * rhat * preRF2_ * *(idat.rij);
663	+	}
664	+
665		} else {
666		// determine the inverse r used if we have split dipoles
667		if (j_is_SplitDipole) {
668	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
668	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
669		ri = 1.0 / BigR;
670	<	scale = idat.rij * ri;
670	>	scale = (idat.rij) ri;
671		} else {
672		ri = riji;
673		scale = 1.0;
#	Line 628 \| Line 677 \| namespace OpenMD {
677
678		if (screeningMethod_ == DAMPED) {
679		// assemble the damping variables
680	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
681	<	erfcVal = res.first;
682	<	derfcVal = res.second;
680	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
681	>	//erfcVal = res.first;
682	>	//derfcVal = res.second;
683	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
684	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
685		c1 = erfcVal * ri;
686		c2 = (-derfcVal + c1) * ri;
687		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 645 \| Line 696 \| namespace OpenMD {
696		// calculate the potential
697		pot_term = scale * c2;
698		vterm = -pref * ct_j * pot_term;
699	<	idat.vpair += vterm;
700	<	epot += idat.sw * vterm;
699	>	vpair += vterm;
700	>	epot += (idat.sw) vterm;
701
702		// calculate derivatives for forces and torques
703
#	Line 661 \| Line 712 \| namespace OpenMD {
712		cx2 = cx_j * cx_j;
713		cy2 = cy_j * cy_j;
714		cz2 = cz_j * cz_j;
715	<	pref = idat.electroMult * pre14_ * q_i * one_third_;
715	>	pref = (idat.electroMult) pre14_ * q_i * one_third_;
716
717		if (screeningMethod_ == DAMPED) {
718		// assemble the damping variables
719	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
720	<	erfcVal = res.first;
721	<	derfcVal = res.second;
719	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
720	>	//erfcVal = res.first;
721	>	//derfcVal = res.second;
722	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
723	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
724		c1 = erfcVal * riji;
725		c2 = (-derfcVal + c1) * riji;
726		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 680 \| Line 733 \| namespace OpenMD {
733		}
734
735		// precompute variables for convenience
736	<	preSw = idat.sw * pref;
736	>	preSw = (idat.sw) pref;
737		c2ri = c2 * riji;
738		c3ri = c3 * riji;
739	<	c4rij = c4 * idat.rij;
739	>	c4rij = c4 * *(idat.rij) ;
740		rhatdot2 = 2.0 * rhat * c3;
741		rhatc4 = rhat * c4rij;
742
#	Line 692 \| Line 745 \| namespace OpenMD {
745		qyy_j * (cy2*c3 - c2ri) +
746		qzz_j * (cz2*c3 - c2ri) );
747		vterm = pref * pot_term;
748	<	idat.vpair += vterm;
749	<	epot += idat.sw * vterm;
748	>	vpair += vterm;
749	>	epot += (idat.sw) vterm;
750
751		// calculate derivatives for the forces and torques
752
#	Line 711 \| Line 764 \| namespace OpenMD {
764
765		if (j_is_Charge) {
766		// variables used by all the methods
767	<	pref = idat.electroMult * pre12_ * q_j * mu_i;
768	<	preSw = idat.sw * pref;
767	>	pref = (idat.electroMult) pre12_ * q_j * mu_i;
768	>	preSw = (idat.sw) pref;
769
770		if (summationMethod_ == esm_REACTION_FIELD) {
771
772		ri2 = riji * riji;
773		ri3 = ri2 * riji;
774
775	<	vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
776	<	idat.vpair += vterm;
777	<	epot += idat.sw * vterm;
775	>	vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
776	>	vpair += vterm;
777	>	epot += (idat.sw) vterm;
778
779		dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
780
781	<	duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
781	>	duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
782	>
783	>	// Even if we excluded this pair from direct interactions,
784	>	// we still have the reaction-field-mediated charge-dipole
785	>	// interaction:
786	>
787	>	if (idat.excluded) {
788	>	indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
789	>	indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
790	>	indirect_dVdr += -preSw * preRF2_ * uz_i;
791	>	indirect_duduz_i += -preSw * rhat * preRF2_ * *(idat.rij);
792	>	}
793
794		} else {
795
796		// determine inverse r if we are using split dipoles
797		if (i_is_SplitDipole) {
798	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
798	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
799		ri = 1.0 / BigR;
800	<	scale = idat.rij * ri;
800	>	scale = (idat.rij) ri;
801		} else {
802		ri = riji;
803		scale = 1.0;
#	Line 743 \| Line 807 \| namespace OpenMD {
807
808		if (screeningMethod_ == DAMPED) {
809		// assemble the damping variables
810	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
811	<	erfcVal = res.first;
812	<	derfcVal = res.second;
810	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
811	>	//erfcVal = res.first;
812	>	//derfcVal = res.second;
813	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
814	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
815		c1 = erfcVal * ri;
816		c2 = (-derfcVal + c1) * ri;
817		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 760 \| Line 826 \| namespace OpenMD {
826		// calculate the potential
827		pot_term = c2 * scale;
828		vterm = pref * ct_i * pot_term;
829	<	idat.vpair += vterm;
830	<	epot += idat.sw * vterm;
829	>	vpair += vterm;
830	>	epot += (idat.sw) vterm;
831
832		// calculate derivatives for the forces and torques
833		dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
#	Line 773 \| Line 839 \| namespace OpenMD {
839		// variables used by all methods
840		ct_ij = dot(uz_i, uz_j);
841
842	<	pref = idat.electroMult * pre22_ * mu_i * mu_j;
843	<	preSw = idat.sw * pref;
842	>	pref = (idat.electroMult) pre22_ * mu_i * mu_j;
843	>	preSw = (idat.sw) pref;
844
845		if (summationMethod_ == esm_REACTION_FIELD) {
846		ri2 = riji * riji;
#	Line 783 \| Line 849 \| namespace OpenMD {
849
850		vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
851		preRF2_ * ct_ij );
852	<	idat.vpair += vterm;
853	<	epot += idat.sw * vterm;
852	>	vpair += vterm;
853	>	epot += (idat.sw) vterm;
854
855		a1 = 5.0 * ct_i * ct_j - ct_ij;
856
#	Line 793 \| Line 859 \| namespace OpenMD {
859		duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
860		duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
861
862	+	if (idat.excluded) {
863	+	indirect_vpair += - pref * preRF2_ * ct_ij;
864	+	indirect_Pot += - preSw * preRF2_ * ct_ij;
865	+	indirect_duduz_i += -preSw * preRF2_ * uz_j;
866	+	indirect_duduz_j += -preSw * preRF2_ * uz_i;
867	+	}
868	+
869		} else {
870
871		if (i_is_SplitDipole) {
872		if (j_is_SplitDipole) {
873	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
873	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i + 0.25 * d_j * d_j);
874		} else {
875	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
875	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
876		}
877		ri = 1.0 / BigR;
878	<	scale = idat.rij * ri;
878	>	scale = (idat.rij) ri;
879		} else {
880		if (j_is_SplitDipole) {
881	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
881	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
882		ri = 1.0 / BigR;
883	<	scale = idat.rij * ri;
883	>	scale = (idat.rij) ri;
884		} else {
885		ri = riji;
886		scale = 1.0;
#	Line 815 \| Line 888 \| namespace OpenMD {
888		}
889		if (screeningMethod_ == DAMPED) {
890		// assemble damping variables
891	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
892	<	erfcVal = res.first;
893	<	derfcVal = res.second;
891	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
892	>	//erfcVal = res.first;
893	>	//derfcVal = res.second;
894	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
895	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
896		c1 = erfcVal * ri;
897		c2 = (-derfcVal + c1) * ri;
898		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 837 \| Line 912 \| namespace OpenMD {
912		preSwSc = preSw * scale;
913		c2ri = c2 * ri;
914		c3ri = c3 * ri;
915	<	c4rij = c4 * idat.rij;
915	>	c4rij = c4 * *(idat.rij) ;
916
917		// calculate the potential
918		pot_term = (ct_ij * c2ri - ctidotj * c3);
919		vterm = pref * pot_term;
920	<	idat.vpair += vterm;
921	<	epot += idat.sw * vterm;
920	>	vpair += vterm;
921	>	epot += (idat.sw) vterm;
922
923		// calculate derivatives for the forces and torques
924		dVdr += preSwSc * ( ctidotj * rhat * c4rij -
#	Line 862 \| Line 937 \| namespace OpenMD {
937		cy2 = cy_i * cy_i;
938		cz2 = cz_i * cz_i;
939
940	<	pref = idat.electroMult * pre14_ * q_j * one_third_;
940	>	pref = (idat.electroMult) pre14_ * q_j * one_third_;
941
942		if (screeningMethod_ == DAMPED) {
943		// assemble the damping variables
944	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
945	<	erfcVal = res.first;
946	<	derfcVal = res.second;
944	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
945	>	//erfcVal = res.first;
946	>	//derfcVal = res.second;
947	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
948	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
949		c1 = erfcVal * riji;
950		c2 = (-derfcVal + c1) * riji;
951		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 881 \| Line 958 \| namespace OpenMD {
958		}
959
960		// precompute some variables for convenience
961	<	preSw = idat.sw * pref;
961	>	preSw = (idat.sw) pref;
962		c2ri = c2 * riji;
963		c3ri = c3 * riji;
964	<	c4rij = c4 * idat.rij;
964	>	c4rij = c4 * *(idat.rij) ;
965		rhatdot2 = 2.0 * rhat * c3;
966		rhatc4 = rhat * c4rij;
967
#	Line 894 \| Line 971 \| namespace OpenMD {
971		qzz_i * (cz2 * c3 - c2ri) );
972
973		vterm = pref * pot_term;
974	<	idat.vpair += vterm;
975	<	epot += idat.sw * vterm;
974	>	vpair += vterm;
975	>	epot += (idat.sw) vterm;
976
977		// calculate the derivatives for the forces and torques
978
#	Line 909 \| Line 986 \| namespace OpenMD {
986		}
987		}
988
912	–	idat.pot += epot;
913	–	idat.f1 += dVdr;
989
990	<	if (i_is_Dipole \|\| i_is_Quadrupole)
991	<	idat.t1 -= cross(uz_i, duduz_i);
992	<	if (i_is_Quadrupole) {
993	<	idat.t1 -= cross(ux_i, dudux_i);
994	<	idat.t1 -= cross(uy_i, duduy_i);
995	<	}
990	>	if (!idat.excluded) {
991	>	*(idat.vpair) += vpair;
992	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
993	>	*(idat.f1) += dVdr;
994	>
995	>	if (i_is_Dipole \|\| i_is_Quadrupole)
996	>	*(idat.t1) -= cross(uz_i, duduz_i);
997	>	if (i_is_Quadrupole) {
998	>	*(idat.t1) -= cross(ux_i, dudux_i);
999	>	*(idat.t1) -= cross(uy_i, duduy_i);
1000	>	}
1001	>
1002	>	if (j_is_Dipole \|\| j_is_Quadrupole)
1003	>	*(idat.t2) -= cross(uz_j, duduz_j);
1004	>	if (j_is_Quadrupole) {
1005	>	*(idat.t2) -= cross(uz_j, dudux_j);
1006	>	*(idat.t2) -= cross(uz_j, duduy_j);
1007	>	}
1008
1009	<	if (j_is_Dipole \|\| j_is_Quadrupole)
923	<	idat.t2 -= cross(uz_j, duduz_j);
924	<	if (j_is_Quadrupole) {
925	<	idat.t2 -= cross(uz_j, dudux_j);
926	<	idat.t2 -= cross(uz_j, duduy_j);
927	<	}
1009	>	} else {
1010
1011	<	return;
1012	<	}
1011	>	// only accumulate the forces and torques resulting from the
1012	>	// indirect reaction field terms.
1013
1014	<	void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
1014	>	*(idat.vpair) += indirect_vpair;
1015	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1016	>	*(idat.f1) += indirect_dVdr;
1017	>
1018	>	if (i_is_Dipole)
1019	>	*(idat.t1) -= cross(uz_i, indirect_duduz_i);
1020	>	if (j_is_Dipole)
1021	>	*(idat.t2) -= cross(uz_j, indirect_duduz_j);
1022	>	}
1023
934	–	if (!initialized_) initialize();
935	–
936	–	ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
937	–	ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
938	–
939	–	// logicals
1024
1025	<	bool i_is_Charge = data1.is_Charge;
1026	<	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;
1025	>	return;
1026	>	}
1027
1028	<	// 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	<	skdat.skippedCharge2 += q_i;
955	<	}
956	<
957	<	if (j_is_Charge) {
958	<	q_j = data2.charge;
959	<	skdat.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 q_i, mu_i, ct_i;
967	<	RealType q_j, mu_j, ct_j;
968	<	RealType preVal, rfVal, vterm, dudr, pref, myPot;
969	<	Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
970	<
971	<	// some variables we'll need independent of electrostatic type:
972	<
973	<	riji = 1.0 / skdat.rij;
974	<	rhat = skdat.d * riji;
975	<
976	<	if (i_is_Dipole) {
977	<	mu_i = data1.dipole_moment;
978	<	uz_i = skdat.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 = skdat.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 = skdat.electroMult * pre11_ * q_i * q_j;
993	<	rfVal = preRF_ * skdat.rij * skdat.rij;
994	<	vterm = preVal * rfVal;
995	<	myPot += skdat.sw * vterm;
996	<	dudr = skdat.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 = skdat.electroMult * pre12_ * q_i * mu_j;
1004	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1005	<	myPot += skdat.sw * vterm;
1006	<	dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1007	<	duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1008	<	}
1009	<	}
1010	<	if (i_is_Dipole) {
1011	<	if (j_is_Charge) {
1012	<	ri2 = riji * riji;
1013	<	ri3 = ri2 * riji;
1014	<	pref = skdat.electroMult * pre12_ * q_j * mu_i;
1015	<	vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1016	<	myPot += skdat.sw * vterm;
1017	<	dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
1018	<	duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1019	<	}
1020	<	}
1021	<
1022	<	// accumulate the forces and torques resulting from the self term
1023	<	skdat.pot += myPot;
1024	<	skdat.f1 += dVdr;
1025	<
1026	<	if (i_is_Dipole)
1027	<	skdat.t1 -= cross(uz_i, duduz_i);
1028	<	if (j_is_Dipole)
1029	<	skdat.t2 -= cross(uz_j, duduz_j);
1030	<	}
1031	<	}
1032	<
1033	<	void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1028	>	void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1029		RealType mu1, preVal, chg1, self;
1030
1031		if (!initialized_) initialize();
1032	<
1033	<	ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1032	>
1033	>	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1034
1035		// logicals
1041	–
1036		bool i_is_Charge = data.is_Charge;
1037		bool i_is_Dipole = data.is_Dipole;
1038
#	Line 1046 \| Line 1040 \| namespace OpenMD {
1040		if (i_is_Dipole) {
1041		mu1 = data.dipole_moment;
1042		preVal = pre22_ * preRF2_ * mu1 * mu1;
1043	<	scdat.pot -= 0.5 * preVal;
1043	>	((sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 preVal;
1044
1045		// The self-correction term adds into the reaction field vector
1046	<	Vector3d uz_i = scdat.eFrame.getColumn(2);
1046	>	Vector3d uz_i = sdat.eFrame->getColumn(2);
1047		Vector3d ei = preVal * uz_i;
1048
1049		// This looks very wrong. A vector crossed with itself is zero.
1050	<	scdat.t -= cross(uz_i, ei);
1050	>	*(sdat.t) -= cross(uz_i, ei);
1051		}
1052		} else if (summationMethod_ == esm_SHIFTED_FORCE \|\| summationMethod_ == esm_SHIFTED_POTENTIAL) {
1053		if (i_is_Charge) {
1054		chg1 = data.charge;
1055		if (screeningMethod_ == DAMPED) {
1056	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1056	>	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1057		} else {
1058	<	self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1058	>	self = - 0.5 * rcuti_ * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1059		}
1060	<	scdat.pot += self;
1060	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1061		}
1062		}
1063		}
1064
1065	<	RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1065	>	RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
1066		// This seems to work moderately well as a default. There's no
1067		// inherent scale for 1/r interactions that we can standardize.
1068		// 12 angstroms seems to be a reasonably good guess for most

Diff Legend

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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents): Revision 1535 by gezelter, Fri Dec 31 18:31:56 2010 UTC vs. Revision 1616 by gezelter, Tue Aug 30 15:45:35 2011 UTC

Diff Legend

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1535 by gezelter, Fri Dec 31 18:31:56 2010 UTC vs.
Revision 1616 by gezelter, Tue Aug 30 15:45:35 2011 UTC