ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/Electrostatic.cpp
(Generate patch)

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1601 by gezelter, Thu Aug 4 20:04:35 2011 UTC vs.
Revision 1668 by gezelter, Fri Jan 6 19:03:05 2012 UTC

# Line 36 | Line 36
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
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 64 | Line 65 | namespace OpenMD {
65      Globals* simParams_ = info_->getSimParams();
66  
67      summationMap_["HARD"]               = esm_HARD;
68 +    summationMap_["NONE"]               = esm_HARD;
69      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
70      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
71      summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
# Line 116 | Line 118 | namespace OpenMD {
118          sprintf( painCave.errMsg,
119                   "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
120                   "\t(Input file specified %s .)\n"
121 <                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
121 >                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
122                   "\t\"shifted_potential\", \"shifted_force\", or \n"
123                   "\t\"reaction_field\".\n", myMethod.c_str() );
124          painCave.isFatal = 1;
# Line 249 | Line 251 | namespace OpenMD {
251        preRF2_ = 2.0 * preRF_;
252      }
253      
254 <    RealType dx = cutoffRadius_ / RealType(np_ - 1);
254 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
255 >    // have charged atoms longer than the cutoffRadius away from each
256 >    // other.  Splining may not be the best choice here.  Direct calls
257 >    // to erfc might be preferrable.
258 >
259 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
260      RealType rval;
261      vector<RealType> rvals;
262      vector<RealType> yvals;
# Line 454 | Line 461 | namespace OpenMD {
461      RealType c1, c2, c3, c4;
462      RealType erfcVal(1.0), derfcVal(0.0);
463      RealType BigR;
464 +    RealType two(2.0), three(3.0);
465  
466      Vector3d Q_i, Q_j;
467      Vector3d ux_i, uy_i, uz_i;
# Line 578 | Line 586 | namespace OpenMD {
586        if (j_is_Charge) {
587          if (screeningMethod_ == DAMPED) {
588            // assemble the damping variables
589 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
590 <          erfcVal = res.first;
591 <          derfcVal = res.second;
589 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
590 >          //erfcVal = res.first;
591 >          //derfcVal = res.second;
592 >
593 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
594 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
595 >
596            c1 = erfcVal * riji;
597            c2 = (-derfcVal + c1) * riji;
598          } else {
# Line 610 | Line 622 | namespace OpenMD {
622            if (idat.excluded) {
623              indirect_vpair += preVal * rfVal;
624              indirect_Pot += *(idat.sw) * preVal * rfVal;
625 <            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
625 >            indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
626            }
627            
628          } else {
# Line 638 | Line 650 | namespace OpenMD {
650            vpair += vterm;
651            epot +=  *(idat.sw)  * vterm;
652  
653 <          dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
653 >          dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
654            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
655  
656            // Even if we excluded this pair from direct interactions,
# Line 667 | Line 679 | namespace OpenMD {
679  
680            if (screeningMethod_ == DAMPED) {
681              // assemble the damping variables
682 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
683 <            erfcVal = res.first;
684 <            derfcVal = res.second;
682 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
683 >            //erfcVal = res.first;
684 >            //derfcVal = res.second;
685 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
686 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
687              c1 = erfcVal * ri;
688              c2 = (-derfcVal + c1) * ri;
689              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 704 | Line 718 | namespace OpenMD {
718            
719          if (screeningMethod_ == DAMPED) {
720            // assemble the damping variables
721 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
722 <          erfcVal = res.first;
723 <          derfcVal = res.second;
721 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
722 >          //erfcVal = res.first;
723 >          //derfcVal = res.second;
724 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
725 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
726            c1 = erfcVal * riji;
727            c2 = (-derfcVal + c1) * riji;
728            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 723 | Line 739 | namespace OpenMD {
739          c2ri = c2 * riji;
740          c3ri = c3 * riji;
741          c4rij = c4 *  *(idat.rij) ;
742 <        rhatdot2 = 2.0 * rhat * c3;
742 >        rhatdot2 = two * rhat * c3;
743          rhatc4 = rhat * c4rij;
744  
745          // calculate the potential
# Line 736 | Line 752 | namespace OpenMD {
752                  
753          // calculate derivatives for the forces and torques
754  
755 <        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
756 <                           qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
757 <                           qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
755 >        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
756 >                           qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
757 >                           qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
758                            
759          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
760          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
# Line 762 | Line 778 | namespace OpenMD {
778            vpair += vterm;
779            epot +=  *(idat.sw)  * vterm;
780            
781 <          dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
781 >          dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
782            
783            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
784  
# Line 793 | Line 809 | namespace OpenMD {
809              
810            if (screeningMethod_ == DAMPED) {
811              // assemble the damping variables
812 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
813 <            erfcVal = res.first;
814 <            derfcVal = res.second;
812 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
813 >            //erfcVal = res.first;
814 >            //derfcVal = res.second;
815 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
816 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
817              c1 = erfcVal * ri;
818              c2 = (-derfcVal + c1) * ri;
819              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 838 | Line 856 | namespace OpenMD {
856              
857            a1 = 5.0 * ct_i * ct_j - ct_ij;
858              
859 <          dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
859 >          dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
860  
861 <          duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
862 <          duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
861 >          duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
862 >          duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
863  
864            if (idat.excluded) {
865              indirect_vpair +=  - pref * preRF2_ * ct_ij;
# Line 872 | Line 890 | namespace OpenMD {
890            }
891            if (screeningMethod_ == DAMPED) {
892              // assemble damping variables
893 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
894 <            erfcVal = res.first;
895 <            derfcVal = res.second;
893 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
894 >            //erfcVal = res.first;
895 >            //derfcVal = res.second;
896 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
897 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
898              c1 = erfcVal * ri;
899              c2 = (-derfcVal + c1) * ri;
900              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 923 | Line 943 | namespace OpenMD {
943  
944          if (screeningMethod_ == DAMPED) {
945            // assemble the damping variables
946 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
947 <          erfcVal = res.first;
948 <          derfcVal = res.second;
946 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
947 >          //erfcVal = res.first;
948 >          //derfcVal = res.second;
949 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
950 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
951            c1 = erfcVal * riji;
952            c2 = (-derfcVal + c1) * riji;
953            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 942 | Line 964 | namespace OpenMD {
964          c2ri = c2 * riji;
965          c3ri = c3 * riji;
966          c4rij = c4 *  *(idat.rij) ;
967 <        rhatdot2 = 2.0 * rhat * c3;
967 >        rhatdot2 = two * rhat * c3;
968          rhatc4 = rhat * c4rij;
969  
970          // calculate the potential
# Line 956 | Line 978 | namespace OpenMD {
978  
979          // calculate the derivatives for the forces and torques
980  
981 <        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
982 <                          qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
983 <                          qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
981 >        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
982 >                          qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
983 >                          qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
984  
985          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
986          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
# Line 990 | Line 1012 | namespace OpenMD {
1012  
1013        // only accumulate the forces and torques resulting from the
1014        // indirect reaction field terms.
1015 +
1016        *(idat.vpair) += indirect_vpair;
1017        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1018        *(idat.f1) += indirect_dVdr;

Diff Legend

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