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 1601 by gezelter, Thu Aug 4 20:04:35 2011 UTC vs.
Revision 1723 by gezelter, Thu May 24 20:59:54 2012 UTC

# Line 36 | Line 36
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
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 46 | Line 47
47   #include "nonbonded/Electrostatic.hpp"
48   #include "utils/simError.h"
49   #include "types/NonBondedInteractionType.hpp"
50 < #include "types/DirectionalAtomType.hpp"
50 > #include "types/FixedChargeAdapter.hpp"
51 > #include "types/FluctuatingChargeAdapter.hpp"
52 > #include "types/MultipoleAdapter.hpp"
53   #include "io/Globals.hpp"
54 + #include "nonbonded/SlaterIntegrals.hpp"
55 + #include "utils/PhysicalConstants.hpp"
56  
57 +
58   namespace OpenMD {
59    
60    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
# Line 64 | Line 70 | namespace OpenMD {
70      Globals* simParams_ = info_->getSimParams();
71  
72      summationMap_["HARD"]               = esm_HARD;
73 +    summationMap_["NONE"]               = esm_HARD;
74      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
75      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
76      summationMap_["SHIFTED_FORCE"]      = esm_SHIFTED_FORCE;    
# Line 116 | Line 123 | namespace OpenMD {
123          sprintf( painCave.errMsg,
124                   "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
125                   "\t(Input file specified %s .)\n"
126 <                 "\telectrostaticSummationMethod must be one of: \"none\",\n"
126 >                 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
127                   "\t\"shifted_potential\", \"shifted_force\", or \n"
128                   "\t\"reaction_field\".\n", myMethod.c_str() );
129          painCave.isFatal = 1;
# Line 210 | Line 217 | namespace OpenMD {
217          addType(at);
218      }
219      
213 –
220      cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
221      rcuti_ = 1.0 / cutoffRadius_;
222      rcuti2_ = rcuti_ * rcuti_;
# Line 249 | Line 255 | namespace OpenMD {
255        preRF2_ = 2.0 * preRF_;
256      }
257      
258 <    RealType dx = cutoffRadius_ / RealType(np_ - 1);
258 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
259 >    // have charged atoms longer than the cutoffRadius away from each
260 >    // other.  Splining may not be the best choice here.  Direct calls
261 >    // to erfc might be preferrable.
262 >
263 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
264      RealType rval;
265      vector<RealType> rvals;
266      vector<RealType> yvals;
# Line 272 | Line 283 | namespace OpenMD {
283      electrostaticAtomData.is_Dipole = false;
284      electrostaticAtomData.is_SplitDipole = false;
285      electrostaticAtomData.is_Quadrupole = false;
286 +    electrostaticAtomData.is_Fluctuating = false;
287  
288 <    if (atomType->isCharge()) {
277 <      GenericData* data = atomType->getPropertyByName("Charge");
288 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
289  
290 <      if (data == NULL) {
280 <        sprintf( painCave.errMsg, "Electrostatic::addType could not find "
281 <                 "Charge\n"
282 <                 "\tparameters for atomType %s.\n",
283 <                 atomType->getName().c_str());
284 <        painCave.severity = OPENMD_ERROR;
285 <        painCave.isFatal = 1;
286 <        simError();                  
287 <      }
288 <      
289 <      DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
290 <      if (doubleData == NULL) {
291 <        sprintf( painCave.errMsg,
292 <                 "Electrostatic::addType could not convert GenericData to "
293 <                 "Charge for\n"
294 <                 "\tatom type %s\n", atomType->getName().c_str());
295 <        painCave.severity = OPENMD_ERROR;
296 <        painCave.isFatal = 1;
297 <        simError();          
298 <      }
290 >    if (fca.isFixedCharge()) {
291        electrostaticAtomData.is_Charge = true;
292 <      electrostaticAtomData.charge = doubleData->getData();          
292 >      electrostaticAtomData.fixedCharge = fca.getCharge();
293      }
294  
295 <    if (atomType->isDirectional()) {
296 <      DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
297 <      
306 <      if (daType->isDipole()) {
307 <        GenericData* data = daType->getPropertyByName("Dipole");
308 <        
309 <        if (data == NULL) {
310 <          sprintf( painCave.errMsg,
311 <                   "Electrostatic::addType could not find Dipole\n"
312 <                   "\tparameters for atomType %s.\n",
313 <                   daType->getName().c_str());
314 <          painCave.severity = OPENMD_ERROR;
315 <          painCave.isFatal = 1;
316 <          simError();                  
317 <        }
318 <      
319 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
320 <        if (doubleData == NULL) {
321 <          sprintf( painCave.errMsg,
322 <                   "Electrostatic::addType could not convert GenericData to "
323 <                   "Dipole Moment\n"
324 <                   "\tfor atom type %s\n", daType->getName().c_str());
325 <          painCave.severity = OPENMD_ERROR;
326 <          painCave.isFatal = 1;
327 <          simError();          
328 <        }
295 >    MultipoleAdapter ma = MultipoleAdapter(atomType);
296 >    if (ma.isMultipole()) {
297 >      if (ma.isDipole()) {
298          electrostaticAtomData.is_Dipole = true;
299 <        electrostaticAtomData.dipole_moment = doubleData->getData();
299 >        electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
300        }
301 <
333 <      if (daType->isSplitDipole()) {
334 <        GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
335 <        
336 <        if (data == NULL) {
337 <          sprintf(painCave.errMsg,
338 <                  "Electrostatic::addType could not find SplitDipoleDistance\n"
339 <                  "\tparameter for atomType %s.\n",
340 <                  daType->getName().c_str());
341 <          painCave.severity = OPENMD_ERROR;
342 <          painCave.isFatal = 1;
343 <          simError();                  
344 <        }
345 <      
346 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
347 <        if (doubleData == NULL) {
348 <          sprintf( painCave.errMsg,
349 <                   "Electrostatic::addType could not convert GenericData to "
350 <                   "SplitDipoleDistance for\n"
351 <                   "\tatom type %s\n", daType->getName().c_str());
352 <          painCave.severity = OPENMD_ERROR;
353 <          painCave.isFatal = 1;
354 <          simError();          
355 <        }
301 >      if (ma.isSplitDipole()) {
302          electrostaticAtomData.is_SplitDipole = true;
303 <        electrostaticAtomData.split_dipole_distance = doubleData->getData();
303 >        electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
304        }
305 <
360 <      if (daType->isQuadrupole()) {
361 <        GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
362 <        
363 <        if (data == NULL) {
364 <          sprintf( painCave.errMsg,
365 <                   "Electrostatic::addType could not find QuadrupoleMoments\n"
366 <                   "\tparameter for atomType %s.\n",
367 <                   daType->getName().c_str());
368 <          painCave.severity = OPENMD_ERROR;
369 <          painCave.isFatal = 1;
370 <          simError();                  
371 <        }
372 <        
305 >      if (ma.isQuadrupole()) {
306          // Quadrupoles in OpenMD are set as the diagonal elements
307          // of the diagonalized traceless quadrupole moment tensor.
308          // The column vectors of the unitary matrix that diagonalizes
309          // the quadrupole moment tensor become the eFrame (or the
310          // electrostatic version of the body-fixed frame.
378 –
379 –        Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
380 –        if (v3dData == NULL) {
381 –          sprintf( painCave.errMsg,
382 –                   "Electrostatic::addType could not convert GenericData to "
383 –                   "Quadrupole Moments for\n"
384 –                   "\tatom type %s\n", daType->getName().c_str());
385 –          painCave.severity = OPENMD_ERROR;
386 –          painCave.isFatal = 1;
387 –          simError();          
388 –        }
311          electrostaticAtomData.is_Quadrupole = true;
312 <        electrostaticAtomData.quadrupole_moments = v3dData->getData();
312 >        electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
313        }
314      }
315      
316 <    AtomTypeProperties atp = atomType->getATP();    
316 >    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
317  
318 +    if (fqa.isFluctuatingCharge()) {
319 +      electrostaticAtomData.is_Fluctuating = true;
320 +      electrostaticAtomData.electronegativity = fqa.getElectronegativity();
321 +      electrostaticAtomData.hardness = fqa.getHardness();
322 +      electrostaticAtomData.slaterN = fqa.getSlaterN();
323 +      electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
324 +    }
325 +
326      pair<map<int,AtomType*>::iterator,bool> ret;    
327 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
327 >    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
328 >                                                        atomType) );
329      if (ret.second == false) {
330        sprintf( painCave.errMsg,
331                 "Electrostatic already had a previous entry with ident %d\n",
332 <               atp.ident);
332 >               atomType->getIdent() );
333        painCave.severity = OPENMD_INFO;
334        painCave.isFatal = 0;
335        simError();        
336      }
337      
338 <    ElectrostaticMap[atomType] = electrostaticAtomData;    
338 >    ElectrostaticMap[atomType] = electrostaticAtomData;  
339 >
340 >    // Now, iterate over all known types and add to the mixing map:
341 >    
342 >    map<AtomType*, ElectrostaticAtomData>::iterator it;
343 >    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
344 >      AtomType* atype2 = (*it).first;
345 >      ElectrostaticAtomData eaData2 = (*it).second;
346 >      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
347 >        
348 >        RealType a = electrostaticAtomData.slaterZeta;
349 >        RealType b = eaData2.slaterZeta;
350 >        int m = electrostaticAtomData.slaterN;
351 >        int n = eaData2.slaterN;
352 >
353 >        // Create the spline of the coulombic integral for s-type
354 >        // Slater orbitals.  Add a 2 angstrom safety window to deal
355 >        // with cutoffGroups that have charged atoms longer than the
356 >        // cutoffRadius away from each other.
357 >
358 >        RealType rval;
359 >        RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
360 >        vector<RealType> rvals;
361 >        vector<RealType> J1vals;
362 >        vector<RealType> J2vals;
363 >        for (int i = 0; i < np_; i++) {
364 >          rval = RealType(i) * dr;
365 >          rvals.push_back(rval);
366 >          J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
367 >          // may not be necessary if Slater coulomb integral is symmetric
368 >          J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
369 >        }
370 >
371 >        CubicSpline* J1 = new CubicSpline();
372 >        J1->addPoints(rvals, J1vals);
373 >        CubicSpline* J2 = new CubicSpline();
374 >        J2->addPoints(rvals, J2vals);
375 >        
376 >        pair<AtomType*, AtomType*> key1, key2;
377 >        key1 = make_pair(atomType, atype2);
378 >        key2 = make_pair(atype2, atomType);
379 >        
380 >        Jij[key1] = J1;
381 >        Jij[key2] = J2;
382 >      }
383 >    }
384 >
385      return;
386    }
387    
# Line 454 | Line 431 | namespace OpenMD {
431      RealType c1, c2, c3, c4;
432      RealType erfcVal(1.0), derfcVal(0.0);
433      RealType BigR;
434 +    RealType two(2.0), three(3.0);
435  
436      Vector3d Q_i, Q_j;
437      Vector3d ux_i, uy_i, uz_i;
# Line 469 | Line 447 | namespace OpenMD {
447      Vector3d indirect_dVdr(V3Zero);
448      Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
449  
450 +    RealType coulInt, vFluc1(0.0), vFluc2(0.0);
451      pair<RealType, RealType> res;
452      
453 +    // splines for coulomb integrals
454 +    CubicSpline* J1;
455 +    CubicSpline* J2;
456 +    
457      if (!initialized_) initialize();
458      
459      ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
# Line 487 | Line 470 | namespace OpenMD {
470      bool i_is_Dipole = data1.is_Dipole;
471      bool i_is_SplitDipole = data1.is_SplitDipole;
472      bool i_is_Quadrupole = data1.is_Quadrupole;
473 +    bool i_is_Fluctuating = data1.is_Fluctuating;
474  
475      bool j_is_Charge = data2.is_Charge;
476      bool j_is_Dipole = data2.is_Dipole;
477      bool j_is_SplitDipole = data2.is_SplitDipole;
478      bool j_is_Quadrupole = data2.is_Quadrupole;
479 +    bool j_is_Fluctuating = data2.is_Fluctuating;
480      
481      if (i_is_Charge) {
482 <      q_i = data1.charge;
482 >      q_i = data1.fixedCharge;
483 >
484 >      if (i_is_Fluctuating) {
485 >        q_i += *(idat.flucQ1);
486 >      }
487 >      
488        if (idat.excluded) {
489          *(idat.skippedCharge2) += q_i;
490        }
# Line 532 | Line 522 | namespace OpenMD {
522      }
523  
524      if (j_is_Charge) {
525 <      q_j = data2.charge;
525 >      q_j = data2.fixedCharge;
526 >
527 >      if (i_is_Fluctuating)
528 >        q_j += *(idat.flucQ2);
529 >
530        if (idat.excluded) {
531          *(idat.skippedCharge1) += q_j;
532        }
# Line 570 | Line 564 | namespace OpenMD {
564        duduz_j = V3Zero;
565      }
566      
567 +    if (i_is_Fluctuating && j_is_Fluctuating) {
568 +      J1 = Jij[idat.atypes];
569 +      J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
570 +    }
571 +
572      epot = 0.0;
573      dVdr = V3Zero;
574      
# Line 578 | Line 577 | namespace OpenMD {
577        if (j_is_Charge) {
578          if (screeningMethod_ == DAMPED) {
579            // assemble the damping variables
580 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
581 <          erfcVal = res.first;
582 <          derfcVal = res.second;
580 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
581 >          //erfcVal = res.first;
582 >          //derfcVal = res.second;
583 >
584 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
585 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
586 >
587            c1 = erfcVal * riji;
588            c2 = (-derfcVal + c1) * riji;
589          } else {
# Line 610 | Line 613 | namespace OpenMD {
613            if (idat.excluded) {
614              indirect_vpair += preVal * rfVal;
615              indirect_Pot += *(idat.sw) * preVal * rfVal;
616 <            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
616 >            indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
617            }
618            
619          } else {
620  
621            vterm = preVal * riji * erfcVal;          
622            dudr  = -  *(idat.sw)  * preVal * c2;
623 <
623 >          
624          }
625 <
625 >        
626          vpair += vterm;
627          epot +=  *(idat.sw)  * vterm;
628 <        dVdr += dudr * rhat;                
628 >        dVdr += dudr * rhat;
629 >
630 >        if (i_is_Fluctuating) {
631 >          if (idat.excluded) {
632 >            // vFluc1 is the difference between the direct coulomb integral
633 >            // and the normal 1/r-like  interaction between point charges.
634 >            coulInt = J1->getValueAt( *(idat.rij) );
635 >            vFluc1 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
636 >          } else {
637 >            vFluc1 = 0.0;
638 >          }
639 >          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) / q_i;
640 >        }
641 >
642 >        if (j_is_Fluctuating) {
643 >          if (idat.excluded) {
644 >            // vFluc2 is the difference between the direct coulomb integral
645 >            // and the normal 1/r-like  interaction between point charges.
646 >            coulInt = J2->getValueAt( *(idat.rij) );
647 >            vFluc2 = pre11_ * coulInt * q_i * q_j  - (*(idat.sw) * vterm);
648 >          } else {
649 >            vFluc2 = 0.0;
650 >          }
651 >          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) / q_j;
652 >        }
653 >          
654 >
655        }
656  
657        if (j_is_Dipole) {
# Line 638 | Line 667 | namespace OpenMD {
667            vpair += vterm;
668            epot +=  *(idat.sw)  * vterm;
669  
670 <          dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
670 >          dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
671            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
672  
673            // Even if we excluded this pair from direct interactions,
# Line 667 | Line 696 | namespace OpenMD {
696  
697            if (screeningMethod_ == DAMPED) {
698              // assemble the damping variables
699 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
700 <            erfcVal = res.first;
701 <            derfcVal = res.second;
699 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
700 >            //erfcVal = res.first;
701 >            //derfcVal = res.second;
702 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
703 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
704              c1 = erfcVal * ri;
705              c2 = (-derfcVal + c1) * ri;
706              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 693 | Line 724 | namespace OpenMD {
724            duduz_j += -preSw * pot_term * rhat;
725  
726          }
727 +        if (i_is_Fluctuating) {
728 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
729 +        }
730        }
731  
732        if (j_is_Quadrupole) {
# Line 704 | Line 738 | namespace OpenMD {
738            
739          if (screeningMethod_ == DAMPED) {
740            // assemble the damping variables
741 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
742 <          erfcVal = res.first;
743 <          derfcVal = res.second;
741 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
742 >          //erfcVal = res.first;
743 >          //derfcVal = res.second;
744 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
745 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
746            c1 = erfcVal * riji;
747            c2 = (-derfcVal + c1) * riji;
748            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 723 | Line 759 | namespace OpenMD {
759          c2ri = c2 * riji;
760          c3ri = c3 * riji;
761          c4rij = c4 *  *(idat.rij) ;
762 <        rhatdot2 = 2.0 * rhat * c3;
762 >        rhatdot2 = two * rhat * c3;
763          rhatc4 = rhat * c4rij;
764  
765          // calculate the potential
# Line 736 | Line 772 | namespace OpenMD {
772                  
773          // calculate derivatives for the forces and torques
774  
775 <        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
776 <                           qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
777 <                           qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
775 >        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
776 >                           qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
777 >                           qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
778                            
779          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
780          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
781          duduz_j += preSw * qzz_j * cz_j * rhatdot2;
782 +        if (i_is_Fluctuating) {
783 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
784 +        }
785 +
786        }
787      }
788      
# Line 762 | Line 802 | namespace OpenMD {
802            vpair += vterm;
803            epot +=  *(idat.sw)  * vterm;
804            
805 <          dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
805 >          dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
806            
807            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
808  
# Line 793 | Line 833 | namespace OpenMD {
833              
834            if (screeningMethod_ == DAMPED) {
835              // assemble the damping variables
836 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
837 <            erfcVal = res.first;
838 <            derfcVal = res.second;
836 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
837 >            //erfcVal = res.first;
838 >            //derfcVal = res.second;
839 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
840 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
841              c1 = erfcVal * ri;
842              c2 = (-derfcVal + c1) * ri;
843              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 816 | Line 858 | namespace OpenMD {
858            // calculate derivatives for the forces and torques
859            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
860            duduz_i += preSw * pot_term * rhat;
861 +        }
862 +
863 +        if (j_is_Fluctuating) {
864 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
865          }
866 +
867        }
868  
869        if (j_is_Dipole) {
# Line 838 | Line 885 | namespace OpenMD {
885              
886            a1 = 5.0 * ct_i * ct_j - ct_ij;
887              
888 <          dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
888 >          dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
889  
890 <          duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
891 <          duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
890 >          duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
891 >          duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
892  
893            if (idat.excluded) {
894              indirect_vpair +=  - pref * preRF2_ * ct_ij;
# Line 872 | Line 919 | namespace OpenMD {
919            }
920            if (screeningMethod_ == DAMPED) {
921              // assemble damping variables
922 <            res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
923 <            erfcVal = res.first;
924 <            derfcVal = res.second;
922 >            //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
923 >            //erfcVal = res.first;
924 >            //derfcVal = res.second;
925 >            erfcVal = erfc(dampingAlpha_ * *(idat.rij));
926 >            derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
927              c1 = erfcVal * ri;
928              c2 = (-derfcVal + c1) * ri;
929              c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
# Line 923 | Line 972 | namespace OpenMD {
972  
973          if (screeningMethod_ == DAMPED) {
974            // assemble the damping variables
975 <          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
976 <          erfcVal = res.first;
977 <          derfcVal = res.second;
975 >          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
976 >          //erfcVal = res.first;
977 >          //derfcVal = res.second;
978 >          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
979 >          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
980            c1 = erfcVal * riji;
981            c2 = (-derfcVal + c1) * riji;
982            c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
# Line 942 | Line 993 | namespace OpenMD {
993          c2ri = c2 * riji;
994          c3ri = c3 * riji;
995          c4rij = c4 *  *(idat.rij) ;
996 <        rhatdot2 = 2.0 * rhat * c3;
996 >        rhatdot2 = two * rhat * c3;
997          rhatc4 = rhat * c4rij;
998  
999          // calculate the potential
# Line 956 | Line 1007 | namespace OpenMD {
1007  
1008          // calculate the derivatives for the forces and torques
1009  
1010 <        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
1011 <                          qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
1012 <                          qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
1010 >        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
1011 >                          qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
1012 >                          qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
1013  
1014          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1015          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1016          duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1017 +
1018 +        if (j_is_Fluctuating) {
1019 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1020 +        }
1021 +
1022        }
1023      }
1024  
# Line 990 | Line 1046 | namespace OpenMD {
1046  
1047        // only accumulate the forces and torques resulting from the
1048        // indirect reaction field terms.
1049 +
1050        *(idat.vpair) += indirect_vpair;
1051        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1052        *(idat.f1) += indirect_dVdr;
# Line 1000 | Line 1057 | namespace OpenMD {
1057          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1058      }
1059  
1003 –
1060      return;
1061    }  
1062      
# Line 1030 | Line 1086 | namespace OpenMD {
1086        }
1087      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1088        if (i_is_Charge) {        
1089 <        chg1 = data.charge;
1089 >        chg1 = data.fixedCharge;
1090          if (screeningMethod_ == DAMPED) {
1091            self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1092          } else {        

Diff Legend

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