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 1528 by gezelter, Fri Dec 17 20:11:05 2010 UTC vs.
Revision 1545 by gezelter, Fri Apr 8 21:25:19 2011 UTC

# Line 47 | Line 47
47   #include "utils/simError.h"
48   #include "types/NonBondedInteractionType.hpp"
49   #include "types/DirectionalAtomType.hpp"
50 + #include "io/Globals.hpp"
51  
51 –
52   namespace OpenMD {
53    
54    Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55 <                                  forceField_(NULL), info_(NULL) {}
55 >                                  forceField_(NULL) {}
56    
57    void Electrostatic::initialize() {
58  
59 +    Globals* simParams_;
60 +
61      summationMap_["HARD"]               = esm_HARD;
62      summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
63      summationMap_["SHIFTED_POTENTIAL"]  = esm_SHIFTED_POTENTIAL;
# Line 111 | Line 113 | namespace OpenMD {
113        } else {
114          // throw error
115          sprintf( painCave.errMsg,
116 <                 "SimInfo error: Unknown electrostaticSummationMethod.\n"
116 >                 "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
117                   "\t(Input file specified %s .)\n"
118                   "\telectrostaticSummationMethod must be one of: \"none\",\n"
119                   "\t\"shifted_potential\", \"shifted_force\", or \n"
# Line 427 | Line 429 | namespace OpenMD {
429      haveDielectric_ = true;
430    }
431  
432 <  void Electrostatic::calcForce(InteractionData idat) {
432 >  void Electrostatic::calcForce(InteractionData &idat) {
433  
434      // utility variables.  Should clean these up and use the Vector3d and
435      // Mat3x3d to replace as many as we can in future versions:
# Line 461 | Line 463 | namespace OpenMD {
463      
464      if (!initialized_) initialize();
465      
466 <    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
467 <    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
466 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
467 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
468      
469      // some variables we'll need independent of electrostatic type:
470  
# Line 589 | Line 591 | namespace OpenMD {
591  
592          }
593  
594 <        idat.vpair += vterm;
594 >        idat.vpair[2] += vterm;
595          epot += idat.sw * vterm;
596  
597          dVdr += dudr * rhat;      
# Line 605 | Line 607 | namespace OpenMD {
607            ri3 = ri2 * riji;
608      
609            vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
610 <          idat.vpair += vterm;
610 >          idat.vpair[2] += vterm;
611            epot += idat.sw * vterm;
612  
613            dVdr +=  -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
# Line 643 | Line 645 | namespace OpenMD {
645            // calculate the potential
646            pot_term =  scale * c2;
647            vterm = -pref * ct_j * pot_term;
648 <          idat.vpair += vterm;
648 >          idat.vpair[2] += vterm;
649            epot += idat.sw * vterm;
650              
651            // calculate derivatives for forces and torques
# Line 690 | Line 692 | namespace OpenMD {
692                       qyy_j * (cy2*c3 - c2ri) +
693                       qzz_j * (cz2*c3 - c2ri) );
694          vterm = pref * pot_term;
695 <        idat.vpair += vterm;
695 >        idat.vpair[2] += vterm;
696          epot += idat.sw * vterm;
697                  
698          // calculate derivatives for the forces and torques
# Line 718 | Line 720 | namespace OpenMD {
720            ri3 = ri2 * riji;
721  
722            vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
723 <          idat.vpair += vterm;
723 >          idat.vpair[2] += vterm;
724            epot += idat.sw * vterm;
725            
726            dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
# Line 758 | Line 760 | namespace OpenMD {
760            // calculate the potential
761            pot_term = c2 * scale;
762            vterm = pref * ct_i * pot_term;
763 <          idat.vpair += vterm;
763 >          idat.vpair[2] += vterm;
764            epot += idat.sw * vterm;
765  
766            // calculate derivatives for the forces and torques
# Line 781 | Line 783 | namespace OpenMD {
783  
784            vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
785                             preRF2_ * ct_ij );
786 <          idat.vpair += vterm;
786 >          idat.vpair[2] += vterm;
787            epot += idat.sw * vterm;
788              
789            a1 = 5.0 * ct_i * ct_j - ct_ij;
# Line 840 | Line 842 | namespace OpenMD {
842            // calculate the potential
843            pot_term = (ct_ij * c2ri - ctidotj * c3);
844            vterm = pref * pot_term;
845 <          idat.vpair += vterm;
845 >          idat.vpair[2] += vterm;
846            epot += idat.sw * vterm;
847  
848            // calculate derivatives for the forces and torques
# Line 892 | Line 894 | namespace OpenMD {
894                       qzz_i * (cz2 * c3 - c2ri) );
895          
896          vterm = pref * pot_term;
897 <        idat.vpair += vterm;
897 >        idat.vpair[2] += vterm;
898          epot += idat.sw * vterm;
899  
900          // calculate the derivatives for the forces and torques
# Line 907 | Line 909 | namespace OpenMD {
909        }
910      }
911  
912 <    idat.pot += epot;
912 >    idat.pot[2] += epot;
913      idat.f1 += dVdr;
914  
915      if (i_is_Dipole || i_is_Quadrupole)
# Line 927 | Line 929 | namespace OpenMD {
929      return;
930    }  
931  
932 <  void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
932 >  void Electrostatic::calcSkipCorrection(InteractionData &idat) {
933  
934      if (!initialized_) initialize();
935      
936 <    ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
937 <    ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
936 >    ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
937 >    ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
938      
939      // logicals
940  
# Line 949 | Line 951 | namespace OpenMD {
951  
952      if (i_is_Charge) {
953        q_i = data1.charge;
954 <      skdat.skippedCharge2 += q_i;
954 >      idat.skippedCharge2 += q_i;
955      }
956  
957      if (j_is_Charge) {
958        q_j = data2.charge;
959 <      skdat.skippedCharge1 += q_j;
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 q_i, mu_i, ct_i;
967 <      RealType q_j, mu_j, ct_j;
968 <      RealType preVal, rfVal, vterm, dudr, pref, myPot;
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:
972        
973 <      riji = 1.0 / skdat.rij;
974 <      rhat = skdat.d  * riji;
973 >      riji = 1.0 / idat.rij;
974 >      rhat = idat.d  * riji;
975  
976        if (i_is_Dipole) {
977          mu_i = data1.dipole_moment;
978 <        uz_i = skdat.eFrame1.getColumn(2);      
978 >        uz_i = idat.eFrame1.getColumn(2);      
979          ct_i = dot(uz_i, rhat);
980          duduz_i = V3Zero;
981        }
982              
983        if (j_is_Dipole) {
984          mu_j = data2.dipole_moment;
985 <        uz_j = skdat.eFrame2.getColumn(2);      
985 >        uz_j = idat.eFrame2.getColumn(2);      
986          ct_j = dot(uz_j, rhat);
987          duduz_j = V3Zero;
988        }
989      
990        if (i_is_Charge) {
991          if (j_is_Charge) {
992 <          preVal = skdat.electroMult * pre11_ * q_i * q_j;
993 <          rfVal = preRF_ * skdat.rij * skdat.rij;
992 >          preVal = idat.electroMult * pre11_ * q_i * q_j;
993 >          rfVal = preRF_ * idat.rij * idat.rij;
994            vterm = preVal * rfVal;
995 <          myPot += skdat.sw * vterm;        
996 <          dudr  = skdat.sw * preVal * 2.0 * rfVal * riji;        
995 >          myPot += idat.sw * vterm;        
996 >          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 = skdat.electroMult * pre12_ * q_i * mu_j;
1004 <          vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1005 <          myPot += skdat.sw * vterm;        
1006 <          dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1007 <          duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
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 = skdat.electroMult * pre12_ * q_j * mu_i;
1015 <          vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1016 <          myPot += skdat.sw * vterm;        
1017 <          dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);      
1018 <          duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
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        }
1021        
1022        // accumulate the forces and torques resulting from the self term
1023 <      skdat.pot += myPot;
1024 <      skdat.f1 += dVdr;
1023 >      idat.pot[2] += myPot;
1024 >      idat.f1 += dVdr;
1025        
1026        if (i_is_Dipole)
1027 <        skdat.t1 -= cross(uz_i, duduz_i);
1027 >        idat.t1 -= cross(uz_i, duduz_i);
1028        if (j_is_Dipole)
1029 <        skdat.t2 -= cross(uz_j, duduz_j);
1029 >        idat.t2 -= cross(uz_j, duduz_j);
1030      }
1031    }
1032      
1033 <  void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1033 >  void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1034      RealType mu1, preVal, chg1, self;
1035      
1036      if (!initialized_) initialize();
1037      
1038 <    ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1038 >    ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1039    
1040      // logicals
1041  
# Line 1044 | Line 1046 | namespace OpenMD {
1046        if (i_is_Dipole) {
1047          mu1 = data.dipole_moment;          
1048          preVal = pre22_ * preRF2_ * mu1 * mu1;
1049 <        scdat.pot -= 0.5 * preVal;
1049 >        sdat.pot[2] -= 0.5 * preVal;
1050          
1051          // The self-correction term adds into the reaction field vector
1052 <        Vector3d uz_i = scdat.eFrame.getColumn(2);
1052 >        Vector3d uz_i = sdat.eFrame.getColumn(2);
1053          Vector3d ei = preVal * uz_i;
1054  
1055          // This looks very wrong.  A vector crossed with itself is zero.
1056 <        scdat.t -= cross(uz_i, ei);
1056 >        sdat.t -= cross(uz_i, ei);
1057        }
1058      } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) {
1059        if (i_is_Charge) {        
1060          chg1 = data.charge;
1061          if (screeningMethod_ == DAMPED) {
1062 <          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1062 >          self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + sdat.skippedCharge) * pre11_;
1063          } else {        
1064 <          self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1064 >          self = - 0.5 * rcuti_ * chg1 * (chg1 + sdat.skippedCharge) * pre11_;
1065          }
1066 <        scdat.pot += self;
1066 >        sdat.pot[2] += self;
1067        }
1068      }
1069    }
1070  
1071 <  RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1071 >  RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1072      // This seems to work moderately well as a default.  There's no
1073      // inherent scale for 1/r interactions that we can standardize.
1074      // 12 angstroms seems to be a reasonably good guess for most

Diff Legend

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