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