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 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC vs.
Revision 1668 by gezelter, Fri Jan 6 19:03:05 2012 UTC

# Line 34 | Line 34
34   * work.  Good starting points are:
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 < * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
37 > * [2]  Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <stdio.h>
# Line 47 | Line 48
48   #include "utils/simError.h"
49   #include "types/NonBondedInteractionType.hpp"
50   #include "types/DirectionalAtomType.hpp"
51 + #include "io/Globals.hpp"
52  
51 –
53   namespace OpenMD {
54    
55    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
56 <                                  forceField_(NULL) {}
56 >                                  forceField_(NULL), info_(NULL),
57 >                                  haveCutoffRadius_(false),
58 >                                  haveDampingAlpha_(false),
59 >                                  haveDielectric_(false),
60 >                                  haveElectroSpline_(false)
61 >  {}
62    
63    void Electrostatic::initialize() {
64 +    
65 +    Globals* simParams_ = info_->getSimParams();
66 +
67 +    summationMap_["HARD"]               = esm_HARD;
68 +    summationMap_["NONE"]               = esm_HARD;
69 +    summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
70 +    summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
71 +    summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
72 +    summationMap_["REACTION_FIELD"]     = esm_REACTION_FIELD;    
73 +    summationMap_["EWALD_FULL"]         = esm_EWALD_FULL;        
74 +    summationMap_["EWALD_PME"]          = esm_EWALD_PME;        
75 +    summationMap_["EWALD_SPME"]         = esm_EWALD_SPME;        
76 +    screeningMap_["DAMPED"]             = DAMPED;
77 +    screeningMap_["UNDAMPED"]           = UNDAMPED;
78 +
79      // these prefactors convert the multipole interactions into kcal / mol
80      // all were computed assuming distances are measured in angstroms
81      // Charge-Charge, assuming charges are measured in electrons
# Line 79 | Line 100 | namespace OpenMD {
100      
101      // variables to handle different summation methods for long-range
102      // electrostatics:
103 <    summationMethod_ = NONE;    
103 >    summationMethod_ = esm_HARD;    
104      screeningMethod_ = UNDAMPED;
105      dielectric_ = 1.0;
106      one_third_ = 1.0 / 3.0;
86 –    haveDefaultCutoff_ = false;
87 –    haveDampingAlpha_ = false;
88 –    haveDielectric_ = false;  
89 –    haveElectroSpline_ = false;
107    
108 +    // check the summation method:
109 +    if (simParams_->haveElectrostaticSummationMethod()) {
110 +      string myMethod = simParams_->getElectrostaticSummationMethod();
111 +      toUpper(myMethod);
112 +      map<string, ElectrostaticSummationMethod>::iterator i;
113 +      i = summationMap_.find(myMethod);
114 +      if ( i != summationMap_.end() ) {
115 +        summationMethod_ = (*i).second;
116 +      } else {
117 +        // throw error
118 +        sprintf( painCave.errMsg,
119 +                 "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
120 +                 "\t(Input file specified %s .)\n"
121 +                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
122 +                 "\t\"shifted_potential\", \"shifted_force\", or \n"
123 +                 "\t\"reaction_field\".\n", myMethod.c_str() );
124 +        painCave.isFatal = 1;
125 +        simError();
126 +      }
127 +    } else {
128 +      // set ElectrostaticSummationMethod to the cutoffMethod:
129 +      if (simParams_->haveCutoffMethod()){
130 +        string myMethod = simParams_->getCutoffMethod();
131 +        toUpper(myMethod);
132 +        map<string, ElectrostaticSummationMethod>::iterator i;
133 +        i = summationMap_.find(myMethod);
134 +        if ( i != summationMap_.end() ) {
135 +          summationMethod_ = (*i).second;
136 +        }
137 +      }
138 +    }
139 +    
140 +    if (summationMethod_ == esm_REACTION_FIELD) {        
141 +      if (!simParams_->haveDielectric()) {
142 +        // throw warning
143 +        sprintf( painCave.errMsg,
144 +                 "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
145 +                 "\tA default value of %f will be used for the dielectric.\n", dielectric_);
146 +        painCave.isFatal = 0;
147 +        painCave.severity = OPENMD_INFO;
148 +        simError();
149 +      } else {
150 +        dielectric_ = simParams_->getDielectric();      
151 +      }
152 +      haveDielectric_ = true;
153 +    }
154 +    
155 +    if (simParams_->haveElectrostaticScreeningMethod()) {
156 +      string myScreen = simParams_->getElectrostaticScreeningMethod();
157 +      toUpper(myScreen);
158 +      map<string, ElectrostaticScreeningMethod>::iterator i;
159 +      i = screeningMap_.find(myScreen);
160 +      if ( i != screeningMap_.end()) {
161 +        screeningMethod_ = (*i).second;
162 +      } else {
163 +        sprintf( painCave.errMsg,
164 +                 "SimInfo error: Unknown electrostaticScreeningMethod.\n"
165 +                 "\t(Input file specified %s .)\n"
166 +                 "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
167 +                 "or \"damped\".\n", myScreen.c_str() );
168 +        painCave.isFatal = 1;
169 +        simError();
170 +      }
171 +    }
172 +
173 +    // check to make sure a cutoff value has been set:
174 +    if (!haveCutoffRadius_) {
175 +      sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
176 +               "Cutoff value!\n");
177 +      painCave.severity = OPENMD_ERROR;
178 +      painCave.isFatal = 1;
179 +      simError();
180 +    }
181 +          
182 +    if (screeningMethod_ == DAMPED) {      
183 +      if (!simParams_->haveDampingAlpha()) {
184 +        // first set a cutoff dependent alpha value
185 +        // we assume alpha depends linearly with rcut from 0 to 20.5 ang
186 +        dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
187 +        if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
188 +        
189 +        // throw warning
190 +        sprintf( painCave.errMsg,
191 +                 "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
192 +                 "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
193 +                 dampingAlpha_, cutoffRadius_);
194 +        painCave.severity = OPENMD_INFO;
195 +        painCave.isFatal = 0;
196 +        simError();
197 +      } else {
198 +        dampingAlpha_ = simParams_->getDampingAlpha();
199 +      }
200 +      haveDampingAlpha_ = true;
201 +    }
202 +
203      // find all of the Electrostatic atom Types:
204      ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
205      ForceField::AtomTypeContainer::MapTypeIterator i;
206      AtomType* at;
207 <
207 >    
208      for (at = atomTypes->beginType(i); at != NULL;
209           at = atomTypes->nextType(i)) {
210        
# Line 100 | Line 212 | namespace OpenMD {
212          addType(at);
213      }
214      
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 –    }
215  
216 <    defaultCutoff2_ = defaultCutoff_ * defaultCutoff_;
217 <    rcuti_ = 1.0 / defaultCutoff_;
216 >    cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
217 >    rcuti_ = 1.0 / cutoffRadius_;
218      rcuti2_ = rcuti_ * rcuti_;
219      rcuti3_ = rcuti2_ * rcuti_;
220      rcuti4_ = rcuti2_ * rcuti2_;
221  
222      if (screeningMethod_ == DAMPED) {
223 <      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 <
223 >      
224        alpha2_ = dampingAlpha_ * dampingAlpha_;
225        alpha4_ = alpha2_ * alpha2_;
226        alpha6_ = alpha4_ * alpha2_;
227        alpha8_ = alpha4_ * alpha4_;
228        
229 <      constEXP_ = exp(-alpha2_ * defaultCutoff2_);
229 >      constEXP_ = exp(-alpha2_ * cutoffRadius2_);
230        invRootPi_ = 0.56418958354775628695;
231        alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
232  
233 <      c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * rcuti_;
233 >      c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
234        c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
235        c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
236        c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_;
# Line 148 | Line 245 | namespace OpenMD {
245        c6c_ = 9.0 * c5c_ * rcuti2_;
246      }
247    
248 <    if (summationMethod_ == REACTION_FIELD) {
249 <      if (haveDielectric_) {
250 <        preRF_ = (dielectric_ - 1.0) /
251 <            ((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 <      }
248 >    if (summationMethod_ == esm_REACTION_FIELD) {
249 >      preRF_ = (dielectric_ - 1.0) /
250 >        ((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
251 >      preRF2_ = 2.0 * preRF_;
252      }
253 <                              
254 <    RealType dx = defaultCutoff_ / RealType(np_ - 1);
253 >    
254 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
255 >    // have charged atoms longer than the cutoffRadius away from each
256 >    // other.  Splining may not be the best choice here.  Direct calls
257 >    // to erfc might be preferrable.
258 >
259 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
260      RealType rval;
261      vector<RealType> rvals;
262      vector<RealType> yvals;
# Line 321 | Line 415 | namespace OpenMD {
415      return;
416    }
417    
418 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
419 <                                                    RealType theRSW ) {
420 <    defaultCutoff_ = theECR;
421 <    rrf_ = defaultCutoff_;
328 <    rt_ = theRSW;
329 <    haveDefaultCutoff_ = true;
418 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
419 >    cutoffRadius_ = rCut;
420 >    rrf_ = cutoffRadius_;
421 >    haveCutoffRadius_ = true;
422    }
423 +
424 +  void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
425 +    rt_ = rSwitch;
426 +  }
427    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
428      summationMethod_ = esm;
429    }
# Line 343 | Line 439 | namespace OpenMD {
439      haveDielectric_ = true;
440    }
441  
442 <  void Electrostatic::calcForce(InteractionData idat) {
442 >  void Electrostatic::calcForce(InteractionData &idat) {
443  
444      // utility variables.  Should clean these up and use the Vector3d and
445      // Mat3x3d to replace as many as we can in future versions:
# Line 357 | Line 453 | namespace OpenMD {
453      RealType ct_i, ct_j, ct_ij, a1;
454      RealType riji, ri, ri2, ri3, ri4;
455      RealType pref, vterm, epot, dudr;
456 +    RealType vpair(0.0);
457      RealType scale, sc2;
458      RealType pot_term, preVal, rfVal;
459      RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
460      RealType preSw, preSwSc;
461      RealType c1, c2, c3, c4;
462 <    RealType erfcVal, derfcVal;
462 >    RealType erfcVal(1.0), derfcVal(0.0);
463      RealType BigR;
464 +    RealType two(2.0), three(3.0);
465  
466      Vector3d Q_i, Q_j;
467      Vector3d ux_i, uy_i, uz_i;
# Line 373 | Line 471 | namespace OpenMD {
471      Vector3d rhatdot2, rhatc4;
472      Vector3d dVdr;
473  
474 +    // variables for indirect (reaction field) interactions for excluded pairs:
475 +    RealType indirect_Pot(0.0);
476 +    RealType indirect_vpair(0.0);
477 +    Vector3d indirect_dVdr(V3Zero);
478 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
479 +
480      pair<RealType, RealType> res;
481      
482      if (!initialized_) initialize();
483      
484 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
485 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
484 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
485 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
486      
487      // some variables we'll need independent of electrostatic type:
488  
489 <    riji = 1.0 / idat.rij;
490 <    Vector3d rhat = idat.d  * riji;
489 >    riji = 1.0 /  *(idat.rij) ;
490 >    Vector3d rhat =  *(idat.d)   * riji;
491  
492      // logicals
493  
# Line 397 | Line 501 | namespace OpenMD {
501      bool j_is_SplitDipole = data2.is_SplitDipole;
502      bool j_is_Quadrupole = data2.is_Quadrupole;
503      
504 <    if (i_is_Charge)
504 >    if (i_is_Charge) {
505        q_i = data1.charge;
506 +      if (idat.excluded) {
507 +        *(idat.skippedCharge2) += q_i;
508 +      }
509 +    }
510  
511      if (i_is_Dipole) {
512        mu_i = data1.dipole_moment;
513 <      uz_i = idat.eFrame1.getColumn(2);
513 >      uz_i = idat.eFrame1->getColumn(2);
514        
515        ct_i = dot(uz_i, rhat);
516  
# Line 418 | Line 526 | namespace OpenMD {
526        qyy_i = Q_i.y();
527        qzz_i = Q_i.z();
528        
529 <      ux_i = idat.eFrame1.getColumn(0);
530 <      uy_i = idat.eFrame1.getColumn(1);
531 <      uz_i = idat.eFrame1.getColumn(2);
529 >      ux_i = idat.eFrame1->getColumn(0);
530 >      uy_i = idat.eFrame1->getColumn(1);
531 >      uz_i = idat.eFrame1->getColumn(2);
532  
533        cx_i = dot(ux_i, rhat);
534        cy_i = dot(uy_i, rhat);
# Line 431 | Line 539 | namespace OpenMD {
539        duduz_i = V3Zero;
540      }
541  
542 <    if (j_is_Charge)
542 >    if (j_is_Charge) {
543        q_j = data2.charge;
544 +      if (idat.excluded) {
545 +        *(idat.skippedCharge1) += q_j;
546 +      }
547 +    }
548  
549 +
550      if (j_is_Dipole) {
551        mu_j = data2.dipole_moment;
552 <      uz_j = idat.eFrame2.getColumn(2);
552 >      uz_j = idat.eFrame2->getColumn(2);
553        
554        ct_j = dot(uz_j, rhat);
555  
# Line 452 | Line 565 | namespace OpenMD {
565        qyy_j = Q_j.y();
566        qzz_j = Q_j.z();
567        
568 <      ux_j = idat.eFrame2.getColumn(0);
569 <      uy_j = idat.eFrame2.getColumn(1);
570 <      uz_j = idat.eFrame2.getColumn(2);
568 >      ux_j = idat.eFrame2->getColumn(0);
569 >      uy_j = idat.eFrame2->getColumn(1);
570 >      uz_j = idat.eFrame2->getColumn(2);
571  
572        cx_j = dot(ux_j, rhat);
573        cy_j = dot(uy_j, rhat);
# Line 473 | Line 586 | namespace OpenMD {
586        if (j_is_Charge) {
587          if (screeningMethod_ == DAMPED) {
588            // assemble the damping variables
589 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
590 <          erfcVal = res.first;
591 <          derfcVal = res.second;
589 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
590 >          //erfcVal = res.first;
591 >          //derfcVal = res.second;
592 >
593 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
594 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
595 >
596            c1 = erfcVal * riji;
597            c2 = (-derfcVal + c1) * riji;
598          } else {
# Line 483 | Line 600 | namespace OpenMD {
600            c2 = c1 * riji;
601          }
602  
603 <        preVal = idat.electroMult * pre11_ * q_i * q_j;
603 >        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
604          
605 <        if (summationMethod_ == SHIFTED_POTENTIAL) {
605 >        if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
606            vterm = preVal * (c1 - c1c_);
607 <          dudr  = -idat.sw * preVal * c2;
607 >          dudr  = - *(idat.sw)  * preVal * c2;
608  
609 <        } else if (summationMethod_ == SHIFTED_FORCE)  {
610 <          vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) );
611 <          dudr  = idat.sw * preVal * (c2c_ - c2);
609 >        } else if (summationMethod_ == esm_SHIFTED_FORCE)  {
610 >          vterm = preVal * ( c1 - c1c_ + c2c_*( *(idat.rij)  - cutoffRadius_) );
611 >          dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
612  
613 <        } else if (summationMethod_ == REACTION_FIELD) {
614 <          rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
498 <          vterm = preVal * ( riji + rfVal );            
499 <          dudr  = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
613 >        } else if (summationMethod_ == esm_REACTION_FIELD) {
614 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
615  
616 +          vterm = preVal * ( riji + rfVal );            
617 +          dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
618 +          
619 +          // if this is an excluded pair, there are still indirect
620 +          // interactions via the reaction field we must worry about:
621 +
622 +          if (idat.excluded) {
623 +            indirect_vpair += preVal * rfVal;
624 +            indirect_Pot += *(idat.sw) * preVal * rfVal;
625 +            indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
626 +          }
627 +          
628          } else {
502 –          vterm = preVal * riji * erfcVal;            
629  
630 <          dudr  = - idat.sw * preVal * c2;
630 >          vterm = preVal * riji * erfcVal;          
631 >          dudr  = -  *(idat.sw)  * preVal * c2;
632  
633          }
507 –
508 –        idat.vpair += vterm;
509 –        epot += idat.sw * vterm;
634  
635 <        dVdr += dudr * rhat;      
635 >        vpair += vterm;
636 >        epot +=  *(idat.sw)  * vterm;
637 >        dVdr += dudr * rhat;                
638        }
639  
640        if (j_is_Dipole) {
641          // pref is used by all the possible methods
642 <        pref = idat.electroMult * pre12_ * q_i * mu_j;
643 <        preSw = idat.sw * pref;
642 >        pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
643 >        preSw =  *(idat.sw)  * pref;
644  
645 <        if (summationMethod_ == REACTION_FIELD) {
645 >        if (summationMethod_ == esm_REACTION_FIELD) {
646            ri2 = riji * riji;
647            ri3 = ri2 * riji;
648      
649 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
650 <          idat.vpair += vterm;
651 <          epot += idat.sw * vterm;
649 >          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
650 >          vpair += vterm;
651 >          epot +=  *(idat.sw)  * vterm;
652  
653 <          dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
654 <          duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);  
653 >          dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
654 >          duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
655  
656 +          // Even if we excluded this pair from direct interactions,
657 +          // we still have the reaction-field-mediated charge-dipole
658 +          // interaction:
659 +
660 +          if (idat.excluded) {
661 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
662 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
663 +            indirect_dVdr += preSw * preRF2_ * uz_j;
664 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
665 +          }
666 +                      
667          } else {
668            // determine the inverse r used if we have split dipoles
669            if (j_is_SplitDipole) {
670 <            BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
670 >            BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
671              ri = 1.0 / BigR;
672 <            scale = idat.rij * ri;
672 >            scale =  *(idat.rij)  * ri;
673            } else {
674              ri = riji;
675              scale = 1.0;
# Line 542 | Line 679 | namespace OpenMD {
679  
680            if (screeningMethod_ == DAMPED) {
681              // assemble the damping variables
682 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
683 <            erfcVal = res.first;
684 <            derfcVal = res.second;
682 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
683 >            //erfcVal = res.first;
684 >            //derfcVal = res.second;
685 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
686 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
687              c1 = erfcVal * ri;
688              c2 = (-derfcVal + c1) * ri;
689              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 559 | Line 698 | namespace OpenMD {
698            // calculate the potential
699            pot_term =  scale * c2;
700            vterm = -pref * ct_j * pot_term;
701 <          idat.vpair += vterm;
702 <          epot += idat.sw * vterm;
701 >          vpair += vterm;
702 >          epot +=  *(idat.sw)  * vterm;
703              
704            // calculate derivatives for forces and torques
705  
# Line 575 | Line 714 | namespace OpenMD {
714          cx2 = cx_j * cx_j;
715          cy2 = cy_j * cy_j;
716          cz2 = cz_j * cz_j;
717 <        pref =  idat.electroMult * pre14_ * q_i * one_third_;
717 >        pref =   *(idat.electroMult) * pre14_ * q_i * one_third_;
718            
719          if (screeningMethod_ == DAMPED) {
720            // assemble the damping variables
721 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
722 <          erfcVal = res.first;
723 <          derfcVal = res.second;
721 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
722 >          //erfcVal = res.first;
723 >          //derfcVal = res.second;
724 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
725 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
726            c1 = erfcVal * riji;
727            c2 = (-derfcVal + c1) * riji;
728            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 594 | Line 735 | namespace OpenMD {
735          }
736  
737          // precompute variables for convenience
738 <        preSw = idat.sw * pref;
738 >        preSw =  *(idat.sw)  * pref;
739          c2ri = c2 * riji;
740          c3ri = c3 * riji;
741 <        c4rij = c4 * idat.rij;
742 <        rhatdot2 = 2.0 * rhat * c3;
741 >        c4rij = c4 *  *(idat.rij) ;
742 >        rhatdot2 = two * rhat * c3;
743          rhatc4 = rhat * c4rij;
744  
745          // calculate the potential
# Line 606 | Line 747 | namespace OpenMD {
747                       qyy_j * (cy2*c3 - c2ri) +
748                       qzz_j * (cz2*c3 - c2ri) );
749          vterm = pref * pot_term;
750 <        idat.vpair += vterm;
751 <        epot += idat.sw * vterm;
750 >        vpair += vterm;
751 >        epot +=  *(idat.sw)  * vterm;
752                  
753          // calculate derivatives for the forces and torques
754  
755 <        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
756 <                           qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
757 <                           qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
755 >        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
756 >                           qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
757 >                           qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
758                            
759          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
760          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
# Line 625 | Line 766 | namespace OpenMD {
766  
767        if (j_is_Charge) {
768          // variables used by all the methods
769 <        pref = idat.electroMult * pre12_ * q_j * mu_i;
770 <        preSw = idat.sw * pref;
769 >        pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
770 >        preSw =  *(idat.sw)  * pref;
771  
772 <        if (summationMethod_ == REACTION_FIELD) {
772 >        if (summationMethod_ == esm_REACTION_FIELD) {
773  
774            ri2 = riji * riji;
775            ri3 = ri2 * riji;
776  
777 <          vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
778 <          idat.vpair += vterm;
779 <          epot += idat.sw * vterm;
777 >          vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
778 >          vpair += vterm;
779 >          epot +=  *(idat.sw)  * vterm;
780            
781 <          dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
781 >          dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
782            
783 <          duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
783 >          duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
784 >
785 >          // Even if we excluded this pair from direct interactions,
786 >          // we still have the reaction-field-mediated charge-dipole
787 >          // interaction:
788 >
789 >          if (idat.excluded) {
790 >            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
791 >            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
792 >            indirect_dVdr += -preSw * preRF2_ * uz_i;
793 >            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
794 >          }
795              
796          } else {
797            
798            // determine inverse r if we are using split dipoles
799            if (i_is_SplitDipole) {
800 <            BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
800 >            BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
801              ri = 1.0 / BigR;
802 <            scale = idat.rij * ri;
802 >            scale =  *(idat.rij)  * ri;
803            } else {
804              ri = riji;
805              scale = 1.0;
# Line 657 | Line 809 | namespace OpenMD {
809              
810            if (screeningMethod_ == DAMPED) {
811              // assemble the damping variables
812 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
813 <            erfcVal = res.first;
814 <            derfcVal = res.second;
812 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
813 >            //erfcVal = res.first;
814 >            //derfcVal = res.second;
815 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
816 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
817              c1 = erfcVal * ri;
818              c2 = (-derfcVal + c1) * ri;
819              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 674 | Line 828 | namespace OpenMD {
828            // calculate the potential
829            pot_term = c2 * scale;
830            vterm = pref * ct_i * pot_term;
831 <          idat.vpair += vterm;
832 <          epot += idat.sw * vterm;
831 >          vpair += vterm;
832 >          epot +=  *(idat.sw)  * vterm;
833  
834            // calculate derivatives for the forces and torques
835            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
# Line 687 | Line 841 | namespace OpenMD {
841          // variables used by all methods
842          ct_ij = dot(uz_i, uz_j);
843  
844 <        pref = idat.electroMult * pre22_ * mu_i * mu_j;
845 <        preSw = idat.sw * pref;
844 >        pref =  *(idat.electroMult) * pre22_ * mu_i * mu_j;
845 >        preSw =  *(idat.sw)  * pref;
846  
847 <        if (summationMethod_ == REACTION_FIELD) {
847 >        if (summationMethod_ == esm_REACTION_FIELD) {
848            ri2 = riji * riji;
849            ri3 = ri2 * riji;
850            ri4 = ri2 * ri2;
851  
852            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
853                             preRF2_ * ct_ij );
854 <          idat.vpair += vterm;
855 <          epot += idat.sw * vterm;
854 >          vpair += vterm;
855 >          epot +=  *(idat.sw)  * vterm;
856              
857            a1 = 5.0 * ct_i * ct_j - ct_ij;
858              
859 <          dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
859 >          dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
860  
861 <          duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
862 <          duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
861 >          duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
862 >          duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
863  
864 +          if (idat.excluded) {
865 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
866 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
867 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
868 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
869 +          }
870 +
871          } else {
872            
873            if (i_is_SplitDipole) {
874              if (j_is_SplitDipole) {
875 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
875 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
876              } else {
877 <              BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
877 >              BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i);
878              }
879              ri = 1.0 / BigR;
880 <            scale = idat.rij * ri;
880 >            scale =  *(idat.rij)  * ri;
881            } else {
882              if (j_is_SplitDipole) {
883 <              BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
883 >              BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j);
884                ri = 1.0 / BigR;
885 <              scale = idat.rij * ri;
885 >              scale =  *(idat.rij)  * ri;
886              } else {
887                ri = riji;
888                scale = 1.0;
# Line 729 | Line 890 | namespace OpenMD {
890            }
891            if (screeningMethod_ == DAMPED) {
892              // assemble damping variables
893 <            res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
894 <            erfcVal = res.first;
895 <            derfcVal = res.second;
893 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
894 >            //erfcVal = res.first;
895 >            //derfcVal = res.second;
896 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
897 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
898              c1 = erfcVal * ri;
899              c2 = (-derfcVal + c1) * ri;
900              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 751 | Line 914 | namespace OpenMD {
914            preSwSc = preSw * scale;
915            c2ri = c2 * ri;
916            c3ri = c3 * ri;
917 <          c4rij = c4 * idat.rij;
917 >          c4rij = c4 *  *(idat.rij) ;
918  
919            // calculate the potential
920            pot_term = (ct_ij * c2ri - ctidotj * c3);
921            vterm = pref * pot_term;
922 <          idat.vpair += vterm;
923 <          epot += idat.sw * vterm;
922 >          vpair += vterm;
923 >          epot +=  *(idat.sw)  * vterm;
924  
925            // calculate derivatives for the forces and torques
926            dVdr += preSwSc * ( ctidotj * rhat * c4rij  -
# Line 776 | Line 939 | namespace OpenMD {
939          cy2 = cy_i * cy_i;
940          cz2 = cz_i * cz_i;
941  
942 <        pref = idat.electroMult * pre14_ * q_j * one_third_;
942 >        pref =  *(idat.electroMult) * pre14_ * q_j * one_third_;
943  
944          if (screeningMethod_ == DAMPED) {
945            // assemble the damping variables
946 <          res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
947 <          erfcVal = res.first;
948 <          derfcVal = res.second;
946 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
947 >          //erfcVal = res.first;
948 >          //derfcVal = res.second;
949 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
950 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
951            c1 = erfcVal * riji;
952            c2 = (-derfcVal + c1) * riji;
953            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 795 | Line 960 | namespace OpenMD {
960          }
961            
962          // precompute some variables for convenience
963 <        preSw = idat.sw * pref;
963 >        preSw =  *(idat.sw)  * pref;
964          c2ri = c2 * riji;
965          c3ri = c3 * riji;
966 <        c4rij = c4 * idat.rij;
967 <        rhatdot2 = 2.0 * rhat * c3;
966 >        c4rij = c4 *  *(idat.rij) ;
967 >        rhatdot2 = two * rhat * c3;
968          rhatc4 = rhat * c4rij;
969  
970          // calculate the potential
# Line 808 | Line 973 | namespace OpenMD {
973                       qzz_i * (cz2 * c3 - c2ri) );
974          
975          vterm = pref * pot_term;
976 <        idat.vpair += vterm;
977 <        epot += idat.sw * vterm;
976 >        vpair += vterm;
977 >        epot +=  *(idat.sw)  * vterm;
978  
979          // calculate the derivatives for the forces and torques
980  
981 <        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
982 <                          qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
983 <                          qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
981 >        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
982 >                          qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
983 >                          qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
984  
985          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
986          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
# Line 823 | Line 988 | namespace OpenMD {
988        }
989      }
990  
826 –    idat.pot += epot;
827 –    idat.f1 += dVdr;
991  
992 <    if (i_is_Dipole || i_is_Quadrupole)
993 <      idat.t1 -= cross(uz_i, duduz_i);
994 <    if (i_is_Quadrupole) {
995 <      idat.t1 -= cross(ux_i, dudux_i);
996 <      idat.t1 -= cross(uy_i, duduy_i);
997 <    }
992 >    if (!idat.excluded) {
993 >      *(idat.vpair) += vpair;
994 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
995 >      *(idat.f1) += dVdr;
996 >      
997 >      if (i_is_Dipole || i_is_Quadrupole)
998 >        *(idat.t1) -= cross(uz_i, duduz_i);
999 >      if (i_is_Quadrupole) {
1000 >        *(idat.t1) -= cross(ux_i, dudux_i);
1001 >        *(idat.t1) -= cross(uy_i, duduy_i);
1002 >      }
1003 >      
1004 >      if (j_is_Dipole || j_is_Quadrupole)
1005 >        *(idat.t2) -= cross(uz_j, duduz_j);
1006 >      if (j_is_Quadrupole) {
1007 >        *(idat.t2) -= cross(uz_j, dudux_j);
1008 >        *(idat.t2) -= cross(uz_j, duduy_j);
1009 >      }
1010  
1011 <    if (j_is_Dipole || j_is_Quadrupole)
837 <      idat.t2 -= cross(uz_j, duduz_j);
838 <    if (j_is_Quadrupole) {
839 <      idat.t2 -= cross(uz_j, dudux_j);
840 <      idat.t2 -= cross(uz_j, duduy_j);
841 <    }
1011 >    } else {
1012  
1013 <    return;
1014 <  }  
1013 >      // only accumulate the forces and torques resulting from the
1014 >      // indirect reaction field terms.
1015  
1016 <  void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
1017 <
1018 <    if (!initialized_) initialize();
849 <    
850 <    ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
851 <    ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
852 <    
853 <    // logicals
854 <
855 <    bool i_is_Charge = data1.is_Charge;
856 <    bool i_is_Dipole = data1.is_Dipole;
857 <
858 <    bool j_is_Charge = data2.is_Charge;
859 <    bool j_is_Dipole = data2.is_Dipole;
860 <
861 <    RealType q_i, q_j;
862 <    
863 <    // The skippedCharge computation is needed by the real-space cutoff methods
864 <    // (i.e. shifted force and shifted potential)
865 <
866 <    if (i_is_Charge) {
867 <      q_i = data1.charge;
868 <      skdat.skippedCharge2 += q_i;
869 <    }
870 <
871 <    if (j_is_Charge) {
872 <      q_j = data2.charge;
873 <      skdat.skippedCharge1 += q_j;
874 <    }
875 <
876 <    // the rest of this function should only be necessary for reaction field.
877 <
878 <    if (summationMethod_ == REACTION_FIELD) {
879 <      RealType riji, ri2, ri3;
880 <      RealType q_i, mu_i, ct_i;
881 <      RealType q_j, mu_j, ct_j;
882 <      RealType preVal, rfVal, vterm, dudr, pref, myPot;
883 <      Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
884 <
885 <      // some variables we'll need independent of electrostatic type:
1016 >      *(idat.vpair) += indirect_vpair;
1017 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1018 >      *(idat.f1) += indirect_dVdr;
1019        
887 –      riji = 1.0 / skdat.rij;
888 –      rhat = skdat.d  * riji;
889 –
890 –      if (i_is_Dipole) {
891 –        mu_i = data1.dipole_moment;
892 –        uz_i = skdat.eFrame1.getColumn(2);      
893 –        ct_i = dot(uz_i, rhat);
894 –        duduz_i = V3Zero;
895 –      }
896 –            
897 –      if (j_is_Dipole) {
898 –        mu_j = data2.dipole_moment;
899 –        uz_j = skdat.eFrame2.getColumn(2);      
900 –        ct_j = dot(uz_j, rhat);
901 –        duduz_j = V3Zero;
902 –      }
903 –    
904 –      if (i_is_Charge) {
905 –        if (j_is_Charge) {
906 –          preVal = skdat.electroMult * pre11_ * q_i * q_j;
907 –          rfVal = preRF_ * skdat.rij * skdat.rij;
908 –          vterm = preVal * rfVal;
909 –          myPot += skdat.sw * vterm;        
910 –          dudr  = skdat.sw * preVal * 2.0 * rfVal * riji;        
911 –          dVdr += dudr * rhat;
912 –        }
913 –        
914 –        if (j_is_Dipole) {
915 –          ri2 = riji * riji;
916 –          ri3 = ri2 * riji;        
917 –          pref = skdat.electroMult * pre12_ * q_i * mu_j;
918 –          vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
919 –          myPot += skdat.sw * vterm;        
920 –          dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
921 –          duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
922 –        }
923 –      }
924 –      if (i_is_Dipole) {
925 –        if (j_is_Charge) {
926 –          ri2 = riji * riji;
927 –          ri3 = ri2 * riji;        
928 –          pref = skdat.electroMult * pre12_ * q_j * mu_i;
929 –          vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
930 –          myPot += skdat.sw * vterm;        
931 –          dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
932 –          duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
933 –        }
934 –      }
935 –      
936 –      // accumulate the forces and torques resulting from the self term
937 –      skdat.pot += myPot;
938 –      skdat.f1 += dVdr;
939 –      
1020        if (i_is_Dipole)
1021 <        skdat.t1 -= cross(uz_i, duduz_i);
1021 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
1022        if (j_is_Dipole)
1023 <        skdat.t2 -= cross(uz_j, duduz_j);
1023 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1024      }
1025 <  }
1025 >
1026 >
1027 >    return;
1028 >  }  
1029      
1030 <  void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1030 >  void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1031      RealType mu1, preVal, chg1, self;
1032      
1033      if (!initialized_) initialize();
1034 <    
1035 <    ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1034 >
1035 >    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1036    
1037      // logicals
955 –
1038      bool i_is_Charge = data.is_Charge;
1039      bool i_is_Dipole = data.is_Dipole;
1040  
1041 <    if (summationMethod_ == REACTION_FIELD) {
1041 >    if (summationMethod_ == esm_REACTION_FIELD) {
1042        if (i_is_Dipole) {
1043          mu1 = data.dipole_moment;          
1044          preVal = pre22_ * preRF2_ * mu1 * mu1;
1045 <        scdat.pot -= 0.5 * preVal;
1045 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1046          
1047          // The self-correction term adds into the reaction field vector
1048 <        Vector3d uz_i = scdat.eFrame.getColumn(2);
1048 >        Vector3d uz_i = sdat.eFrame->getColumn(2);
1049          Vector3d ei = preVal * uz_i;
1050  
1051          // This looks very wrong.  A vector crossed with itself is zero.
1052 <        scdat.t -= cross(uz_i, ei);
1052 >        *(sdat.t) -= cross(uz_i, ei);
1053        }
1054 <    } else if (summationMethod_ == SHIFTED_FORCE || summationMethod_ == SHIFTED_POTENTIAL) {
1054 >    } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1055        if (i_is_Charge) {        
1056          chg1 = data.charge;
1057          if (screeningMethod_ == DAMPED) {
1058 <          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1058 >          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1059          } else {        
1060 <          self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1060 >          self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1061          }
1062 <        scdat.pot += self;
1062 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1063        }
1064      }
1065    }
1066  
1067 <  RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1067 >  RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1068      // This seems to work moderately well as a default.  There's no
1069      // inherent scale for 1/r interactions that we can standardize.
1070      // 12 angstroms seems to be a reasonably good guess for most

Diff Legend

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