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 1723 by gezelter, Thu May 24 20:59:54 2012 UTC vs.
Revision 1767 by gezelter, Fri Jul 6 22:01:58 2012 UTC

# Line 53 | Line 53
53   #include "io/Globals.hpp"
54   #include "nonbonded/SlaterIntegrals.hpp"
55   #include "utils/PhysicalConstants.hpp"
56 + #include "math/erfc.hpp"
57  
57 –
58   namespace OpenMD {
59    
60    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 193 | Line 193 | namespace OpenMD {
193          
194          // throw warning
195          sprintf( painCave.errMsg,
196 <                 "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
197 <                 "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
196 >                 "Electrostatic::initialize: dampingAlpha was not specified in the\n"
197 >                 "\tinput file.  A default value of %f (1/ang) will be used for the\n"
198 >                 "\tcutoff of %f (ang).\n",
199                   dampingAlpha_, cutoffRadius_);
200          painCave.severity = OPENMD_INFO;
201          painCave.isFatal = 0;
# Line 360 | Line 361 | namespace OpenMD {
361          vector<RealType> rvals;
362          vector<RealType> J1vals;
363          vector<RealType> J2vals;
364 <        for (int i = 0; i < np_; i++) {
364 >        // don't start at i = 0, as rval = 0 is undefined for the slater overlap integrals.
365 >        for (int i = 1; i < np_+1; i++) {
366            rval = RealType(i) * dr;
367            rvals.push_back(rval);
368 <          J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
368 >          J1vals.push_back(sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
369            // may not be necessary if Slater coulomb integral is symmetric
370 <          J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
370 >          J2vals.push_back(sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
371          }
372  
373          CubicSpline* J1 = new CubicSpline();
# Line 524 | Line 526 | namespace OpenMD {
526      if (j_is_Charge) {
527        q_j = data2.fixedCharge;
528  
529 <      if (i_is_Fluctuating)
529 >      if (j_is_Fluctuating)
530          q_j += *(idat.flucQ2);
531  
532        if (idat.excluded) {
# Line 577 | Line 579 | namespace OpenMD {
579        if (j_is_Charge) {
580          if (screeningMethod_ == DAMPED) {
581            // assemble the damping variables
582 <          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
583 <          //erfcVal = res.first;
584 <          //derfcVal = res.second;
582 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
583 >          erfcVal = res.first;
584 >          derfcVal = res.second;
585  
586 <          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
587 <          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
586 >          //erfcVal = erfc(dampingAlpha_ * *(idat.rij));
587 >          //derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
588  
589            c1 = erfcVal * riji;
590            c2 = (-derfcVal + c1) * riji;
# Line 591 | Line 593 | namespace OpenMD {
593            c2 = c1 * riji;
594          }
595  
596 <        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
596 >        preVal =  *(idat.electroMult) * pre11_;
597          
598          if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
599            vterm = preVal * (c1 - c1c_);
# Line 623 | Line 625 | namespace OpenMD {
625            
626          }
627          
628 <        vpair += vterm;
629 <        epot +=  *(idat.sw)  * vterm;
630 <        dVdr += dudr * rhat;
628 >        vpair += vterm * q_i * q_j;
629 >        epot +=  *(idat.sw)  * vterm * q_i * q_j;
630 >        dVdr += dudr * rhat * q_i * q_j;
631  
632          if (i_is_Fluctuating) {
633            if (idat.excluded) {
634              // vFluc1 is the difference between the direct coulomb integral
635              // and the normal 1/r-like  interaction between point charges.
636              coulInt = J1->getValueAt( *(idat.rij) );
637 <            vFluc1 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
637 >            vFluc1 = coulInt - (*(idat.sw) * vterm);
638            } else {
639              vFluc1 = 0.0;
640            }
641 <          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) / q_i;
641 >          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j;
642          }
643  
644          if (j_is_Fluctuating) {
# Line 644 | Line 646 | namespace OpenMD {
646              // vFluc2 is the difference between the direct coulomb integral
647              // and the normal 1/r-like  interaction between point charges.
648              coulInt = J2->getValueAt( *(idat.rij) );
649 <            vFluc2 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
649 >            vFluc2 = coulInt - (*(idat.sw) * vterm);
650            } else {
651              vFluc2 = 0.0;
652            }
653 <          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) / q_j;
653 >          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i;
654          }
655            
656  
# Line 1048 | Line 1050 | namespace OpenMD {
1050        // indirect reaction field terms.
1051  
1052        *(idat.vpair) += indirect_vpair;
1053 +      
1054 +      (*(idat.excludedPot))[ELECTROSTATIC_FAMILY] +=   (*(idat.sw) * vterm +
1055 +                                                        vFluc1 ) * q_i * q_j;
1056        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1057        *(idat.f1) += indirect_dVdr;
1058        
# Line 1061 | Line 1066 | namespace OpenMD {
1066    }  
1067      
1068    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1069 <    RealType mu1, preVal, chg1, self;
1065 <    
1069 >    RealType mu1, preVal, self;
1070      if (!initialized_) initialize();
1071  
1072      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 1070 | Line 1074 | namespace OpenMD {
1074      // logicals
1075      bool i_is_Charge = data.is_Charge;
1076      bool i_is_Dipole = data.is_Dipole;
1077 +    bool i_is_Fluctuating = data.is_Fluctuating;
1078 +    RealType chg1 = data.fixedCharge;  
1079 +    
1080 +    if (i_is_Fluctuating) {
1081 +      chg1 += *(sdat.flucQ);
1082 +      // dVdFQ is really a force, so this is negative the derivative
1083 +      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1084 +      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1085 +        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1086 +    }
1087  
1088      if (summationMethod_ == esm_REACTION_FIELD) {
1089        if (i_is_Dipole) {
# Line 1086 | Line 1100 | namespace OpenMD {
1100        }
1101      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1102        if (i_is_Charge) {        
1089 –        chg1 = data.fixedCharge;
1103          if (screeningMethod_ == DAMPED) {
1104            self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1105          } else {        

Diff Legend

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