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 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs.
Revision 1528 by gezelter, Fri Dec 17 20:11:05 2010 UTC

# Line 52 | Line 52 | namespace OpenMD {
52   namespace OpenMD {
53    
54    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55 <                                  forceField_(NULL) {}
55 >                                  forceField_(NULL), info_(NULL) {}
56    
57    void Electrostatic::initialize() {
58 +
59 +    summationMap_["HARD"]               = esm_HARD;
60 +    summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
61 +    summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
62 +    summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
63 +    summationMap_["REACTION_FIELD"]     = esm_REACTION_FIELD;    
64 +    summationMap_["EWALD_FULL"]         = esm_EWALD_FULL;        
65 +    summationMap_["EWALD_PME"]          = esm_EWALD_PME;        
66 +    summationMap_["EWALD_SPME"]         = esm_EWALD_SPME;        
67 +    screeningMap_["DAMPED"]             = DAMPED;
68 +    screeningMap_["UNDAMPED"]           = UNDAMPED;
69 +
70      // these prefactors convert the multipole interactions into kcal / mol
71      // all were computed assuming distances are measured in angstroms
72      // Charge-Charge, assuming charges are measured in electrons
# Line 79 | Line 91 | namespace OpenMD {
91      
92      // variables to handle different summation methods for long-range
93      // electrostatics:
94 <    summationMethod_ = NONE;    
94 >    summationMethod_ = esm_HARD;    
95      screeningMethod_ = UNDAMPED;
96      dielectric_ = 1.0;
97      one_third_ = 1.0 / 3.0;
98 <    haveDefaultCutoff_ = false;
98 >    haveCutoffRadius_ = false;
99      haveDampingAlpha_ = false;
100      haveDielectric_ = false;  
101      haveElectroSpline_ = false;
102    
103 +    // check the summation method:
104 +    if (simParams_->haveElectrostaticSummationMethod()) {
105 +      string myMethod = simParams_->getElectrostaticSummationMethod();
106 +      toUpper(myMethod);
107 +      map<string, ElectrostaticSummationMethod>::iterator i;
108 +      i = summationMap_.find(myMethod);
109 +      if ( i != summationMap_.end() ) {
110 +        summationMethod_ = (*i).second;
111 +      } else {
112 +        // throw error
113 +        sprintf( painCave.errMsg,
114 +                 "SimInfo error: Unknown electrostaticSummationMethod.\n"
115 +                 "\t(Input file specified %s .)\n"
116 +                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
117 +                 "\t\"shifted_potential\", \"shifted_force\", or \n"
118 +                 "\t\"reaction_field\".\n", myMethod.c_str() );
119 +        painCave.isFatal = 1;
120 +        simError();
121 +      }
122 +    } else {
123 +      // set ElectrostaticSummationMethod to the cutoffMethod:
124 +      if (simParams_->haveCutoffMethod()){
125 +        string myMethod = simParams_->getCutoffMethod();
126 +        toUpper(myMethod);
127 +        map<string, ElectrostaticSummationMethod>::iterator i;
128 +        i = summationMap_.find(myMethod);
129 +        if ( i != summationMap_.end() ) {
130 +          summationMethod_ = (*i).second;
131 +        }
132 +      }
133 +    }
134 +    
135 +    if (summationMethod_ == esm_REACTION_FIELD) {        
136 +      if (!simParams_->haveDielectric()) {
137 +        // throw warning
138 +        sprintf( painCave.errMsg,
139 +                 "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
140 +                 "\tA default value of %f will be used for the dielectric.\n", dielectric_);
141 +        painCave.isFatal = 0;
142 +        painCave.severity = OPENMD_INFO;
143 +        simError();
144 +      } else {
145 +        dielectric_ = simParams_->getDielectric();      
146 +      }
147 +      haveDielectric_ = true;
148 +    }
149 +    
150 +    if (simParams_->haveElectrostaticScreeningMethod()) {
151 +      string myScreen = simParams_->getElectrostaticScreeningMethod();
152 +      toUpper(myScreen);
153 +      map<string, ElectrostaticScreeningMethod>::iterator i;
154 +      i = screeningMap_.find(myScreen);
155 +      if ( i != screeningMap_.end()) {
156 +        screeningMethod_ = (*i).second;
157 +      } else {
158 +        sprintf( painCave.errMsg,
159 +                 "SimInfo error: Unknown electrostaticScreeningMethod.\n"
160 +                 "\t(Input file specified %s .)\n"
161 +                 "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
162 +                 "or \"damped\".\n", myScreen.c_str() );
163 +        painCave.isFatal = 1;
164 +        simError();
165 +      }
166 +    }
167 +
168 +    // check to make sure a cutoff value has been set:
169 +    if (!haveCutoffRadius_) {
170 +      sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
171 +               "Cutoff value!\n");
172 +      painCave.severity = OPENMD_ERROR;
173 +      painCave.isFatal = 1;
174 +      simError();
175 +    }
176 +          
177 +    if (screeningMethod_ == DAMPED) {      
178 +      if (!simParams_->haveDampingAlpha()) {
179 +        // first set a cutoff dependent alpha value
180 +        // we assume alpha depends linearly with rcut from 0 to 20.5 ang
181 +        dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
182 +        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
183 +        
184 +        // throw warning
185 +        sprintf( painCave.errMsg,
186 +                 "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
187 +                 "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
188 +                 dampingAlpha_, cutoffRadius_);
189 +        painCave.severity = OPENMD_INFO;
190 +        painCave.isFatal = 0;
191 +        simError();
192 +      } else {
193 +        dampingAlpha_ = simParams_->getDampingAlpha();
194 +      }
195 +      haveDampingAlpha_ = true;
196 +    }
197 +
198      // find all of the Electrostatic atom Types:
199      ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
200      ForceField::AtomTypeContainer::MapTypeIterator i;
201      AtomType* at;
202 <
202 >    
203      for (at = atomTypes->beginType(i); at != NULL;
204           at = atomTypes->nextType(i)) {
205        
# Line 100 | Line 207 | namespace OpenMD {
207          addType(at);
208      }
209      
103 –    // check to make sure a cutoff value has been set:
104 –    if (!haveDefaultCutoff_) {
105 –      sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
106 –               "Cutoff value!\n");
107 –      painCave.severity = OPENMD_ERROR;
108 –      painCave.isFatal = 1;
109 –      simError();
110 –    }
210  
211 <    defaultCutoff2_ = defaultCutoff_ * defaultCutoff_;
212 <    rcuti_ = 1.0 / defaultCutoff_;
211 >    cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
212 >    rcuti_ = 1.0 / cutoffRadius_;
213      rcuti2_ = rcuti_ * rcuti_;
214      rcuti3_ = rcuti2_ * rcuti_;
215      rcuti4_ = rcuti2_ * rcuti2_;
216  
217      if (screeningMethod_ == DAMPED) {
218 <      if (!haveDampingAlpha_) {
120 <        sprintf( painCave.errMsg, "Electrostatic::initialize has no "
121 <                 "DampingAlpha value!\n");
122 <        painCave.severity = OPENMD_ERROR;
123 <        painCave.isFatal = 1;
124 <        simError();
125 <      }
126 <
218 >      
219        alpha2_ = dampingAlpha_ * dampingAlpha_;
220        alpha4_ = alpha2_ * alpha2_;
221        alpha6_ = alpha4_ * alpha2_;
222        alpha8_ = alpha4_ * alpha4_;
223        
224 <      constEXP_ = exp(-alpha2_ * defaultCutoff2_);
224 >      constEXP_ = exp(-alpha2_ * cutoffRadius2_);
225        invRootPi_ = 0.56418958354775628695;
226        alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
227  
228 <      c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * rcuti_;
228 >      c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
229        c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
230        c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
231        c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_;
# Line 148 | Line 240 | namespace OpenMD {
240        c6c_ = 9.0 * c5c_ * rcuti2_;
241      }
242    
243 <    if (summationMethod_ == REACTION_FIELD) {
244 <      if (haveDielectric_) {
245 <        preRF_ = (dielectric_ - 1.0) /
246 <            ((2.0 * dielectric_ + 1.0) * defaultCutoff2_ * defaultCutoff_);
155 <        preRF2_ = 2.0 * preRF_;
156 <      } else {
157 <        sprintf( painCave.errMsg, "Electrostatic::initialize has no Dielectric"
158 <                 " value!\n");
159 <        painCave.severity = OPENMD_ERROR;
160 <        painCave.isFatal = 1;
161 <        simError();
162 <      }
243 >    if (summationMethod_ == esm_REACTION_FIELD) {
244 >      preRF_ = (dielectric_ - 1.0) /
245 >        ((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
246 >      preRF2_ = 2.0 * preRF_;
247      }
248 <                              
249 <    RealType dx = defaultCutoff_ / RealType(np_ - 1);
248 >    
249 >    RealType dx = cutoffRadius_ / RealType(np_ - 1);
250      RealType rval;
251      vector<RealType> rvals;
252      vector<RealType> yvals;
# Line 283 | Line 367 | namespace OpenMD {
367            simError();                  
368          }
369          
370 +        // Quadrupoles in OpenMD are set as the diagonal elements
371 +        // of the diagonalized traceless quadrupole moment tensor.
372 +        // The column vectors of the unitary matrix that diagonalizes
373 +        // the quadrupole moment tensor become the eFrame (or the
374 +        // electrostatic version of the body-fixed frame.
375 +
376          Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
377          if (v3dData == NULL) {
378            sprintf( painCave.errMsg,
# Line 317 | Line 407 | namespace OpenMD {
407    
408    void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
409                                                      RealType theRSW ) {
410 <    defaultCutoff_ = theECR;
411 <    rrf_ = defaultCutoff_;
410 >    cutoffRadius_ = theECR;
411 >    rrf_ = cutoffRadius_;
412      rt_ = theRSW;
413 <    haveDefaultCutoff_ = true;
413 >    haveCutoffRadius_ = true;
414    }
415    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
416      summationMethod_ = esm;
# Line 479 | Line 569 | namespace OpenMD {
569  
570          preVal = idat.electroMult * pre11_ * q_i * q_j;
571          
572 <        if (summationMethod_ == SHIFTED_POTENTIAL) {
572 >        if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
573            vterm = preVal * (c1 - c1c_);
574            dudr  = -idat.sw * preVal * c2;
575  
576 <        } else if (summationMethod_ == SHIFTED_FORCE)  {
577 <          vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) );
576 >        } else if (summationMethod_ == esm_SHIFTED_FORCE)  {
577 >          vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
578            dudr  = idat.sw * preVal * (c2c_ - c2);
579  
580 <        } else if (summationMethod_ == REACTION_FIELD) {
580 >        } else if (summationMethod_ == esm_REACTION_FIELD) {
581            rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
582            vterm = preVal * ( riji + rfVal );            
583            dudr  = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
# Line 510 | Line 600 | namespace OpenMD {
600          pref = idat.electroMult * pre12_ * q_i * mu_j;
601          preSw = idat.sw * pref;
602  
603 <        if (summationMethod_ == REACTION_FIELD) {
603 >        if (summationMethod_ == esm_REACTION_FIELD) {
604            ri2 = riji * riji;
605            ri3 = ri2 * riji;
606      
# Line 622 | Line 712 | namespace OpenMD {
712          pref = idat.electroMult * pre12_ * q_j * mu_i;
713          preSw = idat.sw * pref;
714  
715 <        if (summationMethod_ == REACTION_FIELD) {
715 >        if (summationMethod_ == esm_REACTION_FIELD) {
716  
717            ri2 = riji * riji;
718            ri3 = ri2 * riji;
# Line 684 | Line 774 | namespace OpenMD {
774          pref = idat.electroMult * pre22_ * mu_i * mu_j;
775          preSw = idat.sw * pref;
776  
777 <        if (summationMethod_ == REACTION_FIELD) {
777 >        if (summationMethod_ == esm_REACTION_FIELD) {
778            ri2 = riji * riji;
779            ri3 = ri2 * riji;
780            ri4 = ri2 * ri2;
# Line 869 | Line 959 | namespace OpenMD {
959  
960      // the rest of this function should only be necessary for reaction field.
961  
962 <    if (summationMethod_ == REACTION_FIELD) {
962 >    if (summationMethod_ == esm_REACTION_FIELD) {
963        RealType riji, ri2, ri3;
964        RealType q_i, mu_i, ct_i;
965        RealType q_j, mu_j, ct_j;
# Line 950 | Line 1040 | namespace OpenMD {
1040      bool i_is_Charge = data.is_Charge;
1041      bool i_is_Dipole = data.is_Dipole;
1042  
1043 <    if (summationMethod_ == REACTION_FIELD) {
1043 >    if (summationMethod_ == esm_REACTION_FIELD) {
1044        if (i_is_Dipole) {
1045          mu1 = data.dipole_moment;          
1046          preVal = pre22_ * preRF2_ * mu1 * mu1;
# Line 963 | Line 1053 | namespace OpenMD {
1053          // This looks very wrong.  A vector crossed with itself is zero.
1054          scdat.t -= cross(uz_i, ei);
1055        }
1056 <    } else if (summationMethod_ == SHIFTED_FORCE || summationMethod_ == SHIFTED_POTENTIAL) {
1056 >    } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1057        if (i_is_Charge) {        
1058          chg1 = data.charge;
1059          if (screeningMethod_ == DAMPED) {
# Line 975 | Line 1065 | namespace OpenMD {
1065        }
1066      }
1067    }
1068 +
1069 +  RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1070 +    // This seems to work moderately well as a default.  There's no
1071 +    // inherent scale for 1/r interactions that we can standardize.
1072 +    // 12 angstroms seems to be a reasonably good guess for most
1073 +    // cases.
1074 +    return 12.0;
1075 +  }
1076   }

Diff Legend

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