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 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs.
Revision 1613 by gezelter, Thu Aug 18 20:18:19 2011 UTC

# Line 34 | Line 34
34   * work.  Good starting points are:
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 < * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
37 > * [2]  Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39   * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
# Line 52 | Line 52 | namespace OpenMD {
52   namespace OpenMD {
53    
54    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55 <                                  forceField_(NULL) {}
55 >                                  forceField_(NULL), info_(NULL),
56 >                                  haveCutoffRadius_(false),
57 >                                  haveDampingAlpha_(false),
58 >                                  haveDielectric_(false),
59 >                                  haveElectroSpline_(false)
60 >  {}
61    
62    void Electrostatic::initialize() {
63 +    
64 +    Globals* simParams_ = info_->getSimParams();
65  
59 –    Globals* simParams_;
60 –
66      summationMap_["HARD"]               = esm_HARD;
67      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
68      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
# Line 97 | Line 102 | namespace OpenMD {
102      screeningMethod_ = UNDAMPED;
103      dielectric_ = 1.0;
104      one_third_ = 1.0 / 3.0;
100 –    haveCutoffRadius_ = false;
101 –    haveDampingAlpha_ = false;
102 –    haveDielectric_ = false;  
103 –    haveElectroSpline_ = false;
105    
106      // check the summation method:
107      if (simParams_->haveElectrostaticSummationMethod()) {
# Line 248 | Line 249 | namespace OpenMD {
249        preRF2_ = 2.0 * preRF_;
250      }
251      
252 <    RealType dx = cutoffRadius_ / RealType(np_ - 1);
252 >    // Add a 2 angstrom safety window to deal with cutoffGroups that
253 >    // have charged atoms longer than the cutoffRadius away from each
254 >    // other.  Splining may not be the best choice here.  Direct calls
255 >    // to erfc might be preferrable.
256 >
257 >    RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
258      RealType rval;
259      vector<RealType> rvals;
260      vector<RealType> yvals;
# Line 407 | Line 413 | namespace OpenMD {
413      return;
414    }
415    
416 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
417 <                                                    RealType theRSW ) {
412 <    cutoffRadius_ = theECR;
416 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
417 >    cutoffRadius_ = rCut;
418      rrf_ = cutoffRadius_;
414 –    rt_ = theRSW;
419      haveCutoffRadius_ = true;
420    }
421 +
422 +  void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
423 +    rt_ = rSwitch;
424 +  }
425    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
426      summationMethod_ = esm;
427    }
# Line 443 | Line 451 | namespace OpenMD {
451      RealType ct_i, ct_j, ct_ij, a1;
452      RealType riji, ri, ri2, ri3, ri4;
453      RealType pref, vterm, epot, dudr;
454 +    RealType vpair(0.0);
455      RealType scale, sc2;
456      RealType pot_term, preVal, rfVal;
457      RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
458      RealType preSw, preSwSc;
459      RealType c1, c2, c3, c4;
460 <    RealType erfcVal, derfcVal;
460 >    RealType erfcVal(1.0), derfcVal(0.0);
461      RealType BigR;
462  
463      Vector3d Q_i, Q_j;
# Line 459 | Line 468 | namespace OpenMD {
468      Vector3d rhatdot2, rhatc4;
469      Vector3d dVdr;
470  
471 +    // variables for indirect (reaction field) interactions for excluded pairs:
472 +    RealType indirect_Pot(0.0);
473 +    RealType indirect_vpair(0.0);
474 +    Vector3d indirect_dVdr(V3Zero);
475 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
476 +
477      pair<RealType, RealType> res;
478      
479      if (!initialized_) initialize();
# Line 483 | Line 498 | namespace OpenMD {
498      bool j_is_SplitDipole = data2.is_SplitDipole;
499      bool j_is_Quadrupole = data2.is_Quadrupole;
500      
501 <    if (i_is_Charge)
501 >    if (i_is_Charge) {
502        q_i = data1.charge;
503 +      if (idat.excluded) {
504 +        *(idat.skippedCharge2) += q_i;
505 +      }
506 +    }
507  
508      if (i_is_Dipole) {
509        mu_i = data1.dipole_moment;
# Line 517 | Line 536 | namespace OpenMD {
536        duduz_i = V3Zero;
537      }
538  
539 <    if (j_is_Charge)
539 >    if (j_is_Charge) {
540        q_j = data2.charge;
541 +      if (idat.excluded) {
542 +        *(idat.skippedCharge1) += q_j;
543 +      }
544 +    }
545  
546 +
547      if (j_is_Dipole) {
548        mu_j = data2.dipole_moment;
549        uz_j = idat.eFrame2->getColumn(2);
# Line 580 | Line 604 | namespace OpenMD {
604            dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
605  
606          } else if (summationMethod_ == esm_REACTION_FIELD) {
607 <          rfVal =  *(idat.electroMult) * preRF_ *  *(idat.rij)  *  *(idat.rij) ;
607 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
608 >
609            vterm = preVal * ( riji + rfVal );            
610            dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
611 +          
612 +          // if this is an excluded pair, there are still indirect
613 +          // interactions via the reaction field we must worry about:
614  
615 +          if (idat.excluded) {
616 +            indirect_vpair += preVal * rfVal;
617 +            indirect_Pot += *(idat.sw) * preVal * rfVal;
618 +            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
619 +          }
620 +          
621          } else {
588 –          vterm = preVal * riji * erfcVal;            
622  
623 +          vterm = preVal * riji * erfcVal;          
624            dudr  = -  *(idat.sw)  * preVal * c2;
625  
626          }
627 <
628 <        *(idat.vpair) += vterm;
627 >
628 >        vpair += vterm;
629          epot +=  *(idat.sw)  * vterm;
630 <
597 <        dVdr += dudr * rhat;      
630 >        dVdr += dudr * rhat;                
631        }
632  
633        if (j_is_Dipole) {
# Line 607 | Line 640 | namespace OpenMD {
640            ri3 = ri2 * riji;
641      
642            vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
643 <          *(idat.vpair) += vterm;
643 >          vpair += vterm;
644            epot +=  *(idat.sw)  * vterm;
645  
646            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
647            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
648  
649 +          // Even if we excluded this pair from direct interactions,
650 +          // we still have the reaction-field-mediated charge-dipole
651 +          // interaction:
652 +
653 +          if (idat.excluded) {
654 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
655 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
656 +            indirect_dVdr += preSw * preRF2_ * uz_j;
657 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
658 +          }
659 +                      
660          } else {
661            // determine the inverse r used if we have split dipoles
662            if (j_is_SplitDipole) {
# Line 645 | Line 689 | namespace OpenMD {
689            // calculate the potential
690            pot_term =  scale * c2;
691            vterm = -pref * ct_j * pot_term;
692 <          *(idat.vpair) += vterm;
692 >          vpair += vterm;
693            epot +=  *(idat.sw)  * vterm;
694              
695            // calculate derivatives for forces and torques
# Line 692 | Line 736 | namespace OpenMD {
736                       qyy_j * (cy2*c3 - c2ri) +
737                       qzz_j * (cz2*c3 - c2ri) );
738          vterm = pref * pot_term;
739 <        *(idat.vpair) += vterm;
739 >        vpair += vterm;
740          epot +=  *(idat.sw)  * vterm;
741                  
742          // calculate derivatives for the forces and torques
# Line 720 | Line 764 | namespace OpenMD {
764            ri3 = ri2 * riji;
765  
766            vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
767 <          *(idat.vpair) += vterm;
767 >          vpair += vterm;
768            epot +=  *(idat.sw)  * vterm;
769            
770            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
771            
772            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
773 +
774 +          // Even if we excluded this pair from direct interactions,
775 +          // we still have the reaction-field-mediated charge-dipole
776 +          // interaction:
777 +
778 +          if (idat.excluded) {
779 +            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
780 +            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
781 +            indirect_dVdr += -preSw * preRF2_ * uz_i;
782 +            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
783 +          }
784              
785          } else {
786            
# Line 760 | Line 815 | namespace OpenMD {
815            // calculate the potential
816            pot_term = c2 * scale;
817            vterm = pref * ct_i * pot_term;
818 <          *(idat.vpair) += vterm;
818 >          vpair += vterm;
819            epot +=  *(idat.sw)  * vterm;
820  
821            // calculate derivatives for the forces and torques
# Line 783 | Line 838 | namespace OpenMD {
838  
839            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
840                             preRF2_ * ct_ij );
841 <          *(idat.vpair) += vterm;
841 >          vpair += vterm;
842            epot +=  *(idat.sw)  * vterm;
843              
844            a1 = 5.0 * ct_i * ct_j - ct_ij;
# Line 793 | Line 848 | namespace OpenMD {
848            duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
849            duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
850  
851 +          if (idat.excluded) {
852 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
853 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
854 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
855 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
856 +          }
857 +
858          } else {
859            
860            if (i_is_SplitDipole) {
# Line 842 | Line 904 | namespace OpenMD {
904            // calculate the potential
905            pot_term = (ct_ij * c2ri - ctidotj * c3);
906            vterm = pref * pot_term;
907 <          *(idat.vpair) += vterm;
907 >          vpair += vterm;
908            epot +=  *(idat.sw)  * vterm;
909  
910            // calculate derivatives for the forces and torques
# Line 894 | Line 956 | namespace OpenMD {
956                       qzz_i * (cz2 * c3 - c2ri) );
957          
958          vterm = pref * pot_term;
959 <        *(idat.vpair) += vterm;
959 >        vpair += vterm;
960          epot +=  *(idat.sw)  * vterm;
961  
962          // calculate the derivatives for the forces and torques
# Line 909 | Line 971 | namespace OpenMD {
971        }
972      }
973  
912 –    idat.pot[ELECTROSTATIC_FAMILY] += epot;
913 –    *(idat.f1) += dVdr;
974  
975 <    if (i_is_Dipole || i_is_Quadrupole)
976 <      *(idat.t1) -= cross(uz_i, duduz_i);
977 <    if (i_is_Quadrupole) {
978 <      *(idat.t1) -= cross(ux_i, dudux_i);
919 <      *(idat.t1) -= cross(uy_i, duduy_i);
920 <    }
921 <    
922 <    if (j_is_Dipole || j_is_Quadrupole)
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);
927 <    }
928 <
929 <    return;
930 <  }  
931 <
932 <  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
933 <
934 <    if (!initialized_) initialize();
935 <    
936 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
937 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
938 <    
939 <    // logicals
940 <
941 <    bool i_is_Charge = data1.is_Charge;
942 <    bool i_is_Dipole = data1.is_Dipole;
943 <
944 <    bool j_is_Charge = data2.is_Charge;
945 <    bool j_is_Dipole = data2.is_Dipole;
946 <
947 <    RealType q_i, q_j;
948 <    
949 <    // The skippedCharge computation is needed by the real-space cutoff methods
950 <    // (i.e. shifted force and shifted potential)
951 <
952 <    if (i_is_Charge) {
953 <      q_i = data1.charge;
954 <      *(idat.skippedCharge2) += q_i;
955 <    }
956 <
957 <    if (j_is_Charge) {
958 <      q_j = data2.charge;
959 <      *(idat.skippedCharge1) += q_j;
960 <    }
961 <
962 <    // the rest of this function should only be necessary for reaction field.
963 <
964 <    if (summationMethod_ == esm_REACTION_FIELD) {
965 <      RealType riji, ri2, ri3;
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:
975 >    if (!idat.excluded) {
976 >      *(idat.vpair) += vpair;
977 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
978 >      *(idat.f1) += dVdr;
979        
980 <      riji = 1.0 /  *(idat.rij) ;
981 <      rhat =  *(idat.d)  * riji;
982 <
983 <      if (i_is_Dipole) {
984 <        mu_i = data1.dipole_moment;
978 <        uz_i = idat.eFrame1->getColumn(2);      
979 <        ct_i = dot(uz_i, rhat);
980 <        duduz_i = V3Zero;
980 >      if (i_is_Dipole || i_is_Quadrupole)
981 >        *(idat.t1) -= cross(uz_i, duduz_i);
982 >      if (i_is_Quadrupole) {
983 >        *(idat.t1) -= cross(ux_i, dudux_i);
984 >        *(idat.t1) -= cross(uy_i, duduy_i);
985        }
986 <            
987 <      if (j_is_Dipole) {
988 <        mu_j = data2.dipole_moment;
989 <        uz_j = idat.eFrame2->getColumn(2);      
990 <        ct_j = dot(uz_j, rhat);
991 <        duduz_j = V3Zero;
986 >      
987 >      if (j_is_Dipole || j_is_Quadrupole)
988 >        *(idat.t2) -= cross(uz_j, duduz_j);
989 >      if (j_is_Quadrupole) {
990 >        *(idat.t2) -= cross(uz_j, dudux_j);
991 >        *(idat.t2) -= cross(uz_j, duduy_j);
992        }
993 <    
994 <      if (i_is_Charge) {
995 <        if (j_is_Charge) {
996 <          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
997 <          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
998 <          vterm = preVal * rfVal;
999 <          myPot +=  *(idat.sw)  * vterm;        
1000 <          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 =  *(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 =  *(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 <      }
993 >
994 >    } else {
995 >
996 >      // only accumulate the forces and torques resulting from the
997 >      // indirect reaction field terms.
998 >      *(idat.vpair) += indirect_vpair;
999 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1000 >      *(idat.f1) += indirect_dVdr;
1001        
1022 –      // accumulate the forces and torques resulting from the self term
1023 –      idat.pot[ELECTROSTATIC_FAMILY] += myPot;
1024 –      *(idat.f1) += dVdr;
1025 –      
1002        if (i_is_Dipole)
1003 <        *(idat.t1) -= cross(uz_i, duduz_i);
1003 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
1004        if (j_is_Dipole)
1005 <        *(idat.t2) -= cross(uz_j, duduz_j);
1005 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1006      }
1007 <  }
1007 >
1008 >
1009 >    return;
1010 >  }  
1011      
1012    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1013      RealType mu1, preVal, chg1, self;
1014      
1015      if (!initialized_) initialize();
1016 <    
1016 >
1017      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1018    
1019      // logicals
1041 –
1020      bool i_is_Charge = data.is_Charge;
1021      bool i_is_Dipole = data.is_Dipole;
1022  
# Line 1046 | Line 1024 | namespace OpenMD {
1024        if (i_is_Dipole) {
1025          mu1 = data.dipole_moment;          
1026          preVal = pre22_ * preRF2_ * mu1 * mu1;
1027 <        sdat.pot[2] -= 0.5 * preVal;
1027 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1028          
1029          // The self-correction term adds into the reaction field vector
1030          Vector3d uz_i = sdat.eFrame->getColumn(2);
# Line 1063 | Line 1041 | namespace OpenMD {
1041          } else {        
1042            self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1043          }
1044 <        sdat.pot[ELECTROSTATIC_FAMILY] += self;
1044 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1045        }
1046      }
1047    }

Diff Legend

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