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 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC vs.
Revision 1601 by gezelter, Thu Aug 4 20:04:35 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), info_(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 <
63 >    
64      Globals* simParams_ = info_->getSimParams();
65  
66      summationMap_["HARD"]               = esm_HARD;
# 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 445 | Line 446 | namespace OpenMD {
446      RealType ct_i, ct_j, ct_ij, a1;
447      RealType riji, ri, ri2, ri3, ri4;
448      RealType pref, vterm, epot, dudr;
449 +    RealType vpair(0.0);
450      RealType scale, sc2;
451      RealType pot_term, preVal, rfVal;
452      RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
453      RealType preSw, preSwSc;
454      RealType c1, c2, c3, c4;
455 <    RealType erfcVal, derfcVal;
455 >    RealType erfcVal(1.0), derfcVal(0.0);
456      RealType BigR;
457  
458      Vector3d Q_i, Q_j;
# Line 460 | Line 462 | namespace OpenMD {
462      Vector3d dudux_j, duduy_j, duduz_j;
463      Vector3d rhatdot2, rhatc4;
464      Vector3d dVdr;
465 +
466 +    // variables for indirect (reaction field) interactions for excluded pairs:
467 +    RealType indirect_Pot(0.0);
468 +    RealType indirect_vpair(0.0);
469 +    Vector3d indirect_dVdr(V3Zero);
470 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
471  
472      pair<RealType, RealType> res;
473      
# Line 485 | Line 493 | namespace OpenMD {
493      bool j_is_SplitDipole = data2.is_SplitDipole;
494      bool j_is_Quadrupole = data2.is_Quadrupole;
495      
496 <    if (i_is_Charge)
496 >    if (i_is_Charge) {
497        q_i = data1.charge;
498 +      if (idat.excluded) {
499 +        *(idat.skippedCharge2) += q_i;
500 +      }
501 +    }
502  
503      if (i_is_Dipole) {
504        mu_i = data1.dipole_moment;
# Line 519 | Line 531 | namespace OpenMD {
531        duduz_i = V3Zero;
532      }
533  
534 <    if (j_is_Charge)
534 >    if (j_is_Charge) {
535        q_j = data2.charge;
536 +      if (idat.excluded) {
537 +        *(idat.skippedCharge1) += q_j;
538 +      }
539 +    }
540  
541 +
542      if (j_is_Dipole) {
543        mu_j = data2.dipole_moment;
544        uz_j = idat.eFrame2->getColumn(2);
# Line 582 | Line 599 | namespace OpenMD {
599            dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
600  
601          } else if (summationMethod_ == esm_REACTION_FIELD) {
602 <          rfVal =  *(idat.electroMult) * preRF_ *  *(idat.rij)  *  *(idat.rij) ;
602 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
603 >
604            vterm = preVal * ( riji + rfVal );            
605            dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
606 +          
607 +          // if this is an excluded pair, there are still indirect
608 +          // interactions via the reaction field we must worry about:
609  
610 +          if (idat.excluded) {
611 +            indirect_vpair += preVal * rfVal;
612 +            indirect_Pot += *(idat.sw) * preVal * rfVal;
613 +            indirect_dVdr += *(idat.sw)  * preVal * 2.0 * rfVal  * riji * rhat;
614 +          }
615 +          
616          } else {
590 –          vterm = preVal * riji * erfcVal;            
617  
618 +          vterm = preVal * riji * erfcVal;          
619            dudr  = -  *(idat.sw)  * preVal * c2;
620  
621          }
595 –
596 –        *(idat.vpair) += vterm;
597 –        epot +=  *(idat.sw)  * vterm;
622  
623 <        dVdr += dudr * rhat;      
623 >        vpair += vterm;
624 >        epot +=  *(idat.sw)  * vterm;
625 >        dVdr += dudr * rhat;                
626        }
627  
628        if (j_is_Dipole) {
# Line 609 | Line 635 | namespace OpenMD {
635            ri3 = ri2 * riji;
636      
637            vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
638 <          *(idat.vpair) += vterm;
638 >          vpair += vterm;
639            epot +=  *(idat.sw)  * vterm;
640  
641            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
642            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
643  
644 +          // Even if we excluded this pair from direct interactions,
645 +          // we still have the reaction-field-mediated charge-dipole
646 +          // interaction:
647 +
648 +          if (idat.excluded) {
649 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
650 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
651 +            indirect_dVdr += preSw * preRF2_ * uz_j;
652 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
653 +          }
654 +                      
655          } else {
656            // determine the inverse r used if we have split dipoles
657            if (j_is_SplitDipole) {
# Line 647 | Line 684 | namespace OpenMD {
684            // calculate the potential
685            pot_term =  scale * c2;
686            vterm = -pref * ct_j * pot_term;
687 <          *(idat.vpair) += vterm;
687 >          vpair += vterm;
688            epot +=  *(idat.sw)  * vterm;
689              
690            // calculate derivatives for forces and torques
# Line 694 | Line 731 | namespace OpenMD {
731                       qyy_j * (cy2*c3 - c2ri) +
732                       qzz_j * (cz2*c3 - c2ri) );
733          vterm = pref * pot_term;
734 <        *(idat.vpair) += vterm;
734 >        vpair += vterm;
735          epot +=  *(idat.sw)  * vterm;
736                  
737          // calculate derivatives for the forces and torques
# Line 722 | Line 759 | namespace OpenMD {
759            ri3 = ri2 * riji;
760  
761            vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
762 <          *(idat.vpair) += vterm;
762 >          vpair += vterm;
763            epot +=  *(idat.sw)  * vterm;
764            
765            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
766            
767            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
768 +
769 +          // Even if we excluded this pair from direct interactions,
770 +          // we still have the reaction-field-mediated charge-dipole
771 +          // interaction:
772 +
773 +          if (idat.excluded) {
774 +            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
775 +            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
776 +            indirect_dVdr += -preSw * preRF2_ * uz_i;
777 +            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
778 +          }
779              
780          } else {
781            
# Line 762 | Line 810 | namespace OpenMD {
810            // calculate the potential
811            pot_term = c2 * scale;
812            vterm = pref * ct_i * pot_term;
813 <          *(idat.vpair) += vterm;
813 >          vpair += vterm;
814            epot +=  *(idat.sw)  * vterm;
815  
816            // calculate derivatives for the forces and torques
# Line 785 | Line 833 | namespace OpenMD {
833  
834            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
835                             preRF2_ * ct_ij );
836 <          *(idat.vpair) += vterm;
836 >          vpair += vterm;
837            epot +=  *(idat.sw)  * vterm;
838              
839            a1 = 5.0 * ct_i * ct_j - ct_ij;
# Line 794 | Line 842 | namespace OpenMD {
842  
843            duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
844            duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
845 +
846 +          if (idat.excluded) {
847 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
848 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
849 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
850 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
851 +          }
852  
853          } else {
854            
# Line 844 | Line 899 | namespace OpenMD {
899            // calculate the potential
900            pot_term = (ct_ij * c2ri - ctidotj * c3);
901            vterm = pref * pot_term;
902 <          *(idat.vpair) += vterm;
902 >          vpair += vterm;
903            epot +=  *(idat.sw)  * vterm;
904  
905            // calculate derivatives for the forces and torques
# Line 896 | Line 951 | namespace OpenMD {
951                       qzz_i * (cz2 * c3 - c2ri) );
952          
953          vterm = pref * pot_term;
954 <        *(idat.vpair) += vterm;
954 >        vpair += vterm;
955          epot +=  *(idat.sw)  * vterm;
956  
957          // calculate the derivatives for the forces and torques
# Line 911 | Line 966 | namespace OpenMD {
966        }
967      }
968  
914 –    (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
915 –    *(idat.f1) += dVdr;
969  
970 <    if (i_is_Dipole || i_is_Quadrupole)
971 <      *(idat.t1) -= cross(uz_i, duduz_i);
972 <    if (i_is_Quadrupole) {
973 <      *(idat.t1) -= cross(ux_i, dudux_i);
974 <      *(idat.t1) -= cross(uy_i, duduy_i);
975 <    }
976 <    
977 <    if (j_is_Dipole || j_is_Quadrupole)
978 <      *(idat.t2) -= cross(uz_j, duduz_j);
979 <    if (j_is_Quadrupole) {
980 <      *(idat.t2) -= cross(uz_j, dudux_j);
981 <      *(idat.t2) -= cross(uz_j, duduy_j);
982 <    }
970 >    if (!idat.excluded) {
971 >      *(idat.vpair) += vpair;
972 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
973 >      *(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);
979 >        *(idat.t1) -= cross(uy_i, duduy_i);
980 >      }
981 >      
982 >      if (j_is_Dipole || j_is_Quadrupole)
983 >        *(idat.t2) -= cross(uz_j, duduz_j);
984 >      if (j_is_Quadrupole) {
985 >        *(idat.t2) -= cross(uz_j, dudux_j);
986 >        *(idat.t2) -= cross(uz_j, duduy_j);
987 >      }
988  
989 <    return;
932 <  }  
989 >    } else {
990  
991 <  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
992 <
993 <    if (!initialized_) initialize();
994 <    
995 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
939 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
940 <    
941 <    // logicals
942 <
943 <    bool i_is_Charge = data1.is_Charge;
944 <    bool i_is_Dipole = data1.is_Dipole;
945 <
946 <    bool j_is_Charge = data2.is_Charge;
947 <    bool j_is_Dipole = data2.is_Dipole;
948 <
949 <    RealType q_i, q_j;
950 <    
951 <    // The skippedCharge computation is needed by the real-space cutoff methods
952 <    // (i.e. shifted force and shifted potential)
953 <
954 <    if (i_is_Charge) {
955 <      q_i = data1.charge;
956 <      *(idat.skippedCharge2) += q_i;
957 <    }
958 <
959 <    if (j_is_Charge) {
960 <      q_j = data2.charge;
961 <      *(idat.skippedCharge1) += q_j;
962 <    }
963 <
964 <    // the rest of this function should only be necessary for reaction field.
965 <
966 <    if (summationMethod_ == esm_REACTION_FIELD) {
967 <      RealType riji, ri2, ri3;
968 <      RealType mu_i, ct_i;
969 <      RealType mu_j, ct_j;
970 <      RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
971 <      Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
972 <
973 <      // some variables we'll need independent of electrostatic type:
991 >      // only accumulate the forces and torques resulting from the
992 >      // indirect reaction field terms.
993 >      *(idat.vpair) += indirect_vpair;
994 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
995 >      *(idat.f1) += indirect_dVdr;
996        
975 –      riji = 1.0 /  *(idat.rij) ;
976 –      rhat =  *(idat.d)  * riji;
977 –
978 –      if (i_is_Dipole) {
979 –        mu_i = data1.dipole_moment;
980 –        uz_i = idat.eFrame1->getColumn(2);      
981 –        ct_i = dot(uz_i, rhat);
982 –        duduz_i = V3Zero;
983 –      }
984 –            
985 –      if (j_is_Dipole) {
986 –        mu_j = data2.dipole_moment;
987 –        uz_j = idat.eFrame2->getColumn(2);      
988 –        ct_j = dot(uz_j, rhat);
989 –        duduz_j = V3Zero;
990 –      }
991 –    
992 –      if (i_is_Charge) {
993 –        if (j_is_Charge) {
994 –          preVal =  *(idat.electroMult) * pre11_ * q_i * q_j;
995 –          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij) ;
996 –          vterm = preVal * rfVal;
997 –          myPot +=  *(idat.sw)  * vterm;        
998 –          dudr  =  *(idat.sw)  * preVal * 2.0 * rfVal * riji;        
999 –          dVdr += dudr * rhat;
1000 –        }
1001 –        
1002 –        if (j_is_Dipole) {
1003 –          ri2 = riji * riji;
1004 –          ri3 = ri2 * riji;        
1005 –          pref =  *(idat.electroMult) * pre12_ * q_i * mu_j;
1006 –          vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
1007 –          myPot +=  *(idat.sw)  * vterm;        
1008 –          dVdr += - *(idat.sw)  * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1009 –          duduz_j += - *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij) );
1010 –        }
1011 –      }
1012 –      if (i_is_Dipole) {
1013 –        if (j_is_Charge) {
1014 –          ri2 = riji * riji;
1015 –          ri3 = ri2 * riji;        
1016 –          pref =  *(idat.electroMult) * pre12_ * q_j * mu_i;
1017 –          vterm = - pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
1018 –          myPot +=  *(idat.sw)  * vterm;        
1019 –          dVdr +=  *(idat.sw)  * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1020 –          duduz_i +=  *(idat.sw)  * pref * rhat * (ri2 - preRF2_ *  *(idat.rij));
1021 –        }
1022 –      }
1023 –      
1024 –      // accumulate the forces and torques resulting from the self term
1025 –      (*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1026 –      *(idat.f1) += dVdr;
1027 –      
997        if (i_is_Dipole)
998 <        *(idat.t1) -= cross(uz_i, duduz_i);
998 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
999        if (j_is_Dipole)
1000 <        *(idat.t2) -= cross(uz_j, duduz_j);
1000 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1001      }
1002 <  }
1002 >
1003 >
1004 >    return;
1005 >  }  
1006      
1007    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1008      RealType mu1, preVal, chg1, self;
1009      
1010      if (!initialized_) initialize();
1011 <    
1011 >
1012      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1013    
1014      // logicals
1043 –
1015      bool i_is_Charge = data.is_Charge;
1016      bool i_is_Dipole = data.is_Dipole;
1017  
# Line 1048 | Line 1019 | namespace OpenMD {
1019        if (i_is_Dipole) {
1020          mu1 = data.dipole_moment;          
1021          preVal = pre22_ * preRF2_ * mu1 * mu1;
1022 <        sdat.pot[2] -= 0.5 * preVal;
1022 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1023          
1024          // The self-correction term adds into the reaction field vector
1025          Vector3d uz_i = sdat.eFrame->getColumn(2);
# Line 1065 | Line 1036 | namespace OpenMD {
1036          } else {        
1037            self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1038          }
1039 <        sdat.pot[ELECTROSTATIC_FAMILY] += self;
1039 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1040        }
1041      }
1042    }

Diff Legend

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