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 1665 by gezelter, Tue Nov 22 20:38:56 2011 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 578 | Line 585 | namespace OpenMD {
585        if (j_is_Charge) {
586          if (screeningMethod_ == DAMPED) {
587            // assemble the damping variables
588 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
589 <          erfcVal = res.first;
590 <          derfcVal = res.second;
588 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
589 >          //erfcVal = res.first;
590 >          //derfcVal = res.second;
591 >
592 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
593 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
594 >
595            c1 = erfcVal * riji;
596            c2 = (-derfcVal + c1) * riji;
597          } else {
# Line 667 | Line 678 | namespace OpenMD {
678  
679            if (screeningMethod_ == DAMPED) {
680              // assemble the damping variables
681 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
682 <            erfcVal = res.first;
683 <            derfcVal = res.second;
681 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
682 >            //erfcVal = res.first;
683 >            //derfcVal = res.second;
684 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
685 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
686              c1 = erfcVal * ri;
687              c2 = (-derfcVal + c1) * ri;
688              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 704 | Line 717 | namespace OpenMD {
717            
718          if (screeningMethod_ == DAMPED) {
719            // assemble the damping variables
720 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
721 <          erfcVal = res.first;
722 <          derfcVal = res.second;
720 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
721 >          //erfcVal = res.first;
722 >          //derfcVal = res.second;
723 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
724 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
725            c1 = erfcVal * riji;
726            c2 = (-derfcVal + c1) * riji;
727            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 793 | Line 808 | namespace OpenMD {
808              
809            if (screeningMethod_ == DAMPED) {
810              // assemble the damping variables
811 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
812 <            erfcVal = res.first;
813 <            derfcVal = res.second;
811 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
812 >            //erfcVal = res.first;
813 >            //derfcVal = res.second;
814 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
815 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
816              c1 = erfcVal * ri;
817              c2 = (-derfcVal + c1) * ri;
818              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 872 | Line 889 | namespace OpenMD {
889            }
890            if (screeningMethod_ == DAMPED) {
891              // assemble damping variables
892 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
893 <            erfcVal = res.first;
894 <            derfcVal = res.second;
892 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
893 >            //erfcVal = res.first;
894 >            //derfcVal = res.second;
895 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
896 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
897              c1 = erfcVal * ri;
898              c2 = (-derfcVal + c1) * ri;
899              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 923 | Line 942 | namespace OpenMD {
942  
943          if (screeningMethod_ == DAMPED) {
944            // assemble the damping variables
945 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
946 <          erfcVal = res.first;
947 <          derfcVal = res.second;
945 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
946 >          //erfcVal = res.first;
947 >          //derfcVal = res.second;
948 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
949 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
950            c1 = erfcVal * riji;
951            c2 = (-derfcVal + c1) * riji;
952            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 990 | Line 1011 | namespace OpenMD {
1011  
1012        // only accumulate the forces and torques resulting from the
1013        // indirect reaction field terms.
1014 +
1015        *(idat.vpair) += indirect_vpair;
1016        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1017        *(idat.f1) += indirect_dVdr;

Diff Legend

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