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 1721 by gezelter, Thu May 24 14:17:42 2012 UTC vs.
Revision 1734 by jmichalk, Tue Jun 5 17:49:36 2012 UTC

# Line 217 | Line 217 | namespace OpenMD {
217          addType(at);
218      }
219      
220 –
220      cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
221      rcuti_ = 1.0 / cutoffRadius_;
222      rcuti2_ = rcuti_ * rcuti_;
# Line 284 | Line 283 | namespace OpenMD {
283      electrostaticAtomData.is_Dipole = false;
284      electrostaticAtomData.is_SplitDipole = false;
285      electrostaticAtomData.is_Quadrupole = false;
286 +    electrostaticAtomData.is_Fluctuating = false;
287  
288      FixedChargeAdapter fca = FixedChargeAdapter(atomType);
289  
# Line 321 | Line 321 | namespace OpenMD {
321        electrostaticAtomData.hardness = fqa.getHardness();
322        electrostaticAtomData.slaterN = fqa.getSlaterN();
323        electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
324 –    } else {
325 –      electrostaticAtomData.is_Fluctuating = false;
324      }
325  
326      pair<map<int,AtomType*>::iterator,bool> ret;    
# Line 365 | 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 449 | Line 447 | namespace OpenMD {
447      Vector3d indirect_dVdr(V3Zero);
448      Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
449  
450 +    RealType coulInt, vFluc1(0.0), vFluc2(0.0);
451      pair<RealType, RealType> res;
452      
453      // splines for coulomb integrals
# Line 525 | 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 592 | 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 622 | namespace OpenMD {
622            dudr  = -  *(idat.sw)  * preVal * c2;
623            
624          }
626 –
625          
626 <        if (i_is_Fluctuating) {
627 <          if (!idat.excluded)
628 <            *(idat.dVdFQ1) += *(idat.sw) * vterm / q_i;
629 <          else {
630 <            res = J1->getValueAndDerivativeAt( *(idat.rij) );
631 <            *(idat.dVdFQ1) += pre11_ * res.first * q_j;
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 = coulInt - (*(idat.sw) * vterm);
636 >          } else {
637 >            vFluc1 = 0.0;
638            }
639 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j;
640          }
641 +
642          if (j_is_Fluctuating) {
643 <          if (!idat.excluded)
644 <            *(idat.dVdFQ2) += *(idat.sw) * vterm / q_j;
645 <          else {
646 <            res = J2->getValueAndDerivativeAt( *(idat.rij) );
647 <            *(idat.dVdFQ2) += pre11_ * res.first * q_i;
643 >          if (idat.excluded) {
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 = coulInt - (*(idat.sw) * vterm);
648 >          } else {
649 >            vFluc2 = 0.0;
650            }
651 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i;
652          }
653 <
654 <        vpair += vterm;
646 <        epot +=  *(idat.sw)  * vterm;
647 <        dVdr += dudr * rhat;                
653 >          
654 >
655        }
656  
657        if (j_is_Dipole) {
# Line 717 | Line 724 | namespace OpenMD {
724            duduz_j += -preSw * pot_term * rhat;
725  
726          }
727 +        if (i_is_Fluctuating) {
728 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
729 +        }
730        }
731  
732        if (j_is_Quadrupole) {
# Line 769 | Line 779 | namespace OpenMD {
779          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
780          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
781          duduz_j += preSw * qzz_j * cz_j * rhatdot2;
782 +        if (i_is_Fluctuating) {
783 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
784 +        }
785 +
786        }
787      }
788      
# Line 844 | Line 858 | namespace OpenMD {
858            // calculate derivatives for the forces and torques
859            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
860            duduz_i += preSw * pot_term * rhat;
861 +        }
862 +
863 +        if (j_is_Fluctuating) {
864 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
865          }
866 +
867        }
868  
869        if (j_is_Dipole) {
# Line 995 | Line 1014 | namespace OpenMD {
1014          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1015          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1016          duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1017 +
1018 +        if (j_is_Fluctuating) {
1019 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1020 +        }
1021 +
1022        }
1023      }
1024  
# Line 1033 | Line 1057 | namespace OpenMD {
1057          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1058      }
1059  
1036 –
1060      return;
1061    }  
1062      
1063    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1064 <    RealType mu1, preVal, chg1, self;
1042 <    
1064 >    RealType mu1, preVal, self;
1065      if (!initialized_) initialize();
1066  
1067      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 1047 | 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 1063 | Line 1093 | namespace OpenMD {
1093        }
1094      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1095        if (i_is_Charge) {        
1066 –        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