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 1750 by gezelter, Thu Jun 7 12:53:46 2012 UTC

# 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 363 | Line 364 | namespace OpenMD {
364          for (int i = 0; i < np_; i++) {
365            rval = RealType(i) * dr;
366            rvals.push_back(rval);
367 <          J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
367 >          J1vals.push_back(electrostaticAtomData.hardness * sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
368            // may not be necessary if Slater coulomb integral is symmetric
369 <          J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
369 >          J2vals.push_back(eaData2.hardness *  sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
370          }
371  
372          CubicSpline* J1 = new CubicSpline();
# Line 524 | Line 525 | namespace OpenMD {
525      if (j_is_Charge) {
526        q_j = data2.fixedCharge;
527  
528 <      if (i_is_Fluctuating)
528 >      if (j_is_Fluctuating)
529          q_j += *(idat.flucQ2);
530  
531        if (idat.excluded) {
# Line 591 | Line 592 | namespace OpenMD {
592            c2 = c1 * riji;
593          }
594  
595 <        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
595 >        preVal =  *(idat.electroMult) * pre11_;
596          
597          if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
598            vterm = preVal * (c1 - c1c_);
# Line 623 | Line 624 | namespace OpenMD {
624            
625          }
626          
627 <        vpair += vterm;
628 <        epot +=  *(idat.sw)  * vterm;
629 <        dVdr += dudr * rhat;
627 >        vpair += vterm * q_i * q_j;
628 >        epot +=  *(idat.sw)  * vterm * q_i * q_j;
629 >        dVdr += dudr * rhat * q_i * q_j;
630  
631          if (i_is_Fluctuating) {
632            if (idat.excluded) {
633              // vFluc1 is the difference between the direct coulomb integral
634              // and the normal 1/r-like  interaction between point charges.
635              coulInt = J1->getValueAt( *(idat.rij) );
636 <            vFluc1 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
636 >            vFluc1 = coulInt - (*(idat.sw) * vterm);
637            } else {
638              vFluc1 = 0.0;
639            }
640 <          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) / q_i;
640 >          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j;
641          }
642  
643          if (j_is_Fluctuating) {
# Line 644 | Line 645 | namespace OpenMD {
645              // vFluc2 is the difference between the direct coulomb integral
646              // and the normal 1/r-like  interaction between point charges.
647              coulInt = J2->getValueAt( *(idat.rij) );
648 <            vFluc2 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
648 >            vFluc2 = coulInt - (*(idat.sw) * vterm);
649            } else {
650              vFluc2 = 0.0;
651            }
652 <          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) / q_j;
652 >          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i;
653          }
654            
655  
# Line 1061 | Line 1062 | namespace OpenMD {
1062    }  
1063      
1064    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1065 <    RealType mu1, preVal, chg1, self;
1065 <    
1065 >    RealType mu1, preVal, self;
1066      if (!initialized_) initialize();
1067  
1068      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 1070 | Line 1070 | namespace OpenMD {
1070      // logicals
1071      bool i_is_Charge = data.is_Charge;
1072      bool i_is_Dipole = data.is_Dipole;
1073 +    bool i_is_Fluctuating = data.is_Fluctuating;
1074 +    RealType chg1 = data.fixedCharge;  
1075 +    
1076 +    if (i_is_Fluctuating) {
1077 +      chg1 += *(sdat.flucQ);
1078 +      // dVdFQ is really a force, so this is negative the derivative
1079 +      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1080 +    }
1081  
1082      if (summationMethod_ == esm_REACTION_FIELD) {
1083        if (i_is_Dipole) {
# Line 1086 | Line 1094 | namespace OpenMD {
1094        }
1095      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1096        if (i_is_Charge) {        
1089 –        chg1 = data.fixedCharge;
1097          if (screeningMethod_ == DAMPED) {
1098            self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1099          } else {        

Diff Legend

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