ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/devel_omp/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.
branches/devel_omp/src/nonbonded/Electrostatic.cpp (file contents), Revision 1614 by mciznick, Tue Aug 23 20:55:51 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 407 | Line 408 | namespace OpenMD {
408      return;
409    }
410    
411 <  void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
412 <                                                    RealType theRSW ) {
412 <    cutoffRadius_ = theECR;
411 >  void Electrostatic::setCutoffRadius( RealType rCut ) {
412 >    cutoffRadius_ = rCut;
413      rrf_ = cutoffRadius_;
414 –    rt_ = theRSW;
414      haveCutoffRadius_ = true;
415 +  }
416 +
417 +  void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
418 +    rt_ = rSwitch;
419    }
420    void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
421      summationMethod_ = esm;
# Line 429 | Line 432 | namespace OpenMD {
432      haveDielectric_ = true;
433    }
434  
435 +  void Electrostatic::initForce() {
436 +        if (!initialized_) initialize();
437 +  }
438 +
439    void Electrostatic::calcForce(InteractionData &idat) {
440  
441      // utility variables.  Should clean these up and use the Vector3d and
# Line 443 | Line 450 | namespace OpenMD {
450      RealType ct_i, ct_j, ct_ij, a1;
451      RealType riji, ri, ri2, ri3, ri4;
452      RealType pref, vterm, epot, dudr;
453 +    RealType vpair(0.0);
454      RealType scale, sc2;
455      RealType pot_term, preVal, rfVal;
456      RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
457      RealType preSw, preSwSc;
458      RealType c1, c2, c3, c4;
459 <    RealType erfcVal, derfcVal;
459 >    RealType erfcVal(1.0), derfcVal(0.0);
460      RealType BigR;
461  
462      Vector3d Q_i, Q_j;
# Line 459 | Line 467 | namespace OpenMD {
467      Vector3d rhatdot2, rhatc4;
468      Vector3d dVdr;
469  
470 +    // variables for indirect (reaction field) interactions for excluded pairs:
471 +    RealType indirect_Pot(0.0);
472 +    RealType indirect_vpair(0.0);
473 +    Vector3d indirect_dVdr(V3Zero);
474 +    Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
475 +
476      pair<RealType, RealType> res;
477      
478      if (!initialized_) initialize();
# Line 483 | Line 497 | namespace OpenMD {
497      bool j_is_SplitDipole = data2.is_SplitDipole;
498      bool j_is_Quadrupole = data2.is_Quadrupole;
499      
500 <    if (i_is_Charge)
500 >    if (i_is_Charge) {
501        q_i = data1.charge;
502 +      if (idat.excluded) {
503 +        *(idat.skippedCharge2) += q_i;
504 +      }
505 +    }
506  
507      if (i_is_Dipole) {
508        mu_i = data1.dipole_moment;
# Line 517 | Line 535 | namespace OpenMD {
535        duduz_i = V3Zero;
536      }
537  
538 <    if (j_is_Charge)
538 >    if (j_is_Charge) {
539        q_j = data2.charge;
540 +      if (idat.excluded) {
541 +        *(idat.skippedCharge1) += q_j;
542 +      }
543 +    }
544  
545 +
546      if (j_is_Dipole) {
547        mu_j = data2.dipole_moment;
548        uz_j = idat.eFrame2->getColumn(2);
# Line 580 | Line 603 | namespace OpenMD {
603            dudr  =  *(idat.sw)  * preVal * (c2c_ - c2);
604  
605          } else if (summationMethod_ == esm_REACTION_FIELD) {
606 <          rfVal =  *(idat.electroMult) * preRF_ *  *(idat.rij)  *  *(idat.rij) ;
606 >          rfVal = preRF_ *  *(idat.rij)  *  *(idat.rij);
607 >
608            vterm = preVal * ( riji + rfVal );            
609            dudr  =  *(idat.sw)  * preVal * ( 2.0 * rfVal - riji ) * riji;
610 +          
611 +          // if this is an excluded pair, there are still indirect
612 +          // interactions via the reaction field we must worry about:
613  
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;
618 +          }
619 +          
620          } else {
588 –          vterm = preVal * riji * erfcVal;            
621  
622 +          vterm = preVal * riji * erfcVal;          
623            dudr  = -  *(idat.sw)  * preVal * c2;
624  
625          }
593 –
594 –        *(idat.vpair) += vterm;
595 –        epot +=  *(idat.sw)  * vterm;
626  
627 <        dVdr += dudr * rhat;      
627 >        vpair += vterm;
628 >        epot +=  *(idat.sw)  * vterm;
629 >        dVdr += dudr * rhat;                
630        }
631  
632        if (j_is_Dipole) {
# Line 607 | Line 639 | namespace OpenMD {
639            ri3 = ri2 * riji;
640      
641            vterm = - pref * ct_j * ( ri2 - preRF2_ *  *(idat.rij)  );
642 <          *(idat.vpair) += vterm;
642 >          vpair += vterm;
643            epot +=  *(idat.sw)  * vterm;
644  
645            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
646            duduz_j += -preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );  
647  
648 +          // Even if we excluded this pair from direct interactions,
649 +          // we still have the reaction-field-mediated charge-dipole
650 +          // interaction:
651 +
652 +          if (idat.excluded) {
653 +            indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
654 +            indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
655 +            indirect_dVdr += preSw * preRF2_ * uz_j;
656 +            indirect_duduz_j += preSw * rhat * preRF2_ *  *(idat.rij);
657 +          }
658 +                      
659          } else {
660            // determine the inverse r used if we have split dipoles
661            if (j_is_SplitDipole) {
# Line 645 | Line 688 | namespace OpenMD {
688            // calculate the potential
689            pot_term =  scale * c2;
690            vterm = -pref * ct_j * pot_term;
691 <          *(idat.vpair) += vterm;
691 >          vpair += vterm;
692            epot +=  *(idat.sw)  * vterm;
693              
694            // calculate derivatives for forces and torques
# Line 692 | Line 735 | namespace OpenMD {
735                       qyy_j * (cy2*c3 - c2ri) +
736                       qzz_j * (cz2*c3 - c2ri) );
737          vterm = pref * pot_term;
738 <        *(idat.vpair) += vterm;
738 >        vpair += vterm;
739          epot +=  *(idat.sw)  * vterm;
740                  
741          // calculate derivatives for the forces and torques
# Line 720 | Line 763 | namespace OpenMD {
763            ri3 = ri2 * riji;
764  
765            vterm = pref * ct_i * ( ri2 - preRF2_ *  *(idat.rij)  );
766 <          *(idat.vpair) += vterm;
766 >          vpair += vterm;
767            epot +=  *(idat.sw)  * vterm;
768            
769            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
770            
771            duduz_i += preSw * rhat * (ri2 - preRF2_ *  *(idat.rij) );
772 +
773 +          // Even if we excluded this pair from direct interactions,
774 +          // we still have the reaction-field-mediated charge-dipole
775 +          // interaction:
776 +
777 +          if (idat.excluded) {
778 +            indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
779 +            indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
780 +            indirect_dVdr += -preSw * preRF2_ * uz_i;
781 +            indirect_duduz_i += -preSw * rhat * preRF2_ *  *(idat.rij);
782 +          }
783              
784          } else {
785            
# Line 760 | Line 814 | namespace OpenMD {
814            // calculate the potential
815            pot_term = c2 * scale;
816            vterm = pref * ct_i * pot_term;
817 <          *(idat.vpair) += vterm;
817 >          vpair += vterm;
818            epot +=  *(idat.sw)  * vterm;
819  
820            // calculate derivatives for the forces and torques
# Line 783 | Line 837 | namespace OpenMD {
837  
838            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
839                             preRF2_ * ct_ij );
840 <          *(idat.vpair) += vterm;
840 >          vpair += vterm;
841            epot +=  *(idat.sw)  * vterm;
842              
843            a1 = 5.0 * ct_i * ct_j - ct_ij;
# Line 793 | Line 847 | namespace OpenMD {
847            duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
848            duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
849  
850 +          if (idat.excluded) {
851 +            indirect_vpair +=  - pref * preRF2_ * ct_ij;
852 +            indirect_Pot +=    - preSw * preRF2_ * ct_ij;
853 +            indirect_duduz_i += -preSw * preRF2_ * uz_j;
854 +            indirect_duduz_j += -preSw * preRF2_ * uz_i;
855 +          }
856 +
857          } else {
858            
859            if (i_is_SplitDipole) {
# Line 842 | Line 903 | namespace OpenMD {
903            // calculate the potential
904            pot_term = (ct_ij * c2ri - ctidotj * c3);
905            vterm = pref * pot_term;
906 <          *(idat.vpair) += vterm;
906 >          vpair += vterm;
907            epot +=  *(idat.sw)  * vterm;
908  
909            // calculate derivatives for the forces and torques
# Line 894 | Line 955 | namespace OpenMD {
955                       qzz_i * (cz2 * c3 - c2ri) );
956          
957          vterm = pref * pot_term;
958 <        *(idat.vpair) += vterm;
958 >        vpair += vterm;
959          epot +=  *(idat.sw)  * vterm;
960  
961          // calculate the derivatives for the forces and torques
# Line 909 | Line 970 | namespace OpenMD {
970        }
971      }
972  
912 –    idat.pot[ELECTROSTATIC_FAMILY] += epot;
913 –    *(idat.f1) += dVdr;
973  
974 <    if (i_is_Dipole || i_is_Quadrupole)
975 <      *(idat.t1) -= cross(uz_i, duduz_i);
976 <    if (i_is_Quadrupole) {
977 <      *(idat.t1) -= cross(ux_i, dudux_i);
978 <      *(idat.t1) -= cross(uy_i, duduy_i);
979 <    }
980 <    
981 <    if (j_is_Dipole || j_is_Quadrupole)
982 <      *(idat.t2) -= cross(uz_j, duduz_j);
983 <    if (j_is_Quadrupole) {
984 <      *(idat.t2) -= cross(uz_j, dudux_j);
985 <      *(idat.t2) -= cross(uz_j, duduy_j);
974 >    if (!idat.excluded) {
975 >      *(idat.vpair) += vpair;
976 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
977 >      *(idat.f1) += dVdr;
978 >      
979 >      if (i_is_Dipole || i_is_Quadrupole)
980 >        *(idat.t1) -= cross(uz_i, duduz_i);
981 >      if (i_is_Quadrupole) {
982 >        *(idat.t1) -= cross(ux_i, dudux_i);
983 >        *(idat.t1) -= cross(uy_i, duduy_i);
984 >      }
985 >      
986 >      if (j_is_Dipole || j_is_Quadrupole)
987 >        *(idat.t2) -= cross(uz_j, duduz_j);
988 >      if (j_is_Quadrupole) {
989 >        *(idat.t2) -= cross(uz_j, dudux_j);
990 >        *(idat.t2) -= cross(uz_j, duduy_j);
991 >      }
992 >
993 >    } else {
994 >
995 >      // only accumulate the forces and torques resulting from the
996 >      // indirect reaction field terms.
997 >      *(idat.vpair) += indirect_vpair;
998 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
999 >      *(idat.f1) += indirect_dVdr;
1000 >      
1001 >      if (i_is_Dipole)
1002 >        *(idat.t1) -= cross(uz_i, indirect_duduz_i);
1003 >      if (j_is_Dipole)
1004 >        *(idat.t2) -= cross(uz_j, indirect_duduz_j);
1005      }
1006  
1007 +
1008      return;
1009    }  
1010  
# Line 945 | Line 1024 | namespace OpenMD {
1024      bool j_is_Dipole = data2.is_Dipole;
1025  
1026      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)
1027  
1028 +    // The skippedCharge computation is needed by the real-space
1029 +    // cutoff methods (i.e. shifted force and shifted potential)
1030 +      
1031      if (i_is_Charge) {
1032        q_i = data1.charge;
1033        *(idat.skippedCharge2) += q_i;
1034      }
1035 <
1035 >    
1036      if (j_is_Charge) {
1037        q_j = data2.charge;
1038        *(idat.skippedCharge1) += q_j;
# Line 1020 | Line 1099 | namespace OpenMD {
1099        }
1100        
1101        // accumulate the forces and torques resulting from the self term
1102 <      idat.pot[ELECTROSTATIC_FAMILY] += myPot;
1102 >      (*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1103        *(idat.f1) += dVdr;
1104        
1105        if (i_is_Dipole)
# Line 1034 | Line 1113 | namespace OpenMD {
1113      RealType mu1, preVal, chg1, self;
1114      
1115      if (!initialized_) initialize();
1116 <    
1116 >
1117      ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1118    
1119      // logicals
1041 –
1120      bool i_is_Charge = data.is_Charge;
1121      bool i_is_Dipole = data.is_Dipole;
1122  
# Line 1046 | Line 1124 | namespace OpenMD {
1124        if (i_is_Dipole) {
1125          mu1 = data.dipole_moment;          
1126          preVal = pre22_ * preRF2_ * mu1 * mu1;
1127 <        sdat.pot[2] -= 0.5 * preVal;
1127 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal;
1128          
1129          // The self-correction term adds into the reaction field vector
1130          Vector3d uz_i = sdat.eFrame->getColumn(2);
# Line 1063 | Line 1141 | namespace OpenMD {
1141          } else {        
1142            self = - 0.5 * rcuti_ * chg1 * (chg1 +  *(sdat.skippedCharge)) * pre11_;
1143          }
1144 <        sdat.pot[ELECTROSTATIC_FAMILY] += self;
1144 >        (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1145        }
1146      }
1147    }

Diff Legend

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