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 1587 by gezelter, Fri Jul 8 20:25:32 2011 UTC vs.
Revision 1616 by gezelter, Tue Aug 30 15:45:35 2011 UTC

# Line 64 | Line 64 | namespace OpenMD {
64      Globals* simParams_ = info_->getSimParams();
65  
66      summationMap_["HARD"]               = esm_HARD;
67 +    summationMap_["NONE"]               = esm_HARD;
68      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
69      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
70      summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
# Line 116 | Line 117 | namespace OpenMD {
117          sprintf( painCave.errMsg,
118                   "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
119                   "\t(Input file specified %s .)\n"
120 <                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
120 >                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
121                   "\t\"shifted_potential\", \"shifted_force\", or \n"
122                   "\t\"reaction_field\".\n", myMethod.c_str() );
123          painCave.isFatal = 1;
# Line 249 | Line 250 | namespace OpenMD {
250        preRF2_ = 2.0 * preRF_;
251      }
252      
253 <    RealType dx = cutoffRadius_ / RealType(np_ - 1);
253 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
254 >    // have charged atoms longer than the cutoffRadius away from each
255 >    // other.  Splining may not be the best choice here.  Direct calls
256 >    // to erfc might be preferrable.
257 >
258 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
259      RealType rval;
260      vector<RealType> rvals;
261      vector<RealType> yvals;
# Line 578 | Line 584 | namespace OpenMD {
584        if (j_is_Charge) {
585          if (screeningMethod_ == DAMPED) {
586            // assemble the damping variables
587 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
588 <          erfcVal = res.first;
589 <          derfcVal = res.second;
587 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
588 >          //erfcVal = res.first;
589 >          //derfcVal = res.second;
590 >
591 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
592 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
593 >
594            c1 = erfcVal * riji;
595            c2 = (-derfcVal + c1) * riji;
596          } else {
# Line 667 | Line 677 | namespace OpenMD {
677  
678            if (screeningMethod_ == DAMPED) {
679              // assemble the damping variables
680 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
681 <            erfcVal = res.first;
682 <            derfcVal = res.second;
680 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
681 >            //erfcVal = res.first;
682 >            //derfcVal = res.second;
683 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
684 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
685              c1 = erfcVal * ri;
686              c2 = (-derfcVal + c1) * ri;
687              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 704 | Line 716 | namespace OpenMD {
716            
717          if (screeningMethod_ == DAMPED) {
718            // assemble the damping variables
719 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
720 <          erfcVal = res.first;
721 <          derfcVal = res.second;
719 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
720 >          //erfcVal = res.first;
721 >          //derfcVal = res.second;
722 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
723 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
724            c1 = erfcVal * riji;
725            c2 = (-derfcVal + c1) * riji;
726            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 793 | Line 807 | namespace OpenMD {
807              
808            if (screeningMethod_ == DAMPED) {
809              // assemble the damping variables
810 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
811 <            erfcVal = res.first;
812 <            derfcVal = res.second;
810 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
811 >            //erfcVal = res.first;
812 >            //derfcVal = res.second;
813 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
814 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
815              c1 = erfcVal * ri;
816              c2 = (-derfcVal + c1) * ri;
817              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 872 | Line 888 | namespace OpenMD {
888            }
889            if (screeningMethod_ == DAMPED) {
890              // assemble damping variables
891 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
892 <            erfcVal = res.first;
893 <            derfcVal = res.second;
891 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
892 >            //erfcVal = res.first;
893 >            //derfcVal = res.second;
894 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
895 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
896              c1 = erfcVal * ri;
897              c2 = (-derfcVal + c1) * ri;
898              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 923 | Line 941 | namespace OpenMD {
941  
942          if (screeningMethod_ == DAMPED) {
943            // assemble the damping variables
944 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
945 <          erfcVal = res.first;
946 <          derfcVal = res.second;
944 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
945 >          //erfcVal = res.first;
946 >          //derfcVal = res.second;
947 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
948 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
949            c1 = erfcVal * riji;
950            c2 = (-derfcVal + c1) * riji;
951            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 990 | Line 1010 | namespace OpenMD {
1010  
1011        // only accumulate the forces and torques resulting from the
1012        // indirect reaction field terms.
1013 +
1014        *(idat.vpair) += indirect_vpair;
1015        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1016        *(idat.f1) += indirect_dVdr;
# Line 1003 | Line 1024 | namespace OpenMD {
1024  
1025      return;
1026    }  
1006 –
1007 –  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
1008 –
1009 –    if (!initialized_) initialize();
1010 –    
1011 –    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
1012 –    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
1013 –    
1014 –    // logicals
1015 –
1016 –    bool i_is_Charge = data1.is_Charge;
1017 –    bool i_is_Dipole = data1.is_Dipole;
1018 –
1019 –    bool j_is_Charge = data2.is_Charge;
1020 –    bool j_is_Dipole = data2.is_Dipole;
1021 –
1022 –    RealType q_i, q_j;
1023 –
1024 –    // The skippedCharge computation is needed by the real-space
1025 –    // cutoff methods (i.e. shifted force and shifted potential)
1026 –      
1027 –    if (i_is_Charge) {
1028 –      q_i = data1.charge;
1029 –      *(idat.skippedCharge2) += q_i;
1030 –    }
1027      
1032 –    if (j_is_Charge) {
1033 –      q_j = data2.charge;
1034 –      *(idat.skippedCharge1) += q_j;
1035 –    }
1036 –
1037 –    // the rest of this function should only be necessary for reaction field.
1038 –
1039 –    if (summationMethod_ == esm_REACTION_FIELD) {
1040 –      RealType riji, ri2, ri3;
1041 –      RealType mu_i, ct_i;
1042 –      RealType mu_j, ct_j;
1043 –      RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
1044 –      Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
1045 –
1046 –      // some variables we'll need independent of electrostatic type:
1047 –      
1048 –      riji = 1.0 /  *(idat.rij) ;
1049 –      rhat =  *(idat.d)  * riji;
1050 –
1051 –      if (i_is_Dipole) {
1052 –        mu_i = data1.dipole_moment;
1053 –        uz_i = idat.eFrame1->getColumn(2);      
1054 –        ct_i = dot(uz_i, rhat);
1055 –        duduz_i = V3Zero;
1056 –      }
1057 –            
1058 –      if (j_is_Dipole) {
1059 –        mu_j = data2.dipole_moment;
1060 –        uz_j = idat.eFrame2->getColumn(2);      
1061 –        ct_j = dot(uz_j, rhat);
1062 –        duduz_j = V3Zero;
1063 –      }
1064 –    
1065 –      if (i_is_Charge) {
1066 –        if (j_is_Charge) {
1067 –          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
1068 –          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
1069 –          vterm = preVal * rfVal;
1070 –          myPot +=  *(idat.sw)  * vterm;        
1071 –          dudr  =  *(idat.sw)  * preVal * 2.0 * rfVal * riji;        
1072 –          dVdr += dudr * rhat;
1073 –        }
1074 –        
1075 –        if (j_is_Dipole) {
1076 –          ri2 = riji * riji;
1077 –          ri3 = ri2 * riji;        
1078 –          pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
1079 –          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
1080 –          myPot +=  *(idat.sw)  * vterm;        
1081 –          dVdr += - *(idat.sw)  * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1082 –          duduz_j += - *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1083 –        }
1084 –      }
1085 –      if (i_is_Dipole) {
1086 –        if (j_is_Charge) {
1087 –          ri2 = riji * riji;
1088 –          ri3 = ri2 * riji;        
1089 –          pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1090 –          vterm = - pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1091 –          myPot +=  *(idat.sw)  * vterm;        
1092 –          dVdr +=  *(idat.sw)  * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1093 –          duduz_i +=  *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij));
1094 –        }
1095 –      }
1096 –      
1097 –      // accumulate the forces and torques resulting from the self term
1098 –      (*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1099 –      *(idat.f1) += dVdr;
1100 –      
1101 –      if (i_is_Dipole)
1102 –        *(idat.t1) -= cross(uz_i, duduz_i);
1103 –      if (j_is_Dipole)
1104 –        *(idat.t2) -= cross(uz_j, duduz_j);
1105 –    }
1106 –  }
1107 –    
1028    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1029      RealType mu1, preVal, chg1, self;
1030      

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines