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 1586 by gezelter, Tue Jun 21 06:34:35 2011 UTC

# Line 47 | Line 47
47   #include "utils/simError.h"
48   #include "types/NonBondedInteractionType.hpp"
49   #include "types/DirectionalAtomType.hpp"
50 + #include "io/Globals.hpp"
51  
51 –
52   namespace OpenMD {
53    
54    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55 <                                  forceField_(NULL) {}
55 >                                  forceField_(NULL), info_(NULL), haveCutoffRadius_(false),
56 >                                  haveDampingAlpha_(false), haveDielectric_(false),
57 >                                  haveElectroSpline_(false)
58 > {}
59    
60    void Electrostatic::initialize() {
61 +
62 +    Globals* simParams_ = info_->getSimParams();
63 +
64 +    summationMap_["HARD"]               = esm_HARD;
65 +    summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
66 +    summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
67 +    summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
68 +    summationMap_["REACTION_FIELD"]     = esm_REACTION_FIELD;    
69 +    summationMap_["EWALD_FULL"]         = esm_EWALD_FULL;        
70 +    summationMap_["EWALD_PME"]          = esm_EWALD_PME;        
71 +    summationMap_["EWALD_SPME"]         = esm_EWALD_SPME;        
72 +    screeningMap_["DAMPED"]             = DAMPED;
73 +    screeningMap_["UNDAMPED"]           = UNDAMPED;
74 +
75      // these prefactors convert the multipole interactions into kcal / mol
76      // all were computed assuming distances are measured in angstroms
77      // Charge-Charge, assuming charges are measured in electrons
# Line 79 | Line 96 | namespace OpenMD {
96      
97      // variables to handle different summation methods for long-range
98      // electrostatics:
99 <    summationMethod_ = NONE;    
99 >    summationMethod_ = esm_HARD;    
100      screeningMethod_ = UNDAMPED;
101      dielectric_ = 1.0;
102      one_third_ = 1.0 / 3.0;
86 –    haveDefaultCutoff_ = false;
87 –    haveDampingAlpha_ = false;
88 –    haveDielectric_ = false;  
89 –    haveElectroSpline_ = false;
103    
104 +    // check the summation method:
105 +    if (simParams_->haveElectrostaticSummationMethod()) {
106 +      string myMethod = simParams_->getElectrostaticSummationMethod();
107 +      toUpper(myMethod);
108 +      map<string, ElectrostaticSummationMethod>::iterator i;
109 +      i = summationMap_.find(myMethod);
110 +      if ( i != summationMap_.end() ) {
111 +        summationMethod_ = (*i).second;
112 +      } else {
113 +        // throw error
114 +        sprintf( painCave.errMsg,
115 +                 "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
116 +                 "\t(Input file specified %s .)\n"
117 +                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
118 +                 "\t\"shifted_potential\", \"shifted_force\", or \n"
119 +                 "\t\"reaction_field\".\n", myMethod.c_str() );
120 +        painCave.isFatal = 1;
121 +        simError();
122 +      }
123 +    } else {
124 +      // set ElectrostaticSummationMethod to the cutoffMethod:
125 +      if (simParams_->haveCutoffMethod()){
126 +        string myMethod = simParams_->getCutoffMethod();
127 +        toUpper(myMethod);
128 +        map<string, ElectrostaticSummationMethod>::iterator i;
129 +        i = summationMap_.find(myMethod);
130 +        if ( i != summationMap_.end() ) {
131 +          summationMethod_ = (*i).second;
132 +        }
133 +      }
134 +    }
135 +    
136 +    if (summationMethod_ == esm_REACTION_FIELD) {        
137 +      if (!simParams_->haveDielectric()) {
138 +        // throw warning
139 +        sprintf( painCave.errMsg,
140 +                 "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
141 +                 "\tA default value of %f will be used for the dielectric.\n", dielectric_);
142 +        painCave.isFatal = 0;
143 +        painCave.severity = OPENMD_INFO;
144 +        simError();
145 +      } else {
146 +        dielectric_ = simParams_->getDielectric();      
147 +      }
148 +      haveDielectric_ = true;
149 +    }
150 +    
151 +    if (simParams_->haveElectrostaticScreeningMethod()) {
152 +      string myScreen = simParams_->getElectrostaticScreeningMethod();
153 +      toUpper(myScreen);
154 +      map<string, ElectrostaticScreeningMethod>::iterator i;
155 +      i = screeningMap_.find(myScreen);
156 +      if ( i != screeningMap_.end()) {
157 +        screeningMethod_ = (*i).second;
158 +      } else {
159 +        sprintf( painCave.errMsg,
160 +                 "SimInfo error: Unknown electrostaticScreeningMethod.\n"
161 +                 "\t(Input file specified %s .)\n"
162 +                 "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
163 +                 "or \"damped\".\n", myScreen.c_str() );
164 +        painCave.isFatal = 1;
165 +        simError();
166 +      }
167 +    }
168 +
169 +    // check to make sure a cutoff value has been set:
170 +    if (!haveCutoffRadius_) {
171 +      sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
172 +               "Cutoff value!\n");
173 +      painCave.severity = OPENMD_ERROR;
174 +      painCave.isFatal = 1;
175 +      simError();
176 +    }
177 +          
178 +    if (screeningMethod_ == DAMPED) {      
179 +      if (!simParams_->haveDampingAlpha()) {
180 +        // first set a cutoff dependent alpha value
181 +        // we assume alpha depends linearly with rcut from 0 to 20.5 ang
182 +        dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
183 +        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
184 +        
185 +        // throw warning
186 +        sprintf( painCave.errMsg,
187 +                 "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
188 +                 "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
189 +                 dampingAlpha_, cutoffRadius_);
190 +        painCave.severity = OPENMD_INFO;
191 +        painCave.isFatal = 0;
192 +        simError();
193 +      } else {
194 +        dampingAlpha_ = simParams_->getDampingAlpha();
195 +      }
196 +      haveDampingAlpha_ = true;
197 +    }
198 +
199      // find all of the Electrostatic atom Types:
200      ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
201      ForceField::AtomTypeContainer::MapTypeIterator i;
202      AtomType* at;
203 <
203 >    
204      for (at = atomTypes->beginType(i); at != NULL;
205           at = atomTypes->nextType(i)) {
206        
# Line 100 | Line 208 | namespace OpenMD {
208          addType(at);
209      }
210      
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 –    }
211  
212 <    defaultCutoff2_ = defaultCutoff_ * defaultCutoff_;
213 <    rcuti_ = 1.0 / defaultCutoff_;
212 >    cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
213 >    rcuti_ = 1.0 / cutoffRadius_;
214      rcuti2_ = rcuti_ * rcuti_;
215      rcuti3_ = rcuti2_ * rcuti_;
216      rcuti4_ = rcuti2_ * rcuti2_;
217  
218      if (screeningMethod_ == DAMPED) {
219 <      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 <
219 >      
220        alpha2_ = dampingAlpha_ * dampingAlpha_;
221        alpha4_ = alpha2_ * alpha2_;
222        alpha6_ = alpha4_ * alpha2_;
223        alpha8_ = alpha4_ * alpha4_;
224        
225 <      constEXP_ = exp(-alpha2_ * defaultCutoff2_);
225 >      constEXP_ = exp(-alpha2_ * cutoffRadius2_);
226        invRootPi_ = 0.56418958354775628695;
227        alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
228  
229 <      c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * rcuti_;
229 >      c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
230        c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
231        c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
232        c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_;
# Line 148 | Line 241 | namespace OpenMD {
241        c6c_ = 9.0 * c5c_ * rcuti2_;
242      }
243    
244 <    if (summationMethod_ == REACTION_FIELD) {
245 <      if (haveDielectric_) {
246 <        preRF_ = (dielectric_ - 1.0) /
247 <            ((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 <      }
244 >    if (summationMethod_ == esm_REACTION_FIELD) {
245 >      preRF_ = (dielectric_ - 1.0) /
246 >        ((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
247 >      preRF2_ = 2.0 * preRF_;
248      }
249 <                              
250 <    RealType dx = defaultCutoff_ / RealType(np_ - 1);
249 >    
250 >    RealType dx = cutoffRadius_ / RealType(np_ - 1);
251      RealType rval;
252      vector<RealType> rvals;
253      vector<RealType> yvals;
# Line 283 | Line 368 | namespace OpenMD {
368            simError();                  
369          }
370          
371 +        // Quadrupoles in OpenMD are set as the diagonal elements
372 +        // of the diagonalized traceless quadrupole moment tensor.
373 +        // The column vectors of the unitary matrix that diagonalizes
374 +        // the quadrupole moment tensor become the eFrame (or the
375 +        // electrostatic version of the body-fixed frame.
376 +
377          Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
378          if (v3dData == NULL) {
379            sprintf( painCave.errMsg,
# Line 315 | Line 406 | namespace OpenMD {
406      return;
407    }
408    
409 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
410 <                                                    RealType theRSW ) {
411 <    defaultCutoff_ = theECR;
412 <    rrf_ = defaultCutoff_;
322 <    rt_ = theRSW;
323 <    haveDefaultCutoff_ = true;
409 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
410 >    cutoffRadius_ = rCut;
411 >    rrf_ = cutoffRadius_;
412 >    haveCutoffRadius_ = true;
413    }
414 +
415 +  void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
416 +    rt_ = rSwitch;
417 +  }
418    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
419      summationMethod_ = esm;
420    }
# Line 337 | Line 430 | namespace OpenMD {
430      haveDielectric_ = true;
431    }
432  
433 <  void Electrostatic::calcForce(InteractionData idat) {
433 >  void Electrostatic::calcForce(InteractionData &idat) {
434  
435      // utility variables.  Should clean these up and use the Vector3d and
436      // Mat3x3d to replace as many as we can in future versions:
# Line 371 | Line 464 | namespace OpenMD {
464      
465      if (!initialized_) initialize();
466      
467 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
468 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
467 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
468 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
469      
470      // some variables we'll need independent of electrostatic type:
471  
472 <    riji = 1.0 / idat.rij;
473 <    Vector3d rhat = idat.d  * riji;
472 >    riji = 1.0 /  *(idat.rij) ;
473 >    Vector3d rhat =  *(idat.d)   * riji;
474  
475      // logicals
476  
# Line 396 | Line 489 | namespace OpenMD {
489  
490      if (i_is_Dipole) {
491        mu_i = data1.dipole_moment;
492 <      uz_i = idat.eFrame1.getColumn(2);
492 >      uz_i = idat.eFrame1->getColumn(2);
493        
494        ct_i = dot(uz_i, rhat);
495  
# Line 412 | Line 505 | namespace OpenMD {
505        qyy_i = Q_i.y();
506        qzz_i = Q_i.z();
507        
508 <      ux_i = idat.eFrame1.getColumn(0);
509 <      uy_i = idat.eFrame1.getColumn(1);
510 <      uz_i = idat.eFrame1.getColumn(2);
508 >      ux_i = idat.eFrame1->getColumn(0);
509 >      uy_i = idat.eFrame1->getColumn(1);
510 >      uz_i = idat.eFrame1->getColumn(2);
511  
512        cx_i = dot(ux_i, rhat);
513        cy_i = dot(uy_i, rhat);
# Line 430 | Line 523 | namespace OpenMD {
523  
524      if (j_is_Dipole) {
525        mu_j = data2.dipole_moment;
526 <      uz_j = idat.eFrame2.getColumn(2);
526 >      uz_j = idat.eFrame2->getColumn(2);
527        
528        ct_j = dot(uz_j, rhat);
529  
# Line 446 | Line 539 | namespace OpenMD {
539        qyy_j = Q_j.y();
540        qzz_j = Q_j.z();
541        
542 <      ux_j = idat.eFrame2.getColumn(0);
543 <      uy_j = idat.eFrame2.getColumn(1);
544 <      uz_j = idat.eFrame2.getColumn(2);
542 >      ux_j = idat.eFrame2->getColumn(0);
543 >      uy_j = idat.eFrame2->getColumn(1);
544 >      uz_j = idat.eFrame2->getColumn(2);
545  
546        cx_j = dot(ux_j, rhat);
547        cy_j = dot(uy_j, rhat);
# Line 467 | Line 560 | namespace OpenMD {
560        if (j_is_Charge) {
561          if (screeningMethod_ == DAMPED) {
562            // assemble the damping variables
563 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
563 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
564            erfcVal = res.first;
565            derfcVal = res.second;
566            c1 = erfcVal * riji;
# Line 477 | Line 570 | namespace OpenMD {
570            c2 = c1 * riji;
571          }
572  
573 <        preVal = idat.electroMult * pre11_ * q_i * q_j;
573 >        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
574          
575 <        if (summationMethod_ == SHIFTED_POTENTIAL) {
575 >        if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
576            vterm = preVal * (c1 - c1c_);
577 <          dudr  = -idat.sw * preVal * c2;
577 >          dudr  = - *(idat.sw)  * preVal * c2;
578  
579 <        } else if (summationMethod_ == SHIFTED_FORCE)  {
580 <          vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) );
581 <          dudr  = idat.sw * preVal * (c2c_ - c2);
579 >        } else if (summationMethod_ == esm_SHIFTED_FORCE)  {
580 >          vterm = preVal * ( c1 - c1c_ + c2c_*( *(idat.rij)  - cutoffRadius_) );
581 >          dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
582  
583 <        } else if (summationMethod_ == REACTION_FIELD) {
584 <          rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
583 >        } else if (summationMethod_ == esm_REACTION_FIELD) {
584 >          rfVal =  *(idat.electroMult) * preRF_ *  *(idat.rij)  *  *(idat.rij) ;
585            vterm = preVal * ( riji + rfVal );            
586 <          dudr  = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
586 >          dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
587  
588          } else {
589            vterm = preVal * riji * erfcVal;            
590  
591 <          dudr  = - idat.sw * preVal * c2;
591 >          dudr  = -  *(idat.sw)  * preVal * c2;
592  
593          }
594  
595 <        idat.vpair += vterm;
596 <        epot += idat.sw * vterm;
595 >        *(idat.vpair) += vterm;
596 >        epot +=  *(idat.sw)  * vterm;
597  
598          dVdr += dudr * rhat;      
599        }
600  
601        if (j_is_Dipole) {
602          // pref is used by all the possible methods
603 <        pref = idat.electroMult * pre12_ * q_i * mu_j;
604 <        preSw = idat.sw * pref;
603 >        pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
604 >        preSw =  *(idat.sw)  * pref;
605  
606 <        if (summationMethod_ == REACTION_FIELD) {
606 >        if (summationMethod_ == esm_REACTION_FIELD) {
607            ri2 = riji * riji;
608            ri3 = ri2 * riji;
609      
610 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
611 <          idat.vpair += vterm;
612 <          epot += idat.sw * vterm;
610 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
611 >          *(idat.vpair) += vterm;
612 >          epot +=  *(idat.sw)  * vterm;
613  
614            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
615 <          duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);  
615 >          duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
616  
617          } else {
618            // determine the inverse r used if we have split dipoles
619            if (j_is_SplitDipole) {
620 <            BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
620 >            BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
621              ri = 1.0 / BigR;
622 <            scale = idat.rij * ri;
622 >            scale =  *(idat.rij)  * ri;
623            } else {
624              ri = riji;
625              scale = 1.0;
# Line 536 | Line 629 | namespace OpenMD {
629  
630            if (screeningMethod_ == DAMPED) {
631              // assemble the damping variables
632 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
632 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
633              erfcVal = res.first;
634              derfcVal = res.second;
635              c1 = erfcVal * ri;
# Line 553 | Line 646 | namespace OpenMD {
646            // calculate the potential
647            pot_term =  scale * c2;
648            vterm = -pref * ct_j * pot_term;
649 <          idat.vpair += vterm;
650 <          epot += idat.sw * vterm;
649 >          *(idat.vpair) += vterm;
650 >          epot +=  *(idat.sw)  * vterm;
651              
652            // calculate derivatives for forces and torques
653  
# Line 569 | Line 662 | namespace OpenMD {
662          cx2 = cx_j * cx_j;
663          cy2 = cy_j * cy_j;
664          cz2 = cz_j * cz_j;
665 <        pref =  idat.electroMult * pre14_ * q_i * one_third_;
665 >        pref =   *(idat.electroMult) * pre14_ * q_i * one_third_;
666            
667          if (screeningMethod_ == DAMPED) {
668            // assemble the damping variables
669 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
669 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
670            erfcVal = res.first;
671            derfcVal = res.second;
672            c1 = erfcVal * riji;
# Line 588 | Line 681 | namespace OpenMD {
681          }
682  
683          // precompute variables for convenience
684 <        preSw = idat.sw * pref;
684 >        preSw =  *(idat.sw)  * pref;
685          c2ri = c2 * riji;
686          c3ri = c3 * riji;
687 <        c4rij = c4 * idat.rij;
687 >        c4rij = c4 *  *(idat.rij) ;
688          rhatdot2 = 2.0 * rhat * c3;
689          rhatc4 = rhat * c4rij;
690  
# Line 600 | Line 693 | namespace OpenMD {
693                       qyy_j * (cy2*c3 - c2ri) +
694                       qzz_j * (cz2*c3 - c2ri) );
695          vterm = pref * pot_term;
696 <        idat.vpair += vterm;
697 <        epot += idat.sw * vterm;
696 >        *(idat.vpair) += vterm;
697 >        epot +=  *(idat.sw)  * vterm;
698                  
699          // calculate derivatives for the forces and torques
700  
# Line 619 | Line 712 | namespace OpenMD {
712  
713        if (j_is_Charge) {
714          // variables used by all the methods
715 <        pref = idat.electroMult * pre12_ * q_j * mu_i;
716 <        preSw = idat.sw * pref;
715 >        pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
716 >        preSw =  *(idat.sw)  * pref;
717  
718 <        if (summationMethod_ == REACTION_FIELD) {
718 >        if (summationMethod_ == esm_REACTION_FIELD) {
719  
720            ri2 = riji * riji;
721            ri3 = ri2 * riji;
722  
723 <          vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
724 <          idat.vpair += vterm;
725 <          epot += idat.sw * vterm;
723 >          vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
724 >          *(idat.vpair) += vterm;
725 >          epot +=  *(idat.sw)  * vterm;
726            
727            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
728            
729 <          duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
729 >          duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
730              
731          } else {
732            
733            // determine inverse r if we are using split dipoles
734            if (i_is_SplitDipole) {
735 <            BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
735 >            BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
736              ri = 1.0 / BigR;
737 <            scale = idat.rij * ri;
737 >            scale =  *(idat.rij)  * ri;
738            } else {
739              ri = riji;
740              scale = 1.0;
# Line 651 | Line 744 | namespace OpenMD {
744              
745            if (screeningMethod_ == DAMPED) {
746              // assemble the damping variables
747 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
747 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
748              erfcVal = res.first;
749              derfcVal = res.second;
750              c1 = erfcVal * ri;
# Line 668 | Line 761 | namespace OpenMD {
761            // calculate the potential
762            pot_term = c2 * scale;
763            vterm = pref * ct_i * pot_term;
764 <          idat.vpair += vterm;
765 <          epot += idat.sw * vterm;
764 >          *(idat.vpair) += vterm;
765 >          epot +=  *(idat.sw)  * vterm;
766  
767            // calculate derivatives for the forces and torques
768            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
# Line 681 | Line 774 | namespace OpenMD {
774          // variables used by all methods
775          ct_ij = dot(uz_i, uz_j);
776  
777 <        pref = idat.electroMult * pre22_ * mu_i * mu_j;
778 <        preSw = idat.sw * pref;
777 >        pref =  *(idat.electroMult) * pre22_ * mu_i * mu_j;
778 >        preSw =  *(idat.sw)  * pref;
779  
780 <        if (summationMethod_ == REACTION_FIELD) {
780 >        if (summationMethod_ == esm_REACTION_FIELD) {
781            ri2 = riji * riji;
782            ri3 = ri2 * riji;
783            ri4 = ri2 * ri2;
784  
785            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
786                             preRF2_ * ct_ij );
787 <          idat.vpair += vterm;
788 <          epot += idat.sw * vterm;
787 >          *(idat.vpair) += vterm;
788 >          epot +=  *(idat.sw)  * vterm;
789              
790            a1 = 5.0 * ct_i * ct_j - ct_ij;
791              
# Line 705 | Line 798 | namespace OpenMD {
798            
799            if (i_is_SplitDipole) {
800              if (j_is_SplitDipole) {
801 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
801 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
802              } else {
803 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
803 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
804              }
805              ri = 1.0 / BigR;
806 <            scale = idat.rij * ri;
806 >            scale =  *(idat.rij)  * ri;
807            } else {
808              if (j_is_SplitDipole) {
809 <              BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
809 >              BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
810                ri = 1.0 / BigR;
811 <              scale = idat.rij * ri;
811 >              scale =  *(idat.rij)  * ri;
812              } else {
813                ri = riji;
814                scale = 1.0;
# Line 723 | Line 816 | namespace OpenMD {
816            }
817            if (screeningMethod_ == DAMPED) {
818              // assemble damping variables
819 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
819 >            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
820              erfcVal = res.first;
821              derfcVal = res.second;
822              c1 = erfcVal * ri;
# Line 745 | Line 838 | namespace OpenMD {
838            preSwSc = preSw * scale;
839            c2ri = c2 * ri;
840            c3ri = c3 * ri;
841 <          c4rij = c4 * idat.rij;
841 >          c4rij = c4 *  *(idat.rij) ;
842  
843            // calculate the potential
844            pot_term = (ct_ij * c2ri - ctidotj * c3);
845            vterm = pref * pot_term;
846 <          idat.vpair += vterm;
847 <          epot += idat.sw * vterm;
846 >          *(idat.vpair) += vterm;
847 >          epot +=  *(idat.sw)  * vterm;
848  
849            // calculate derivatives for the forces and torques
850            dVdr += preSwSc * ( ctidotj * rhat * c4rij  -
# Line 770 | Line 863 | namespace OpenMD {
863          cy2 = cy_i * cy_i;
864          cz2 = cz_i * cz_i;
865  
866 <        pref = idat.electroMult * pre14_ * q_j * one_third_;
866 >        pref =  *(idat.electroMult) * pre14_ * q_j * one_third_;
867  
868          if (screeningMethod_ == DAMPED) {
869            // assemble the damping variables
870 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
870 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
871            erfcVal = res.first;
872            derfcVal = res.second;
873            c1 = erfcVal * riji;
# Line 789 | Line 882 | namespace OpenMD {
882          }
883            
884          // precompute some variables for convenience
885 <        preSw = idat.sw * pref;
885 >        preSw =  *(idat.sw)  * pref;
886          c2ri = c2 * riji;
887          c3ri = c3 * riji;
888 <        c4rij = c4 * idat.rij;
888 >        c4rij = c4 *  *(idat.rij) ;
889          rhatdot2 = 2.0 * rhat * c3;
890          rhatc4 = rhat * c4rij;
891  
# Line 802 | Line 895 | namespace OpenMD {
895                       qzz_i * (cz2 * c3 - c2ri) );
896          
897          vterm = pref * pot_term;
898 <        idat.vpair += vterm;
899 <        epot += idat.sw * vterm;
898 >        *(idat.vpair) += vterm;
899 >        epot +=  *(idat.sw)  * vterm;
900  
901          // calculate the derivatives for the forces and torques
902  
# Line 817 | Line 910 | namespace OpenMD {
910        }
911      }
912  
913 <    idat.pot += epot;
914 <    idat.f1 += dVdr;
913 >    (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
914 >    *(idat.f1) += dVdr;
915  
916      if (i_is_Dipole || i_is_Quadrupole)
917 <      idat.t1 -= cross(uz_i, duduz_i);
917 >      *(idat.t1) -= cross(uz_i, duduz_i);
918      if (i_is_Quadrupole) {
919 <      idat.t1 -= cross(ux_i, dudux_i);
920 <      idat.t1 -= cross(uy_i, duduy_i);
919 >      *(idat.t1) -= cross(ux_i, dudux_i);
920 >      *(idat.t1) -= cross(uy_i, duduy_i);
921      }
922 <
922 >    
923      if (j_is_Dipole || j_is_Quadrupole)
924 <      idat.t2 -= cross(uz_j, duduz_j);
924 >      *(idat.t2) -= cross(uz_j, duduz_j);
925      if (j_is_Quadrupole) {
926 <      idat.t2 -= cross(uz_j, dudux_j);
927 <      idat.t2 -= cross(uz_j, duduy_j);
926 >      *(idat.t2) -= cross(uz_j, dudux_j);
927 >      *(idat.t2) -= cross(uz_j, duduy_j);
928      }
929  
930      return;
931    }  
932  
933 <  void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
933 >  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
934  
935      if (!initialized_) initialize();
936      
937 <    ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
938 <    ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
937 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
938 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
939      
940      // logicals
941  
# Line 859 | Line 952 | namespace OpenMD {
952  
953      if (i_is_Charge) {
954        q_i = data1.charge;
955 <      skdat.skippedCharge2 += q_i;
955 >      *(idat.skippedCharge2) += q_i;
956      }
957  
958      if (j_is_Charge) {
959        q_j = data2.charge;
960 <      skdat.skippedCharge1 += q_j;
960 >      *(idat.skippedCharge1) += q_j;
961      }
962  
963      // the rest of this function should only be necessary for reaction field.
964  
965 <    if (summationMethod_ == REACTION_FIELD) {
965 >    if (summationMethod_ == esm_REACTION_FIELD) {
966        RealType riji, ri2, ri3;
967 <      RealType q_i, mu_i, ct_i;
968 <      RealType q_j, mu_j, ct_j;
969 <      RealType preVal, rfVal, vterm, dudr, pref, myPot;
967 >      RealType mu_i, ct_i;
968 >      RealType mu_j, ct_j;
969 >      RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
970        Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
971  
972        // some variables we'll need independent of electrostatic type:
973        
974 <      riji = 1.0 / skdat.rij;
975 <      rhat = skdat.d  * riji;
974 >      riji = 1.0 /  *(idat.rij) ;
975 >      rhat =  *(idat.d)  * riji;
976  
977        if (i_is_Dipole) {
978          mu_i = data1.dipole_moment;
979 <        uz_i = skdat.eFrame1.getColumn(2);      
979 >        uz_i = idat.eFrame1->getColumn(2);      
980          ct_i = dot(uz_i, rhat);
981          duduz_i = V3Zero;
982        }
983              
984        if (j_is_Dipole) {
985          mu_j = data2.dipole_moment;
986 <        uz_j = skdat.eFrame2.getColumn(2);      
986 >        uz_j = idat.eFrame2->getColumn(2);      
987          ct_j = dot(uz_j, rhat);
988          duduz_j = V3Zero;
989        }
990      
991        if (i_is_Charge) {
992          if (j_is_Charge) {
993 <          preVal = skdat.electroMult * pre11_ * q_i * q_j;
994 <          rfVal = preRF_ * skdat.rij * skdat.rij;
993 >          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
994 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
995            vterm = preVal * rfVal;
996 <          myPot += skdat.sw * vterm;        
997 <          dudr  = skdat.sw * preVal * 2.0 * rfVal * riji;        
996 >          myPot +=  *(idat.sw)  * vterm;        
997 >          dudr  =  *(idat.sw)  * preVal * 2.0 * rfVal * riji;        
998            dVdr += dudr * rhat;
999          }
1000          
1001          if (j_is_Dipole) {
1002            ri2 = riji * riji;
1003            ri3 = ri2 * riji;        
1004 <          pref = skdat.electroMult * pre12_ * q_i * mu_j;
1005 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1006 <          myPot += skdat.sw * vterm;        
1007 <          dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1008 <          duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1004 >          pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
1005 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
1006 >          myPot +=  *(idat.sw)  * vterm;        
1007 >          dVdr += - *(idat.sw)  * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1008 >          duduz_j += - *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1009          }
1010        }
1011        if (i_is_Dipole) {
1012          if (j_is_Charge) {
1013            ri2 = riji * riji;
1014            ri3 = ri2 * riji;        
1015 <          pref = skdat.electroMult * pre12_ * q_j * mu_i;
1016 <          vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1017 <          myPot += skdat.sw * vterm;        
1018 <          dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1019 <          duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1015 >          pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1016 >          vterm = - pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1017 >          myPot +=  *(idat.sw)  * vterm;        
1018 >          dVdr +=  *(idat.sw)  * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1019 >          duduz_i +=  *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij));
1020          }
1021        }
1022        
1023        // accumulate the forces and torques resulting from the self term
1024 <      skdat.pot += myPot;
1025 <      skdat.f1 += dVdr;
1024 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1025 >      *(idat.f1) += dVdr;
1026        
1027        if (i_is_Dipole)
1028 <        skdat.t1 -= cross(uz_i, duduz_i);
1028 >        *(idat.t1) -= cross(uz_i, duduz_i);
1029        if (j_is_Dipole)
1030 <        skdat.t2 -= cross(uz_j, duduz_j);
1030 >        *(idat.t2) -= cross(uz_j, duduz_j);
1031      }
1032    }
1033      
1034 <  void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1034 >  void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1035      RealType mu1, preVal, chg1, self;
1036      
1037      if (!initialized_) initialize();
1038 <    
1039 <    ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1038 >
1039 >    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1040    
1041      // logicals
949 –
1042      bool i_is_Charge = data.is_Charge;
1043      bool i_is_Dipole = data.is_Dipole;
1044  
1045 <    if (summationMethod_ == REACTION_FIELD) {
1045 >    if (summationMethod_ == esm_REACTION_FIELD) {
1046        if (i_is_Dipole) {
1047          mu1 = data.dipole_moment;          
1048          preVal = pre22_ * preRF2_ * mu1 * mu1;
1049 <        scdat.pot -= 0.5 * preVal;
1049 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1050          
1051          // The self-correction term adds into the reaction field vector
1052 <        Vector3d uz_i = scdat.eFrame.getColumn(2);
1052 >        Vector3d uz_i = sdat.eFrame->getColumn(2);
1053          Vector3d ei = preVal * uz_i;
1054  
1055          // This looks very wrong.  A vector crossed with itself is zero.
1056 <        scdat.t -= cross(uz_i, ei);
1056 >        *(sdat.t) -= cross(uz_i, ei);
1057        }
1058 <    } else if (summationMethod_ == SHIFTED_FORCE || summationMethod_ == SHIFTED_POTENTIAL) {
1058 >    } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1059        if (i_is_Charge) {        
1060          chg1 = data.charge;
1061          if (screeningMethod_ == DAMPED) {
1062 <          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1062 >          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1063          } else {        
1064 <          self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1064 >          self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1065          }
1066 <        scdat.pot += self;
1066 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1067        }
1068      }
1069    }
1070 +
1071 +  RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1072 +    // This seems to work moderately well as a default.  There's no
1073 +    // inherent scale for 1/r interactions that we can standardize.
1074 +    // 12 angstroms seems to be a reasonably good guess for most
1075 +    // cases.
1076 +    return 12.0;
1077 +  }
1078   }

Diff Legend

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