[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 1613 by gezelter, Thu Aug 18 20:18:19 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_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
68		summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
#	Line 97 \| Line 102 \| namespace OpenMD {
102		screeningMethod_ = UNDAMPED;
103		dielectric_ = 1.0;
104		one_third_ = 1.0 / 3.0;
100	–	haveCutoffRadius_ = false;
101	–	haveDampingAlpha_ = false;
102	–	haveDielectric_ = false;
103	–	haveElectroSpline_ = false;
105
106		// check the summation method:
107		if (simParams_->haveElectrostaticSummationMethod()) {
#	Line 113 \| Line 114 \| namespace OpenMD {
114		} else {
115		// throw error
116		sprintf( painCave.errMsg,
117	<	"SimInfo error: Unknown electrostaticSummationMethod.\n"
117	>	"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
118		"\t(Input file specified %s .)\n"
119		"\telectrostaticSummationMethod must be one of: \"none\",\n"
120		"\t\"shifted_potential\", \"shifted_force\", or \n"
#	Line 248 \| Line 249 \| namespace OpenMD {
249		preRF2_ = 2.0 * preRF_;
250		}
251
252	<	RealType dx = cutoffRadius_ / RealType(np_ - 1);
252	>	// Add a 2 angstrom safety window to deal with cutoffGroups that
253	>	// have charged atoms longer than the cutoffRadius away from each
254	>	// other. Splining may not be the best choice here. Direct calls
255	>	// to erfc might be preferrable.
256	>
257	>	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
258		RealType rval;
259		vector<RealType> rvals;
260		vector<RealType> yvals;
#	Line 407 \| Line 413 \| namespace OpenMD {
413		return;
414		}
415
416	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
417	<	RealType theRSW ) {
412	<	cutoffRadius_ = theECR;
416	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
417	>	cutoffRadius_ = rCut;
418		rrf_ = cutoffRadius_;
414	–	rt_ = theRSW;
419		haveCutoffRadius_ = true;
420	+	}
421	+
422	+	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
423	+	rt_ = rSwitch;
424		}
425		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
426		summationMethod_ = esm;
#	Line 429 \| Line 437 \| namespace OpenMD {
437		haveDielectric_ = true;
438		}
439
440	<	void Electrostatic::calcForce(InteractionData idat) {
440	>	void Electrostatic::calcForce(InteractionData &idat) {
441
442		// utility variables. Should clean these up and use the Vector3d and
443		// Mat3x3d to replace as many as we can in future versions:
#	Line 443 \| Line 451 \| namespace OpenMD {
451		RealType ct_i, ct_j, ct_ij, a1;
452		RealType riji, ri, ri2, ri3, ri4;
453		RealType pref, vterm, epot, dudr;
454	+	RealType vpair(0.0);
455		RealType scale, sc2;
456		RealType pot_term, preVal, rfVal;
457		RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
458		RealType preSw, preSwSc;
459		RealType c1, c2, c3, c4;
460	<	RealType erfcVal, derfcVal;
460	>	RealType erfcVal(1.0), derfcVal(0.0);
461		RealType BigR;
462
463		Vector3d Q_i, Q_j;
#	Line 459 \| Line 468 \| namespace OpenMD {
468		Vector3d rhatdot2, rhatc4;
469		Vector3d dVdr;
470
471	+	// variables for indirect (reaction field) interactions for excluded pairs:
472	+	RealType indirect_Pot(0.0);
473	+	RealType indirect_vpair(0.0);
474	+	Vector3d indirect_dVdr(V3Zero);
475	+	Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
476	+
477		pair<RealType, RealType> res;
478
479		if (!initialized_) initialize();
480
481	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
482	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
481	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
482	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
483
484		// some variables we'll need independent of electrostatic type:
485
486	<	riji = 1.0 / idat.rij;
487	<	Vector3d rhat = idat.d * riji;
486	>	riji = 1.0 / *(idat.rij) ;
487	>	Vector3d rhat = (idat.d) riji;
488
489		// logicals
490
#	Line 483 \| Line 498 \| namespace OpenMD {
498		bool j_is_SplitDipole = data2.is_SplitDipole;
499		bool j_is_Quadrupole = data2.is_Quadrupole;
500
501	<	if (i_is_Charge)
501	>	if (i_is_Charge) {
502		q_i = data1.charge;
503	+	if (idat.excluded) {
504	+	*(idat.skippedCharge2) += q_i;
505	+	}
506	+	}
507
508		if (i_is_Dipole) {
509		mu_i = data1.dipole_moment;
510	<	uz_i = idat.eFrame1.getColumn(2);
510	>	uz_i = idat.eFrame1->getColumn(2);
511
512		ct_i = dot(uz_i, rhat);
513
#	Line 504 \| Line 523 \| namespace OpenMD {
523		qyy_i = Q_i.y();
524		qzz_i = Q_i.z();
525
526	<	ux_i = idat.eFrame1.getColumn(0);
527	<	uy_i = idat.eFrame1.getColumn(1);
528	<	uz_i = idat.eFrame1.getColumn(2);
526	>	ux_i = idat.eFrame1->getColumn(0);
527	>	uy_i = idat.eFrame1->getColumn(1);
528	>	uz_i = idat.eFrame1->getColumn(2);
529
530		cx_i = dot(ux_i, rhat);
531		cy_i = dot(uy_i, rhat);
#	Line 517 \| Line 536 \| namespace OpenMD {
536		duduz_i = V3Zero;
537		}
538
539	<	if (j_is_Charge)
539	>	if (j_is_Charge) {
540		q_j = data2.charge;
541	+	if (idat.excluded) {
542	+	*(idat.skippedCharge1) += q_j;
543	+	}
544	+	}
545
546	+
547		if (j_is_Dipole) {
548		mu_j = data2.dipole_moment;
549	<	uz_j = idat.eFrame2.getColumn(2);
549	>	uz_j = idat.eFrame2->getColumn(2);
550
551		ct_j = dot(uz_j, rhat);
552
#	Line 538 \| Line 562 \| namespace OpenMD {
562		qyy_j = Q_j.y();
563		qzz_j = Q_j.z();
564
565	<	ux_j = idat.eFrame2.getColumn(0);
566	<	uy_j = idat.eFrame2.getColumn(1);
567	<	uz_j = idat.eFrame2.getColumn(2);
565	>	ux_j = idat.eFrame2->getColumn(0);
566	>	uy_j = idat.eFrame2->getColumn(1);
567	>	uz_j = idat.eFrame2->getColumn(2);
568
569		cx_j = dot(ux_j, rhat);
570		cy_j = dot(uy_j, rhat);
#	Line 559 \| Line 583 \| namespace OpenMD {
583		if (j_is_Charge) {
584		if (screeningMethod_ == DAMPED) {
585		// assemble the damping variables
586	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
586	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
587		erfcVal = res.first;
588		derfcVal = res.second;
589		c1 = erfcVal * riji;
#	Line 569 \| Line 593 \| namespace OpenMD {
593		c2 = c1 * riji;
594		}
595
596	<	preVal = idat.electroMult * pre11_ * q_i * q_j;
596	>	preVal = (idat.electroMult) pre11_ * q_i * q_j;
597
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_ == esm_SHIFTED_FORCE) {
603	<	vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
604	<	dudr = idat.sw * preVal * (c2c_ - c2);
603	>	vterm = preVal * ( c1 - c1c_ + c2c_( (idat.rij) - cutoffRadius_) );
604	>	dudr = (idat.sw) preVal * (c2c_ - c2);
605
606		} else if (summationMethod_ == esm_REACTION_FIELD) {
607	<	rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
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 * 2.0 * rfVal * riji * rhat;
619	+	}
620	+
621		} else {
588	–	vterm = preVal * riji * erfcVal;
622
623	<	dudr = - idat.sw * preVal * c2;
623	>	vterm = preVal * riji * erfcVal;
624	>	dudr = - (idat.sw) preVal * c2;
625
626		}
593	–
594	–	idat.vpair += vterm;
595	–	epot += idat.sw * vterm;
627
628	<	dVdr += dudr * rhat;
628	>	vpair += vterm;
629	>	epot += (idat.sw) vterm;
630	>	dVdr += dudr * rhat;
631		}
632
633		if (j_is_Dipole) {
634		// pref is used by all the possible methods
635	<	pref = idat.electroMult * pre12_ * q_i * mu_j;
636	<	preSw = idat.sw * pref;
635	>	pref = (idat.electroMult) pre12_ * q_i * mu_j;
636	>	preSw = (idat.sw) pref;
637
638		if (summationMethod_ == esm_REACTION_FIELD) {
639		ri2 = riji * riji;
640		ri3 = ri2 * riji;
641
642	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
643	<	idat.vpair += vterm;
644	<	epot += idat.sw * vterm;
642	>	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
643	>	vpair += vterm;
644	>	epot += (idat.sw) vterm;
645
646		dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
647	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);
647	>	duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
648
649	+	// Even if we excluded this pair from direct interactions,
650	+	// we still have the reaction-field-mediated charge-dipole
651	+	// interaction:
652	+
653	+	if (idat.excluded) {
654	+	indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
655	+	indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
656	+	indirect_dVdr += preSw * preRF2_ * uz_j;
657	+	indirect_duduz_j += preSw * rhat * preRF2_ * *(idat.rij);
658	+	}
659	+
660		} else {
661		// determine the inverse r used if we have split dipoles
662		if (j_is_SplitDipole) {
663	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
663	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
664		ri = 1.0 / BigR;
665	<	scale = idat.rij * ri;
665	>	scale = (idat.rij) ri;
666		} else {
667		ri = riji;
668		scale = 1.0;
#	Line 628 \| Line 672 \| namespace OpenMD {
672
673		if (screeningMethod_ == DAMPED) {
674		// assemble the damping variables
675	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
675	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
676		erfcVal = res.first;
677		derfcVal = res.second;
678		c1 = erfcVal * ri;
#	Line 645 \| Line 689 \| namespace OpenMD {
689		// calculate the potential
690		pot_term = scale * c2;
691		vterm = -pref * ct_j * pot_term;
692	<	idat.vpair += vterm;
693	<	epot += idat.sw * vterm;
692	>	vpair += vterm;
693	>	epot += (idat.sw) vterm;
694
695		// calculate derivatives for forces and torques
696
#	Line 661 \| Line 705 \| namespace OpenMD {
705		cx2 = cx_j * cx_j;
706		cy2 = cy_j * cy_j;
707		cz2 = cz_j * cz_j;
708	<	pref = idat.electroMult * pre14_ * q_i * one_third_;
708	>	pref = (idat.electroMult) pre14_ * q_i * one_third_;
709
710		if (screeningMethod_ == DAMPED) {
711		// assemble the damping variables
712	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
712	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
713		erfcVal = res.first;
714		derfcVal = res.second;
715		c1 = erfcVal * riji;
#	Line 680 \| Line 724 \| namespace OpenMD {
724		}
725
726		// precompute variables for convenience
727	<	preSw = idat.sw * pref;
727	>	preSw = (idat.sw) pref;
728		c2ri = c2 * riji;
729		c3ri = c3 * riji;
730	<	c4rij = c4 * idat.rij;
730	>	c4rij = c4 * *(idat.rij) ;
731		rhatdot2 = 2.0 * rhat * c3;
732		rhatc4 = rhat * c4rij;
733
#	Line 692 \| Line 736 \| namespace OpenMD {
736		qyy_j * (cy2*c3 - c2ri) +
737		qzz_j * (cz2*c3 - c2ri) );
738		vterm = pref * pot_term;
739	<	idat.vpair += vterm;
740	<	epot += idat.sw * vterm;
739	>	vpair += vterm;
740	>	epot += (idat.sw) vterm;
741
742		// calculate derivatives for the forces and torques
743
#	Line 711 \| Line 755 \| namespace OpenMD {
755
756		if (j_is_Charge) {
757		// variables used by all the methods
758	<	pref = idat.electroMult * pre12_ * q_j * mu_i;
759	<	preSw = idat.sw * pref;
758	>	pref = (idat.electroMult) pre12_ * q_j * mu_i;
759	>	preSw = (idat.sw) pref;
760
761		if (summationMethod_ == esm_REACTION_FIELD) {
762
763		ri2 = riji * riji;
764		ri3 = ri2 * riji;
765
766	<	vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
767	<	idat.vpair += vterm;
768	<	epot += idat.sw * vterm;
766	>	vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
767	>	vpair += vterm;
768	>	epot += (idat.sw) vterm;
769
770		dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
771
772	<	duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
772	>	duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
773	>
774	>	// Even if we excluded this pair from direct interactions,
775	>	// we still have the reaction-field-mediated charge-dipole
776	>	// interaction:
777	>
778	>	if (idat.excluded) {
779	>	indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
780	>	indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
781	>	indirect_dVdr += -preSw * preRF2_ * uz_i;
782	>	indirect_duduz_i += -preSw * rhat * preRF2_ * *(idat.rij);
783	>	}
784
785		} else {
786
787		// determine inverse r if we are using split dipoles
788		if (i_is_SplitDipole) {
789	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
789	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
790		ri = 1.0 / BigR;
791	<	scale = idat.rij * ri;
791	>	scale = (idat.rij) ri;
792		} else {
793		ri = riji;
794		scale = 1.0;
#	Line 743 \| Line 798 \| namespace OpenMD {
798
799		if (screeningMethod_ == DAMPED) {
800		// assemble the damping variables
801	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
801	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
802		erfcVal = res.first;
803		derfcVal = res.second;
804		c1 = erfcVal * ri;
#	Line 760 \| Line 815 \| namespace OpenMD {
815		// calculate the potential
816		pot_term = c2 * scale;
817		vterm = pref * ct_i * pot_term;
818	<	idat.vpair += vterm;
819	<	epot += idat.sw * vterm;
818	>	vpair += vterm;
819	>	epot += (idat.sw) vterm;
820
821		// calculate derivatives for the forces and torques
822		dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
#	Line 773 \| Line 828 \| namespace OpenMD {
828		// variables used by all methods
829		ct_ij = dot(uz_i, uz_j);
830
831	<	pref = idat.electroMult * pre22_ * mu_i * mu_j;
832	<	preSw = idat.sw * pref;
831	>	pref = (idat.electroMult) pre22_ * mu_i * mu_j;
832	>	preSw = (idat.sw) pref;
833
834		if (summationMethod_ == esm_REACTION_FIELD) {
835		ri2 = riji * riji;
#	Line 783 \| Line 838 \| namespace OpenMD {
838
839		vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
840		preRF2_ * ct_ij );
841	<	idat.vpair += vterm;
842	<	epot += idat.sw * vterm;
841	>	vpair += vterm;
842	>	epot += (idat.sw) vterm;
843
844		a1 = 5.0 * ct_i * ct_j - ct_ij;
845
#	Line 792 \| Line 847 \| namespace OpenMD {
847
848		duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
849		duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
850	+
851	+	if (idat.excluded) {
852	+	indirect_vpair += - pref * preRF2_ * ct_ij;
853	+	indirect_Pot += - preSw * preRF2_ * ct_ij;
854	+	indirect_duduz_i += -preSw * preRF2_ * uz_j;
855	+	indirect_duduz_j += -preSw * preRF2_ * uz_i;
856	+	}
857
858		} else {
859
860		if (i_is_SplitDipole) {
861		if (j_is_SplitDipole) {
862	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
862	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i + 0.25 * d_j * d_j);
863		} else {
864	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
864	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
865		}
866		ri = 1.0 / BigR;
867	<	scale = idat.rij * ri;
867	>	scale = (idat.rij) ri;
868		} else {
869		if (j_is_SplitDipole) {
870	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
870	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
871		ri = 1.0 / BigR;
872	<	scale = idat.rij * ri;
872	>	scale = (idat.rij) ri;
873		} else {
874		ri = riji;
875		scale = 1.0;
#	Line 815 \| Line 877 \| namespace OpenMD {
877		}
878		if (screeningMethod_ == DAMPED) {
879		// assemble damping variables
880	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
880	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
881		erfcVal = res.first;
882		derfcVal = res.second;
883		c1 = erfcVal * ri;
#	Line 837 \| Line 899 \| namespace OpenMD {
899		preSwSc = preSw * scale;
900		c2ri = c2 * ri;
901		c3ri = c3 * ri;
902	<	c4rij = c4 * idat.rij;
902	>	c4rij = c4 * *(idat.rij) ;
903
904		// calculate the potential
905		pot_term = (ct_ij * c2ri - ctidotj * c3);
906		vterm = pref * pot_term;
907	<	idat.vpair += vterm;
908	<	epot += idat.sw * vterm;
907	>	vpair += vterm;
908	>	epot += (idat.sw) vterm;
909
910		// calculate derivatives for the forces and torques
911		dVdr += preSwSc * ( ctidotj * rhat * c4rij -
#	Line 862 \| Line 924 \| namespace OpenMD {
924		cy2 = cy_i * cy_i;
925		cz2 = cz_i * cz_i;
926
927	<	pref = idat.electroMult * pre14_ * q_j * one_third_;
927	>	pref = (idat.electroMult) pre14_ * q_j * one_third_;
928
929		if (screeningMethod_ == DAMPED) {
930		// assemble the damping variables
931	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
931	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
932		erfcVal = res.first;
933		derfcVal = res.second;
934		c1 = erfcVal * riji;
#	Line 881 \| Line 943 \| namespace OpenMD {
943		}
944
945		// precompute some variables for convenience
946	<	preSw = idat.sw * pref;
946	>	preSw = (idat.sw) pref;
947		c2ri = c2 * riji;
948		c3ri = c3 * riji;
949	<	c4rij = c4 * idat.rij;
949	>	c4rij = c4 * *(idat.rij) ;
950		rhatdot2 = 2.0 * rhat * c3;
951		rhatc4 = rhat * c4rij;
952
#	Line 894 \| Line 956 \| namespace OpenMD {
956		qzz_i * (cz2 * c3 - c2ri) );
957
958		vterm = pref * pot_term;
959	<	idat.vpair += vterm;
960	<	epot += idat.sw * vterm;
959	>	vpair += vterm;
960	>	epot += (idat.sw) vterm;
961
962		// calculate the derivatives for the forces and torques
963
#	Line 909 \| Line 971 \| namespace OpenMD {
971		}
972		}
973
912	–	idat.pot += epot;
913	–	idat.f1 += dVdr;
974
975	<	if (i_is_Dipole \|\| i_is_Quadrupole)
976	<	idat.t1 -= cross(uz_i, duduz_i);
977	<	if (i_is_Quadrupole) {
978	<	idat.t1 -= cross(ux_i, dudux_i);
919	<	idat.t1 -= cross(uy_i, duduy_i);
920	<	}
921	<
922	<	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	<	}
928	<
929	<	return;
930	<	}
931	<
932	<	void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
933	<
934	<	if (!initialized_) initialize();
935	<
936	<	ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
937	<	ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
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	<	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:
975	>	if (!idat.excluded) {
976	>	*(idat.vpair) += vpair;
977	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
978	>	*(idat.f1) += dVdr;
979
980	<	riji = 1.0 / skdat.rij;
981	<	rhat = skdat.d * riji;
982	<
983	<	if (i_is_Dipole) {
984	<	mu_i = data1.dipole_moment;
978	<	uz_i = skdat.eFrame1.getColumn(2);
979	<	ct_i = dot(uz_i, rhat);
980	<	duduz_i = V3Zero;
980	>	if (i_is_Dipole \|\| i_is_Quadrupole)
981	>	*(idat.t1) -= cross(uz_i, duduz_i);
982	>	if (i_is_Quadrupole) {
983	>	*(idat.t1) -= cross(ux_i, dudux_i);
984	>	*(idat.t1) -= cross(uy_i, duduy_i);
985		}
986	<
987	<	if (j_is_Dipole) {
988	<	mu_j = data2.dipole_moment;
989	<	uz_j = skdat.eFrame2.getColumn(2);
990	<	ct_j = dot(uz_j, rhat);
991	<	duduz_j = V3Zero;
986	>
987	>	if (j_is_Dipole \|\| j_is_Quadrupole)
988	>	*(idat.t2) -= cross(uz_j, duduz_j);
989	>	if (j_is_Quadrupole) {
990	>	*(idat.t2) -= cross(uz_j, dudux_j);
991	>	*(idat.t2) -= cross(uz_j, duduy_j);
992		}
993	<
994	<	if (i_is_Charge) {
995	<	if (j_is_Charge) {
996	<	preVal = skdat.electroMult * pre11_ * q_i * q_j;
997	<	rfVal = preRF_ * skdat.rij * skdat.rij;
998	<	vterm = preVal * rfVal;
999	<	myPot += skdat.sw * vterm;
1000	<	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	<	}
993	>
994	>	} else {
995	>
996	>	// only accumulate the forces and torques resulting from the
997	>	// indirect reaction field terms.
998	>	*(idat.vpair) += indirect_vpair;
999	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1000	>	*(idat.f1) += indirect_dVdr;
1001
1022	–	// accumulate the forces and torques resulting from the self term
1023	–	skdat.pot += myPot;
1024	–	skdat.f1 += dVdr;
1025	–
1002		if (i_is_Dipole)
1003	<	skdat.t1 -= cross(uz_i, duduz_i);
1003	>	*(idat.t1) -= cross(uz_i, indirect_duduz_i);
1004		if (j_is_Dipole)
1005	<	skdat.t2 -= cross(uz_j, duduz_j);
1005	>	*(idat.t2) -= cross(uz_j, indirect_duduz_j);
1006		}
1007	<	}
1007	>
1008	>
1009	>	return;
1010	>	}
1011
1012	<	void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1012	>	void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1013		RealType mu1, preVal, chg1, self;
1014
1015		if (!initialized_) initialize();
1016	<
1017	<	ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1016	>
1017	>	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1018
1019		// logicals
1041	–
1020		bool i_is_Charge = data.is_Charge;
1021		bool i_is_Dipole = data.is_Dipole;
1022
#	Line 1046 \| Line 1024 \| namespace OpenMD {
1024		if (i_is_Dipole) {
1025		mu1 = data.dipole_moment;
1026		preVal = pre22_ * preRF2_ * mu1 * mu1;
1027	<	scdat.pot -= 0.5 * preVal;
1027	>	((sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 preVal;
1028
1029		// The self-correction term adds into the reaction field vector
1030	<	Vector3d uz_i = scdat.eFrame.getColumn(2);
1030	>	Vector3d uz_i = sdat.eFrame->getColumn(2);
1031		Vector3d ei = preVal * uz_i;
1032
1033		// This looks very wrong. A vector crossed with itself is zero.
1034	<	scdat.t -= cross(uz_i, ei);
1034	>	*(sdat.t) -= cross(uz_i, ei);
1035		}
1036		} else if (summationMethod_ == esm_SHIFTED_FORCE \|\| summationMethod_ == esm_SHIFTED_POTENTIAL) {
1037		if (i_is_Charge) {
1038		chg1 = data.charge;
1039		if (screeningMethod_ == DAMPED) {
1040	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1040	>	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1041		} else {
1042	<	self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1042	>	self = - 0.5 * rcuti_ * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1043		}
1044	<	scdat.pot += self;
1044	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1045		}
1046		}
1047		}
1048
1049	<	RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1049	>	RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
1050		// This seems to work moderately well as a default. There's no
1051		// inherent scale for 1/r interactions that we can standardize.
1052		// 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 1613 by gezelter, Thu Aug 18 20:18:19 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 1613 by gezelter, Thu Aug 18 20:18:19 2011 UTC