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 1710 by gezelter, Fri May 18 21:44:02 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  
58   namespace OpenMD {
59    
# Line 189 | 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 213 | Line 218 | namespace OpenMD {
218          addType(at);
219      }
220      
216 –
221      cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
222      rcuti_ = 1.0 / cutoffRadius_;
223      rcuti2_ = rcuti_ * rcuti_;
# Line 280 | 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 309 | Line 314 | namespace OpenMD {
314        }
315      }
316      
317 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
318 +
319 +    if (fqa.isFluctuatingCharge()) {
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;    
328      ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
# Line 322 | Line 336 | namespace OpenMD {
336        simError();        
337      }
338      
339 <    ElectrostaticMap[atomType] = electrostaticAtomData;    
339 >    ElectrostaticMap[atomType] = electrostaticAtomData;  
340 >
341 >    // Now, iterate over all known types and add to the mixing map:
342 >    
343 >    map<AtomType*, ElectrostaticAtomData>::iterator it;
344 >    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
345 >      AtomType* atype2 = (*it).first;
346 >      ElectrostaticAtomData eaData2 = (*it).second;
347 >      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
348 >        
349 >        RealType a = electrostaticAtomData.slaterZeta;
350 >        RealType b = eaData2.slaterZeta;
351 >        int m = electrostaticAtomData.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 >        // 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::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();
374 >        J1->addPoints(rvals, J1vals);
375 >        CubicSpline* J2 = new CubicSpline();
376 >        J2->addPoints(rvals, J2vals);
377 >        
378 >        pair<AtomType*, AtomType*> key1, key2;
379 >        key1 = make_pair(atomType, atype2);
380 >        key2 = make_pair(atype2, atomType);
381 >        
382 >        Jij[key1] = J1;
383 >        Jij[key2] = J2;
384 >      }
385 >    }
386 >
387      return;
388    }
389    
# Line 388 | 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 406 | 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 451 | 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 489 | 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 497 | 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 511 | 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 540 | 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 618 | 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 670 | 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 746 | Line 861 | namespace OpenMD {
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 896 | 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 925 | 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 934 | Line 1062 | namespace OpenMD {
1062          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1063      }
1064  
937 –
1065      return;
1066    }  
1067      
1068    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1069 <    RealType mu1, preVal, chg1, self;
943 <    
1069 >    RealType mu1, preVal, self;
1070      if (!initialized_) initialize();
1071  
1072      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 948 | 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 964 | Line 1100 | namespace OpenMD {
1100        }
1101      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1102        if (i_is_Charge) {        
967 –        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