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