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 1668 by gezelter, Fri Jan 6 19:03:05 2012 UTC vs.
Revision 1767 by gezelter, Fri Jul 6 22:01:58 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 + #include "math/erfc.hpp"
57  
58   namespace OpenMD {
59    
# 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 >        // don't start at i = 0, as rval = 0 is undefined for the slater overlap integrals.
365 >        for (int i = 1; i < np_+1; i++) {
366 >          rval = RealType(i) * dr;
367 >          rvals.push_back(rval);
368 >          J1vals.push_back(sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
369 >          // may not be necessary if Slater coulomb integral is symmetric
370 >          J2vals.push_back(sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal );
371 >        }
372 >
373 >        CubicSpline* J1 = new CubicSpline();
374 >        J1->addPoints(rvals, J1vals);
375 >        CubicSpline* J2 = new CubicSpline();
376 >        J2->addPoints(rvals, J2vals);
377 >        
378 >        pair<AtomType*, AtomType*> key1, key2;
379 >        key1 = make_pair(atomType, atype2);
380 >        key2 = make_pair(atype2, atomType);
381 >        
382 >        Jij[key1] = J1;
383 >        Jij[key2] = J2;
384 >      }
385 >    }
386 >
387      return;
388    }
389    
# Line 477 | Line 449 | namespace OpenMD {
449      Vector3d indirect_dVdr(V3Zero);
450      Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
451  
452 +    RealType coulInt, vFluc1(0.0), vFluc2(0.0);
453      pair<RealType, RealType> res;
454      
455 +    // splines for coulomb integrals
456 +    CubicSpline* J1;
457 +    CubicSpline* J2;
458 +    
459      if (!initialized_) initialize();
460      
461      ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
# Line 495 | Line 472 | namespace OpenMD {
472      bool i_is_Dipole = data1.is_Dipole;
473      bool i_is_SplitDipole = data1.is_SplitDipole;
474      bool i_is_Quadrupole = data1.is_Quadrupole;
475 +    bool i_is_Fluctuating = data1.is_Fluctuating;
476  
477      bool j_is_Charge = data2.is_Charge;
478      bool j_is_Dipole = data2.is_Dipole;
479      bool j_is_SplitDipole = data2.is_SplitDipole;
480      bool j_is_Quadrupole = data2.is_Quadrupole;
481 +    bool j_is_Fluctuating = data2.is_Fluctuating;
482      
483      if (i_is_Charge) {
484 <      q_i = data1.charge;
484 >      q_i = data1.fixedCharge;
485 >
486 >      if (i_is_Fluctuating) {
487 >        q_i += *(idat.flucQ1);
488 >      }
489 >      
490        if (idat.excluded) {
491          *(idat.skippedCharge2) += q_i;
492        }
# Line 540 | Line 524 | namespace OpenMD {
524      }
525  
526      if (j_is_Charge) {
527 <      q_j = data2.charge;
527 >      q_j = data2.fixedCharge;
528 >
529 >      if (j_is_Fluctuating)
530 >        q_j += *(idat.flucQ2);
531 >
532        if (idat.excluded) {
533          *(idat.skippedCharge1) += q_j;
534        }
# Line 578 | Line 566 | namespace OpenMD {
566        duduz_j = V3Zero;
567      }
568      
569 +    if (i_is_Fluctuating && j_is_Fluctuating) {
570 +      J1 = Jij[idat.atypes];
571 +      J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)];
572 +    }
573 +
574      epot = 0.0;
575      dVdr = V3Zero;
576      
# Line 586 | Line 579 | namespace OpenMD {
579        if (j_is_Charge) {
580          if (screeningMethod_ == DAMPED) {
581            // assemble the damping variables
582 <          //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
583 <          //erfcVal = res.first;
584 <          //derfcVal = res.second;
582 >          res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
583 >          erfcVal = res.first;
584 >          derfcVal = res.second;
585  
586 <          erfcVal = erfc(dampingAlpha_ * *(idat.rij));
587 <          derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
586 >          //erfcVal = erfc(dampingAlpha_ * *(idat.rij));
587 >          //derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
588  
589            c1 = erfcVal * riji;
590            c2 = (-derfcVal + c1) * riji;
# Line 600 | Line 593 | namespace OpenMD {
593            c2 = c1 * riji;
594          }
595  
596 <        preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
596 >        preVal =  *(idat.electroMult) * pre11_;
597          
598          if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
599            vterm = preVal * (c1 - c1c_);
# Line 629 | Line 622 | namespace OpenMD {
622  
623            vterm = preVal * riji * erfcVal;          
624            dudr  = -  *(idat.sw)  * preVal * c2;
625 +          
626 +        }
627 +        
628 +        vpair += vterm * q_i * q_j;
629 +        epot +=  *(idat.sw)  * vterm * q_i * q_j;
630 +        dVdr += dudr * rhat * q_i * q_j;
631  
632 +        if (i_is_Fluctuating) {
633 +          if (idat.excluded) {
634 +            // vFluc1 is the difference between the direct coulomb integral
635 +            // and the normal 1/r-like  interaction between point charges.
636 +            coulInt = J1->getValueAt( *(idat.rij) );
637 +            vFluc1 = coulInt - (*(idat.sw) * vterm);
638 +          } else {
639 +            vFluc1 = 0.0;
640 +          }
641 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j;
642          }
643  
644 <        vpair += vterm;
645 <        epot +=  *(idat.sw)  * vterm;
646 <        dVdr += dudr * rhat;                
644 >        if (j_is_Fluctuating) {
645 >          if (idat.excluded) {
646 >            // vFluc2 is the difference between the direct coulomb integral
647 >            // and the normal 1/r-like  interaction between point charges.
648 >            coulInt = J2->getValueAt( *(idat.rij) );
649 >            vFluc2 = coulInt - (*(idat.sw) * vterm);
650 >          } else {
651 >            vFluc2 = 0.0;
652 >          }
653 >          *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i;
654 >        }
655 >          
656 >
657        }
658  
659        if (j_is_Dipole) {
# Line 707 | Line 726 | namespace OpenMD {
726            duduz_j += -preSw * pot_term * rhat;
727  
728          }
729 +        if (i_is_Fluctuating) {
730 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
731 +        }
732        }
733  
734        if (j_is_Quadrupole) {
# Line 759 | Line 781 | namespace OpenMD {
781          dudux_j += preSw * qxx_j * cx_j * rhatdot2;
782          duduy_j += preSw * qyy_j * cy_j * rhatdot2;
783          duduz_j += preSw * qzz_j * cz_j * rhatdot2;
784 +        if (i_is_Fluctuating) {
785 +          *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i;
786 +        }
787 +
788        }
789      }
790      
# Line 835 | Line 861 | namespace OpenMD {
861            dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
862            duduz_i += preSw * pot_term * rhat;
863          }
864 +
865 +        if (j_is_Fluctuating) {
866 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
867 +        }
868 +
869        }
870  
871        if (j_is_Dipole) {
# Line 985 | Line 1016 | namespace OpenMD {
1016          dudux_i += preSw * qxx_i * cx_i *  rhatdot2;
1017          duduy_i += preSw * qyy_i * cy_i *  rhatdot2;
1018          duduz_i += preSw * qzz_i * cz_i *  rhatdot2;
1019 +
1020 +        if (j_is_Fluctuating) {
1021 +          *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j;
1022 +        }
1023 +
1024        }
1025      }
1026  
# Line 1014 | Line 1050 | namespace OpenMD {
1050        // indirect reaction field terms.
1051  
1052        *(idat.vpair) += indirect_vpair;
1053 +      
1054 +      (*(idat.excludedPot))[ELECTROSTATIC_FAMILY] +=   (*(idat.sw) * vterm +
1055 +                                                        vFluc1 ) * q_i * q_j;
1056        (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1057        *(idat.f1) += indirect_dVdr;
1058        
# Line 1023 | Line 1062 | namespace OpenMD {
1062          *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1063      }
1064  
1026 –
1065      return;
1066    }  
1067      
1068    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1069 <    RealType mu1, preVal, chg1, self;
1032 <    
1069 >    RealType mu1, preVal, self;
1070      if (!initialized_) initialize();
1071  
1072      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
# Line 1037 | Line 1074 | namespace OpenMD {
1074      // logicals
1075      bool i_is_Charge = data.is_Charge;
1076      bool i_is_Dipole = data.is_Dipole;
1077 +    bool i_is_Fluctuating = data.is_Fluctuating;
1078 +    RealType chg1 = data.fixedCharge;  
1079 +    
1080 +    if (i_is_Fluctuating) {
1081 +      chg1 += *(sdat.flucQ);
1082 +      // dVdFQ is really a force, so this is negative the derivative
1083 +      *(sdat.dVdFQ) -=  *(sdat.flucQ) * data.hardness + data.electronegativity;
1084 +      (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1085 +        (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1086 +    }
1087  
1088      if (summationMethod_ == esm_REACTION_FIELD) {
1089        if (i_is_Dipole) {
# Line 1053 | Line 1100 | namespace OpenMD {
1100        }
1101      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1102        if (i_is_Charge) {        
1056 –        chg1 = data.charge;
1103          if (screeningMethod_ == DAMPED) {
1104            self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_;
1105          } else {        

Diff Legend

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