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

# Line 48 | Line 48
48   #include "utils/simError.h"
49   #include "types/NonBondedInteractionType.hpp"
50   #include "types/FixedChargeAdapter.hpp"
51 + #include "types/FluctuatingChargeAdapter.hpp"
52   #include "types/MultipoleAdapter.hpp"
53   #include "io/Globals.hpp"
54   #include "nonbonded/SlaterIntegrals.hpp"
55   #include "utils/PhysicalConstants.hpp"
56 + #include "math/erfc.hpp"
57  
56 –
58   namespace OpenMD {
59    
60    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 192 | 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 216 | Line 218 | namespace OpenMD {
218          addType(at);
219      }
220      
219 –
221      cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
222      rcuti_ = 1.0 / cutoffRadius_;
223      rcuti2_ = rcuti_ * rcuti_;
# Line 283 | 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  
291      if (fca.isFixedCharge()) {
292        electrostaticAtomData.is_Charge = true;
293 <      electrostaticAtomData.charge = fca.getCharge();
293 >      electrostaticAtomData.fixedCharge = fca.getCharge();
294      }
295  
296      MultipoleAdapter ma = MultipoleAdapter(atomType);
# Line 315 | Line 317 | namespace OpenMD {
317      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
318  
319      if (fqa.isFluctuatingCharge()) {
320 <      electrostaticAtomData.is_FluctuatingCharge = true;
321 <      electrostaticAtomData.electronegativity = fca.getElectronegativity();
322 <      electrostaticAtomData.hardness = fca.getHardness();
323 <      electrostaticAtomData.slaterN = fca.getSlaterN();
324 <      electrostaticAtomData.slaterZeta = fca.getSlaterZeta();
320 >      electrostaticAtomData.is_Fluctuating = true;
321 >      electrostaticAtomData.electronegativity = fqa.getElectronegativity();
322 >      electrostaticAtomData.hardness = fqa.getHardness();
323 >      electrostaticAtomData.slaterN = fqa.getSlaterN();
324 >      electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
325      }
326  
327      pair<map<int,AtomType*>::iterator,bool> ret;    
# Line 341 | Line 343 | namespace OpenMD {
343      map<AtomType*, ElectrostaticAtomData>::iterator it;
344      for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
345        AtomType* atype2 = (*it).first;
346 <      
347 <      if ((*it).is_FluctuatingCharge && electrostaticAtomData.is_FluctuatingCharge) {
346 >      ElectrostaticAtomData eaData2 = (*it).second;
347 >      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
348          
349          RealType a = electrostaticAtomData.slaterZeta;
350 <        RealType b = (*it).slaterZeta;
350 >        RealType b = eaData2.slaterZeta;
351          int m = electrostaticAtomData.slaterN;
352 <        int n = (*it).slaterN;
352 >        int n = eaData2.slaterN;
353  
354          // Create the spline of the coulombic integral for s-type
355          // Slater orbitals.  Add a 2 angstrom safety window to deal
356          // with cutoffGroups that have charged atoms longer than the
357          // cutoffRadius away from each other.
358  
359 +        RealType rval;
360          RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
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 ) );
369 <          J2vals.push_back( sSTOCoulInt( b, a, n, m, 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::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
371          }
372  
373 <        CubicSpline J1 = new CubicSpline();
373 >        CubicSpline* J1 = new CubicSpline();
374          J1->addPoints(rvals, J1vals);
375 <        CubicSpline J2 = new CubicSpline();
375 >        CubicSpline* J2 = new CubicSpline();
376          J2->addPoints(rvals, J2vals);
377          
378          pair<AtomType*, AtomType*> key1, key2;
# Line 444 | Line 449 | namespace OpenMD {
449      Vector3d indirect_dVdr(V3Zero);
450      Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
451  
452 +    RealType coulInt, vFluc1(0.0), vFluc2(0.0);
453      pair<RealType, RealType> res;
454      
455 +    // splines for coulomb integrals
456 +    CubicSpline* J1;
457 +    CubicSpline* J2;
458 +    
459      if (!initialized_) initialize();
460      
461      ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
# Line 462 | Line 472 | namespace OpenMD {
472      bool i_is_Dipole = data1.is_Dipole;
473      bool i_is_SplitDipole = data1.is_SplitDipole;
474      bool i_is_Quadrupole = data1.is_Quadrupole;
475 +    bool i_is_Fluctuating = data1.is_Fluctuating;
476  
477      bool j_is_Charge = data2.is_Charge;
478      bool j_is_Dipole = data2.is_Dipole;
479      bool j_is_SplitDipole = data2.is_SplitDipole;
480      bool j_is_Quadrupole = data2.is_Quadrupole;
481 +    bool j_is_Fluctuating = data2.is_Fluctuating;
482      
483      if (i_is_Charge) {
484 <      q_i = data1.charge;
484 >      q_i = data1.fixedCharge;
485 >
486 >      if (i_is_Fluctuating) {
487 >        q_i += *(idat.flucQ1);
488 >      }
489 >      
490        if (idat.excluded) {
491          *(idat.skippedCharge2) += q_i;
492        }
# Line 507 | Line 524 | namespace OpenMD {
524      }
525  
526      if (j_is_Charge) {
527 <      q_j = data2.charge;
527 >      q_j = data2.fixedCharge;
528 >
529 >      if (j_is_Fluctuating)
530 >        q_j += *(idat.flucQ2);
531 >
532        if (idat.excluded) {
533          *(idat.skippedCharge1) += q_j;
534        }
# Line 545 | Line 566 | namespace OpenMD {
566        duduz_j = V3Zero;
567      }
568      
569 +    if (i_is_Fluctuating && j_is_Fluctuating) {
570 +      J1 = Jij[idat.atypes];
571 +      J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
572 +    }
573 +
574      epot = 0.0;
575      dVdr = V3Zero;
576      
# Line 553 | 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 567 | 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 596 | Line 622 | namespace OpenMD {
622  
623            vterm = preVal * riji * erfcVal;          
624            dudr  = -  *(idat.sw)  * preVal * c2;
625 +          
626 +        }
627 +        
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 = coulInt - (*(idat.sw) * vterm);
638 +          } else {
639 +            vFluc1 = 0.0;
640 +          }
641 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j;
642          }
643  
644 <        vpair += vterm;
645 <        epot +=  *(idat.sw)  * vterm;
646 <        dVdr += dudr * rhat;                
644 >        if (j_is_Fluctuating) {
645 >          if (idat.excluded) {
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 = coulInt - (*(idat.sw) * vterm);
650 >          } else {
651 >            vFluc2 = 0.0;
652 >          }
653 >          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i;
654 >        }
655 >          
656 >
657        }
658  
659        if (j_is_Dipole) {
# Line 674 | Line 726 | namespace OpenMD {
726            duduz_j += -preSw * pot_term * rhat;
727  
728          }
729 +        if (i_is_Fluctuating) {
730 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
731 +        }
732        }
733  
734        if (j_is_Quadrupole) {
# Line 726 | Line 781 | namespace OpenMD {
781          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
782          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
783          duduz_j += preSw * qzz_j * cz_j * rhatdot2;
784 +        if (i_is_Fluctuating) {
785 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
786 +        }
787 +
788        }
789      }
790      
# Line 801 | Line 860 | namespace OpenMD {
860            // calculate derivatives for the forces and torques
861            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
862            duduz_i += preSw * pot_term * rhat;
863 +        }
864 +
865 +        if (j_is_Fluctuating) {
866 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
867          }
868 +
869        }
870  
871        if (j_is_Dipole) {
# Line 952 | Line 1016 | namespace OpenMD {
1016          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1017          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1018          duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1019 +
1020 +        if (j_is_Fluctuating) {
1021 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1022 +        }
1023 +
1024        }
1025      }
1026  
# Line 981 | 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 990 | Line 1062 | namespace OpenMD {
1062          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1063      }
1064  
993 –
1065      return;
1066    }  
1067      
1068    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1069 <    RealType mu1, preVal, chg1, self;
999 <    
1069 >    RealType mu1, preVal, self;
1070      if (!initialized_) initialize();
1071  
1072      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 1004 | 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 1020 | Line 1100 | namespace OpenMD {
1100        }
1101      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1102        if (i_is_Charge) {        
1023 –        chg1 = data.charge;
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