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 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC vs.
Revision 1750 by gezelter, Thu Jun 7 12:53:46 2012 UTC

# Line 47 | 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 188 | Line 193 | namespace OpenMD {
193          
194          // throw warning
195          sprintf( painCave.errMsg,
196 <                 "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
197 <                 "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
196 >                 "Electrostatic::initialize: dampingAlpha was not specified in the\n"
197 >                 "\tinput file.  A default value of %f (1/ang) will be used for the\n"
198 >                 "\tcutoff of %f (ang).\n",
199                   dampingAlpha_, cutoffRadius_);
200          painCave.severity = OPENMD_INFO;
201          painCave.isFatal = 0;
# Line 212 | Line 218 | namespace OpenMD {
218          addType(at);
219      }
220      
215 –
221      cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
222      rcuti_ = 1.0 / cutoffRadius_;
223      rcuti2_ = rcuti_ * rcuti_;
# Line 279 | Line 284 | namespace OpenMD {
284      electrostaticAtomData.is_Dipole = false;
285      electrostaticAtomData.is_SplitDipole = false;
286      electrostaticAtomData.is_Quadrupole = false;
287 +    electrostaticAtomData.is_Fluctuating = false;
288  
289 <    if (atomType->isCharge()) {
284 <      GenericData* data = atomType->getPropertyByName("Charge");
289 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
290  
291 <      if (data == NULL) {
287 <        sprintf( painCave.errMsg, "Electrostatic::addType could not find "
288 <                 "Charge\n"
289 <                 "\tparameters for atomType %s.\n",
290 <                 atomType->getName().c_str());
291 <        painCave.severity = OPENMD_ERROR;
292 <        painCave.isFatal = 1;
293 <        simError();                  
294 <      }
295 <      
296 <      DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
297 <      if (doubleData == NULL) {
298 <        sprintf( painCave.errMsg,
299 <                 "Electrostatic::addType could not convert GenericData to "
300 <                 "Charge for\n"
301 <                 "\tatom type %s\n", atomType->getName().c_str());
302 <        painCave.severity = OPENMD_ERROR;
303 <        painCave.isFatal = 1;
304 <        simError();          
305 <      }
291 >    if (fca.isFixedCharge()) {
292        electrostaticAtomData.is_Charge = true;
293 <      electrostaticAtomData.charge = doubleData->getData();          
293 >      electrostaticAtomData.fixedCharge = fca.getCharge();
294      }
295  
296 <    if (atomType->isDirectional()) {
297 <      DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
298 <      
313 <      if (daType->isDipole()) {
314 <        GenericData* data = daType->getPropertyByName("Dipole");
315 <        
316 <        if (data == NULL) {
317 <          sprintf( painCave.errMsg,
318 <                   "Electrostatic::addType could not find Dipole\n"
319 <                   "\tparameters for atomType %s.\n",
320 <                   daType->getName().c_str());
321 <          painCave.severity = OPENMD_ERROR;
322 <          painCave.isFatal = 1;
323 <          simError();                  
324 <        }
325 <      
326 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
327 <        if (doubleData == NULL) {
328 <          sprintf( painCave.errMsg,
329 <                   "Electrostatic::addType could not convert GenericData to "
330 <                   "Dipole Moment\n"
331 <                   "\tfor atom type %s\n", daType->getName().c_str());
332 <          painCave.severity = OPENMD_ERROR;
333 <          painCave.isFatal = 1;
334 <          simError();          
335 <        }
296 >    MultipoleAdapter ma = MultipoleAdapter(atomType);
297 >    if (ma.isMultipole()) {
298 >      if (ma.isDipole()) {
299          electrostaticAtomData.is_Dipole = true;
300 <        electrostaticAtomData.dipole_moment = doubleData->getData();
300 >        electrostaticAtomData.dipole_moment = ma.getDipoleMoment();
301        }
302 <
340 <      if (daType->isSplitDipole()) {
341 <        GenericData* data = daType->getPropertyByName("SplitDipoleDistance");
342 <        
343 <        if (data == NULL) {
344 <          sprintf(painCave.errMsg,
345 <                  "Electrostatic::addType could not find SplitDipoleDistance\n"
346 <                  "\tparameter for atomType %s.\n",
347 <                  daType->getName().c_str());
348 <          painCave.severity = OPENMD_ERROR;
349 <          painCave.isFatal = 1;
350 <          simError();                  
351 <        }
352 <      
353 <        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
354 <        if (doubleData == NULL) {
355 <          sprintf( painCave.errMsg,
356 <                   "Electrostatic::addType could not convert GenericData to "
357 <                   "SplitDipoleDistance for\n"
358 <                   "\tatom type %s\n", daType->getName().c_str());
359 <          painCave.severity = OPENMD_ERROR;
360 <          painCave.isFatal = 1;
361 <          simError();          
362 <        }
302 >      if (ma.isSplitDipole()) {
303          electrostaticAtomData.is_SplitDipole = true;
304 <        electrostaticAtomData.split_dipole_distance = doubleData->getData();
304 >        electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance();
305        }
306 <
367 <      if (daType->isQuadrupole()) {
368 <        GenericData* data = daType->getPropertyByName("QuadrupoleMoments");
369 <        
370 <        if (data == NULL) {
371 <          sprintf( painCave.errMsg,
372 <                   "Electrostatic::addType could not find QuadrupoleMoments\n"
373 <                   "\tparameter for atomType %s.\n",
374 <                   daType->getName().c_str());
375 <          painCave.severity = OPENMD_ERROR;
376 <          painCave.isFatal = 1;
377 <          simError();                  
378 <        }
379 <        
306 >      if (ma.isQuadrupole()) {
307          // Quadrupoles in OpenMD are set as the diagonal elements
308          // of the diagonalized traceless quadrupole moment tensor.
309          // The column vectors of the unitary matrix that diagonalizes
310          // the quadrupole moment tensor become the eFrame (or the
311          // electrostatic version of the body-fixed frame.
385 –
386 –        Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
387 –        if (v3dData == NULL) {
388 –          sprintf( painCave.errMsg,
389 –                   "Electrostatic::addType could not convert GenericData to "
390 –                   "Quadrupole Moments for\n"
391 –                   "\tatom type %s\n", daType->getName().c_str());
392 –          painCave.severity = OPENMD_ERROR;
393 –          painCave.isFatal = 1;
394 –          simError();          
395 –        }
312          electrostaticAtomData.is_Quadrupole = true;
313 <        electrostaticAtomData.quadrupole_moments = v3dData->getData();
313 >        electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments();
314        }
315      }
316      
317 <    AtomTypeProperties atp = atomType->getATP();    
317 >    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
318  
319 +    if (fqa.isFluctuatingCharge()) {
320 +      electrostaticAtomData.is_Fluctuating = true;
321 +      electrostaticAtomData.electronegativity = fqa.getElectronegativity();
322 +      electrostaticAtomData.hardness = fqa.getHardness();
323 +      electrostaticAtomData.slaterN = fqa.getSlaterN();
324 +      electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
325 +    }
326 +
327      pair<map<int,AtomType*>::iterator,bool> ret;    
328 <    ret = ElectrostaticList.insert( pair<int,AtomType*>(atp.ident, atomType) );
328 >    ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
329 >                                                        atomType) );
330      if (ret.second == false) {
331        sprintf( painCave.errMsg,
332                 "Electrostatic already had a previous entry with ident %d\n",
333 <               atp.ident);
333 >               atomType->getIdent() );
334        painCave.severity = OPENMD_INFO;
335        painCave.isFatal = 0;
336        simError();        
337      }
338      
339 <    ElectrostaticMap[atomType] = electrostaticAtomData;    
339 >    ElectrostaticMap[atomType] = electrostaticAtomData;  
340 >
341 >    // Now, iterate over all known types and add to the mixing map:
342 >    
343 >    map<AtomType*, ElectrostaticAtomData>::iterator it;
344 >    for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
345 >      AtomType* atype2 = (*it).first;
346 >      ElectrostaticAtomData eaData2 = (*it).second;
347 >      if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
348 >        
349 >        RealType a = electrostaticAtomData.slaterZeta;
350 >        RealType b = eaData2.slaterZeta;
351 >        int m = electrostaticAtomData.slaterN;
352 >        int n = eaData2.slaterN;
353 >
354 >        // Create the spline of the coulombic integral for s-type
355 >        // Slater orbitals.  Add a 2 angstrom safety window to deal
356 >        // with cutoffGroups that have charged atoms longer than the
357 >        // cutoffRadius away from each other.
358 >
359 >        RealType rval;
360 >        RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
361 >        vector<RealType> rvals;
362 >        vector<RealType> J1vals;
363 >        vector<RealType> J2vals;
364 >        for (int i = 0; i < np_; i++) {
365 >          rval = RealType(i) * dr;
366 >          rvals.push_back(rval);
367 >          J1vals.push_back(electrostaticAtomData.hardness * sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) );
368 >          // may not be necessary if Slater coulomb integral is symmetric
369 >          J2vals.push_back(eaData2.hardness *  sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) );
370 >        }
371 >
372 >        CubicSpline* J1 = new CubicSpline();
373 >        J1->addPoints(rvals, J1vals);
374 >        CubicSpline* J2 = new CubicSpline();
375 >        J2->addPoints(rvals, J2vals);
376 >        
377 >        pair<AtomType*, AtomType*> key1, key2;
378 >        key1 = make_pair(atomType, atype2);
379 >        key2 = make_pair(atype2, atomType);
380 >        
381 >        Jij[key1] = J1;
382 >        Jij[key2] = J2;
383 >      }
384 >    }
385 >
386      return;
387    }
388    
# Line 461 | Line 432 | namespace OpenMD {
432      RealType c1, c2, c3, c4;
433      RealType erfcVal(1.0), derfcVal(0.0);
434      RealType BigR;
435 +    RealType two(2.0), three(3.0);
436  
437      Vector3d Q_i, Q_j;
438      Vector3d ux_i, uy_i, uz_i;
# Line 476 | Line 448 | namespace OpenMD {
448      Vector3d indirect_dVdr(V3Zero);
449      Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
450  
451 +    RealType coulInt, vFluc1(0.0), vFluc2(0.0);
452      pair<RealType, RealType> res;
453      
454 +    // splines for coulomb integrals
455 +    CubicSpline* J1;
456 +    CubicSpline* J2;
457 +    
458      if (!initialized_) initialize();
459      
460      ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
# Line 494 | Line 471 | namespace OpenMD {
471      bool i_is_Dipole = data1.is_Dipole;
472      bool i_is_SplitDipole = data1.is_SplitDipole;
473      bool i_is_Quadrupole = data1.is_Quadrupole;
474 +    bool i_is_Fluctuating = data1.is_Fluctuating;
475  
476      bool j_is_Charge = data2.is_Charge;
477      bool j_is_Dipole = data2.is_Dipole;
478      bool j_is_SplitDipole = data2.is_SplitDipole;
479      bool j_is_Quadrupole = data2.is_Quadrupole;
480 +    bool j_is_Fluctuating = data2.is_Fluctuating;
481      
482      if (i_is_Charge) {
483 <      q_i = data1.charge;
483 >      q_i = data1.fixedCharge;
484 >
485 >      if (i_is_Fluctuating) {
486 >        q_i += *(idat.flucQ1);
487 >      }
488 >      
489        if (idat.excluded) {
490          *(idat.skippedCharge2) += q_i;
491        }
# Line 539 | Line 523 | namespace OpenMD {
523      }
524  
525      if (j_is_Charge) {
526 <      q_j = data2.charge;
526 >      q_j = data2.fixedCharge;
527 >
528 >      if (j_is_Fluctuating)
529 >        q_j += *(idat.flucQ2);
530 >
531        if (idat.excluded) {
532          *(idat.skippedCharge1) += q_j;
533        }
# Line 577 | Line 565 | namespace OpenMD {
565        duduz_j = V3Zero;
566      }
567      
568 +    if (i_is_Fluctuating && j_is_Fluctuating) {
569 +      J1 = Jij[idat.atypes];
570 +      J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
571 +    }
572 +
573      epot = 0.0;
574      dVdr = V3Zero;
575      
# Line 599 | Line 592 | namespace OpenMD {
592            c2 = c1 * riji;
593          }
594  
595 <        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
595 >        preVal =  *(idat.electroMult) * pre11_;
596          
597          if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
598            vterm = preVal * (c1 - c1c_);
# Line 621 | Line 614 | namespace OpenMD {
614            if (idat.excluded) {
615              indirect_vpair += preVal * rfVal;
616              indirect_Pot += *(idat.sw) * preVal * rfVal;
617 <            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
617 >            indirect_dVdr += *(idat.sw)  * preVal * two * rfVal  * riji * rhat;
618            }
619            
620          } else {
621  
622            vterm = preVal * riji * erfcVal;          
623            dudr  = -  *(idat.sw)  * preVal * c2;
624 +          
625 +        }
626 +        
627 +        vpair += vterm * q_i * q_j;
628 +        epot +=  *(idat.sw)  * vterm * q_i * q_j;
629 +        dVdr += dudr * rhat * q_i * q_j;
630  
631 +        if (i_is_Fluctuating) {
632 +          if (idat.excluded) {
633 +            // vFluc1 is the difference between the direct coulomb integral
634 +            // and the normal 1/r-like  interaction between point charges.
635 +            coulInt = J1->getValueAt( *(idat.rij) );
636 +            vFluc1 = coulInt - (*(idat.sw) * vterm);
637 +          } else {
638 +            vFluc1 = 0.0;
639 +          }
640 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j;
641          }
642  
643 <        vpair += vterm;
644 <        epot +=  *(idat.sw)  * vterm;
645 <        dVdr += dudr * rhat;                
643 >        if (j_is_Fluctuating) {
644 >          if (idat.excluded) {
645 >            // vFluc2 is the difference between the direct coulomb integral
646 >            // and the normal 1/r-like  interaction between point charges.
647 >            coulInt = J2->getValueAt( *(idat.rij) );
648 >            vFluc2 = coulInt - (*(idat.sw) * vterm);
649 >          } else {
650 >            vFluc2 = 0.0;
651 >          }
652 >          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i;
653 >        }
654 >          
655 >
656        }
657  
658        if (j_is_Dipole) {
# Line 649 | Line 668 | namespace OpenMD {
668            vpair += vterm;
669            epot +=  *(idat.sw)  * vterm;
670  
671 <          dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
671 >          dVdr +=  -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
672            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
673  
674            // Even if we excluded this pair from direct interactions,
# Line 706 | Line 725 | namespace OpenMD {
725            duduz_j += -preSw * pot_term * rhat;
726  
727          }
728 +        if (i_is_Fluctuating) {
729 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
730 +        }
731        }
732  
733        if (j_is_Quadrupole) {
# Line 738 | Line 760 | namespace OpenMD {
760          c2ri = c2 * riji;
761          c3ri = c3 * riji;
762          c4rij = c4 *  *(idat.rij) ;
763 <        rhatdot2 = 2.0 * rhat * c3;
763 >        rhatdot2 = two * rhat * c3;
764          rhatc4 = rhat * c4rij;
765  
766          // calculate the potential
# Line 751 | Line 773 | namespace OpenMD {
773                  
774          // calculate derivatives for the forces and torques
775  
776 <        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) +
777 <                           qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) +
778 <                           qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri));
776 >        dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) +
777 >                           qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) +
778 >                           qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri));
779                            
780          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
781          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
782          duduz_j += preSw * qzz_j * cz_j * rhatdot2;
783 +        if (i_is_Fluctuating) {
784 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
785 +        }
786 +
787        }
788      }
789      
# Line 777 | Line 803 | namespace OpenMD {
803            vpair += vterm;
804            epot +=  *(idat.sw)  * vterm;
805            
806 <          dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
806 >          dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i);
807            
808            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
809  
# Line 833 | Line 859 | namespace OpenMD {
859            // calculate derivatives for the forces and torques
860            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
861            duduz_i += preSw * pot_term * rhat;
862 +        }
863 +
864 +        if (j_is_Fluctuating) {
865 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
866          }
867 +
868        }
869  
870        if (j_is_Dipole) {
# Line 855 | Line 886 | namespace OpenMD {
886              
887            a1 = 5.0 * ct_i * ct_j - ct_ij;
888              
889 <          dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
889 >          dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i);
890  
891 <          duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
892 <          duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
891 >          duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j);
892 >          duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i);
893  
894            if (idat.excluded) {
895              indirect_vpair +=  - pref * preRF2_ * ct_ij;
# Line 963 | Line 994 | namespace OpenMD {
994          c2ri = c2 * riji;
995          c3ri = c3 * riji;
996          c4rij = c4 *  *(idat.rij) ;
997 <        rhatdot2 = 2.0 * rhat * c3;
997 >        rhatdot2 = two * rhat * c3;
998          rhatc4 = rhat * c4rij;
999  
1000          // calculate the potential
# Line 977 | Line 1008 | namespace OpenMD {
1008  
1009          // calculate the derivatives for the forces and torques
1010  
1011 <        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) +
1012 <                          qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) +
1013 <                          qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri));
1011 >        dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) +
1012 >                          qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) +
1013 >                          qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri));
1014  
1015          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1016          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1017          duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1018 +
1019 +        if (j_is_Fluctuating) {
1020 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1021 +        }
1022 +
1023        }
1024      }
1025  
# Line 1022 | Line 1058 | namespace OpenMD {
1058          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1059      }
1060  
1025 –
1061      return;
1062    }  
1063      
1064    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1065 <    RealType mu1, preVal, chg1, self;
1031 <    
1065 >    RealType mu1, preVal, self;
1066      if (!initialized_) initialize();
1067  
1068      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 1036 | Line 1070 | namespace OpenMD {
1070      // logicals
1071      bool i_is_Charge = data.is_Charge;
1072      bool i_is_Dipole = data.is_Dipole;
1073 +    bool i_is_Fluctuating = data.is_Fluctuating;
1074 +    RealType chg1 = data.fixedCharge;  
1075 +    
1076 +    if (i_is_Fluctuating) {
1077 +      chg1 += *(sdat.flucQ);
1078 +      // dVdFQ is really a force, so this is negative the derivative
1079 +      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1080 +    }
1081  
1082      if (summationMethod_ == esm_REACTION_FIELD) {
1083        if (i_is_Dipole) {
# Line 1052 | Line 1094 | namespace OpenMD {
1094        }
1095      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1096        if (i_is_Charge) {        
1055 –        chg1 = data.charge;
1097          if (screeningMethod_ == DAMPED) {
1098            self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1099          } else {        

Diff Legend

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