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 1734 by jmichalk, Tue Jun 5 17:49:36 2012 UTC

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

Diff Legend

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