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 1720 by gezelter, Thu May 24 01:48:29 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 363 | 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 ) );
368 <          J2vals.push_back( sSTOCoulInt( b, a, n, m, 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(eaData2.hardness *  sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
370          }
371  
372          CubicSpline* J1 = new CubicSpline();
# Line 446 | 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
455 +    CubicSpline* J1;
456 +    CubicSpline* J2;
457 +    
458      if (!initialized_) initialize();
459      
460      ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
# Line 464 | Line 471 | namespace OpenMD {
471      bool i_is_Dipole = data1.is_Dipole;
472      bool i_is_SplitDipole = data1.is_SplitDipole;
473      bool i_is_Quadrupole = data1.is_Quadrupole;
474 +    bool i_is_Fluctuating = data1.is_Fluctuating;
475  
476      bool j_is_Charge = data2.is_Charge;
477      bool j_is_Dipole = data2.is_Dipole;
478      bool j_is_SplitDipole = data2.is_SplitDipole;
479      bool j_is_Quadrupole = data2.is_Quadrupole;
480 +    bool j_is_Fluctuating = data2.is_Fluctuating;
481      
482      if (i_is_Charge) {
483        q_i = data1.fixedCharge;
484 +
485 +      if (i_is_Fluctuating) {
486 +        q_i += *(idat.flucQ1);
487 +      }
488 +      
489        if (idat.excluded) {
490          *(idat.skippedCharge2) += q_i;
491        }
# Line 510 | Line 524 | namespace OpenMD {
524  
525      if (j_is_Charge) {
526        q_j = data2.fixedCharge;
527 +
528 +      if (j_is_Fluctuating)
529 +        q_j += *(idat.flucQ2);
530 +
531        if (idat.excluded) {
532          *(idat.skippedCharge1) += q_j;
533        }
# Line 547 | Line 565 | namespace OpenMD {
565        duduz_j = V3Zero;
566      }
567      
568 +    if (i_is_Fluctuating && j_is_Fluctuating) {
569 +      J1 = Jij[idat.atypes];
570 +      J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
571 +    }
572 +
573      epot = 0.0;
574      dVdr = V3Zero;
575      
# Line 569 | 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 598 | Line 621 | namespace OpenMD {
621  
622            vterm = preVal * riji * erfcVal;          
623            dudr  = -  *(idat.sw)  * preVal * c2;
624 +          
625 +        }
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 +            // 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 <        vpair += vterm;
644 <        epot +=  *(idat.sw)  * vterm;
645 <        dVdr += dudr * rhat;                
643 >        if (j_is_Fluctuating) {
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 >
656        }
657  
658        if (j_is_Dipole) {
# Line 676 | 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 728 | 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 804 | Line 860 | namespace OpenMD {
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 954 | 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 992 | Line 1058 | namespace OpenMD {
1058          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1059      }
1060  
995 –
1061      return;
1062    }  
1063      
1064    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1065 <    RealType mu1, preVal, chg1, self;
1001 <    
1065 >    RealType mu1, preVal, self;
1066      if (!initialized_) initialize();
1067  
1068      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 1006 | 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 1022 | Line 1094 | namespace OpenMD {
1094        }
1095      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1096        if (i_is_Charge) {        
1025 –        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