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 1821 by gezelter, Mon Jan 7 20:05:43 2013 UTC vs.
Revision 1842 by gezelter, Tue Jan 29 19:10:04 2013 UTC

# Line 106 | Line 106 | namespace OpenMD {
106      angstromToM_ = 1.0e-10;
107      debyeToCm_ = 3.33564095198e-30;
108      
109 <    // number of points for electrostatic splines
110 <    np_ = 1000;
109 >    // Default number of points for electrostatic splines
110 >    np_ = 140;
111      
112      // variables to handle different summation methods for long-range
113      // electrostatics:
# Line 246 | Line 246 | namespace OpenMD {
246        b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
247        b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
248        b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
249 <      selfMult_ = b0c  +  2.0 * a2 * invArootPi;
249 >      selfMult_ = b0c + a2 * invArootPi;
250      } else {
251        a2 = 0.0;
252        b0c = 1.0 / r;
# Line 279 | Line 279 | namespace OpenMD {
279      // working variables for Taylor expansion:
280      RealType rmRc, rmRc2, rmRc3, rmRc4;
281  
282 +    // Approximate using splines using a maximum of 0.1 Angstroms
283 +    // between points.
284 +    int nptest = int((cutoffRadius_ + 2.0) / 0.1);
285 +    np_ = (np_ > nptest) ? np_ : nptest;
286 +  
287      // Add a 2 angstrom safety window to deal with cutoffGroups that
288      // have charged atoms longer than the cutoffRadius away from each
289      // other.  Splining is almost certainly the best choice here.
# Line 685 | Line 690 | namespace OpenMD {
690      Vector3d D_a, D_b;  // Dipoles (space-fixed)
691      Mat3x3d  Q_a, Q_b;  // Quadrupoles (space-fixed)
692  
693 <    RealType ri, ri2, ri3, ri4;                  // Distance utility scalars
693 >    RealType ri;                                 // Distance utility scalar
694      RealType rdDa, rdDb;                         // Dipole utility scalars
695      Vector3d rxDa, rxDb;                         // Dipole utility vectors
696      RealType rdQar, rdQbr, trQa, trQb;           // Quadrupole utility scalars
# Line 695 | Line 700 | namespace OpenMD {
700      RealType DadDb, trQaQb, DadQbr, DbdQar;       // Cross-interaction scalars
701      RealType rQaQbr;
702      Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors
703 <    Vector3d rQaQb, QaQbr, QaxQb, QbxQa, rQaxQbr, QbQar, rQbxQar;
704 <    Mat3x3d  QaQb, QbQa;                          // Cross-interaction matrices
703 >    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr;
704 >    Mat3x3d  QaQb;                                // Cross-interaction matrices
705  
706      RealType U(0.0);  // Potential
707      Vector3d F(0.0);  // Force
# Line 730 | Line 735 | namespace OpenMD {
735  
736      ri = 1.0 /  *(idat.rij);
737      Vector3d rhat =  *(idat.d)  * ri;
733 –    ri2 = ri * ri;
738        
739      // logicals
740  
# Line 793 | Line 797 | namespace OpenMD {
797        v46 = v46s->getValueAt( *(idat.rij) );
798      }
799  
796 –
800      // calculate the single-site contributions (fields, etc).
801      
802      if (a_is_Charge) {
# Line 805 | Line 808 | namespace OpenMD {
808        
809        if (idat.excluded) {
810          *(idat.skippedCharge2) += C_a;
811 +      } else {
812 +        // only do the field if we're not excluded:
813 +        Eb -= C_a *  pre11_ * v02 * rhat;
814        }
809 –      Eb -= C_a *  pre11_ * v02 * rhat;
815      }
816      
817      if (a_is_Dipole) {
818        D_a = *(idat.dipole1);
819        rdDa = dot(rhat, D_a);
820        rxDa = cross(rhat, D_a);
821 <      Eb -=  pre12_ * (v13 * rdDa * rhat + v12 * D_a);
821 >      if (!idat.excluded)
822 >        Eb -=  pre12_ * (v13 * rdDa * rhat + v12 * D_a);
823      }
824      
825      if (a_is_Quadrupole) {
# Line 823 | Line 829 | namespace OpenMD {
829        rQa = rhat * Q_a;
830        rdQar = dot(rhat, Qar);
831        rxQar = cross(rhat, Qar);
832 <      Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24);
832 >      if (!idat.excluded)
833 >        Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24);
834      }
835      
836      if (b_is_Charge) {
# Line 834 | Line 841 | namespace OpenMD {
841        
842        if (idat.excluded) {
843          *(idat.skippedCharge1) += C_b;
844 +      } else {
845 +        // only do the field if we're not excluded:
846 +        Ea += C_b *  pre11_ * v02 * rhat;
847        }
838 –      Ea += C_b *  pre11_ * v02 * rhat;
848      }
849      
850      if (b_is_Dipole) {
851        D_b = *(idat.dipole2);
852        rdDb = dot(rhat, D_b);
853        rxDb = cross(rhat, D_b);
854 <      Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b);
854 >      if (!idat.excluded)
855 >        Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b);
856      }
857      
858      if (b_is_Quadrupole) {
# Line 852 | Line 862 | namespace OpenMD {
862        rQb = rhat * Q_b;
863        rdQbr = dot(rhat, Qbr);
864        rxQbr = cross(rhat, Qbr);
865 <      Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24);
865 >      if (!idat.excluded)
866 >        Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24);
867      }
868      
869      if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
# Line 1011 | Line 1022 | namespace OpenMD {
1022          Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1023        }
1024        if (b_is_Quadrupole) {
1025 <        pref = pre44_ * *(idat.electroMult);
1025 >        pref = pre44_ * *(idat.electroMult);  // yes
1026          QaQb = Q_a * Q_b;
1016 –        QbQa = Q_b * Q_a;
1027          trQaQb = QaQb.trace();
1028          rQaQb = rhat * QaQb;
1029 <        QbQar = QbQa * rhat;
1020 <        QaQbr = QaQb * rhat;        
1029 >        QaQbr = QaQb * rhat;
1030          QaxQb = cross(Q_a, Q_b);
1022 –        QbxQa = cross(Q_b, Q_a);
1031          rQaQbr = dot(rQa, Qbr);
1032          rQaxQbr = cross(rQa, Qbr);
1025 –        //rQbxQar = cross(rQb, Qar);
1026 –
1027 –
1028 –        // First part of potential and associated forces:
1029 –        // U += pref * (trQa * trQb) * v41;
1030 –        // F += pref * rhat * (trQa * trQb) * v44;
1031 –        // Ta += 0.0;
1032 –        // Tb += 0.0;
1033 –
1034 –        // cerr << "Aa:\n";
1035 –        // cerr << *(idat.A1)  << "\n\n";
1036 –
1037 –        // cerr << "Ab:\n";
1038 –        // cerr << *(idat.A2)  << "\n\n";
1039 –
1040 –        // cerr << "Qa:\n";
1041 –        // cerr << Q_a << "\n\n";
1042 –        // cerr << "Qb:\n";
1043 –        // cerr << Q_b << "\n\n";
1044 –
1045 –        
1046 –        // Second part of potential and associated forces:
1047 –        // Probably suspect (particularly torques):
1048 –        // cerr << "trQaQb = " << trQaQb << "\n";
1049 –
1050 –        // U += pref * 2.0 * trQaQb * v41;
1051 –        // F += pref * rhat * (2.0 * trQaQb) * v44;
1052 –        // Ta += pref * (-4.0*QaxQb) * v41;
1053 –        // Tb += pref * (-4.0*QbxQa) * v41;
1054 –
1055 –        // cerr << "QaxQb = " << QaxQb << "\n";
1056 –        // cerr << "QbxQa = " << QbxQa << "\n";
1057 –
1058 –        // Third part of potential:
1059 –        // U += pref * (trQa * rdQbr) * v42;
1060 –        // F += pref * rhat * (trQa * rdQbr) * v45;
1061 –        // F += pref * 2.0 * (trQa * rQb) * v44;
1062 –        // Ta += 0.0;
1063 –        // Tb += pref * 2.0 * trQa * cross(rhat, Qbr) * v42;
1033          
1034          U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1035          U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
# Line 1083 | Line 1052 | namespace OpenMD {
1052  
1053          Tb += pref * (+ 4.0 * QaxQb * v41);
1054          Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1055 <                      + 4.0 * cross(rhat, QbQar)
1055 >                      - 4.0 * cross(rQaQb, rhat)
1056                        + 4.0 * rQaxQbr) * v42;
1057 +        // Possible replacement for line 2 above:
1058 +        //             + 4.0 * cross(rhat, QbQar)
1059 +
1060          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1061  
1090 –        // Tb += pref * (+ 4.0 * QaxQb * v41);
1091 –        // Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1092 –        //               - 4.0 * cross(rQaQb, rhat)
1093 –        //               + 4.0 * rQaxQbr) * v42;
1094 –        // Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1095 –
1062          //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1063        }
1064      }
# Line 1134 | Line 1100 | namespace OpenMD {
1100          *(idat.t2) += *(idat.sw) * indirect_Tb;
1101      }
1102      return;
1103 <  }  
1103 >  }
1104      
1105    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1106  
# Line 1177 | Line 1143 | namespace OpenMD {
1143        
1144      case esm_SHIFTED_FORCE:
1145      case esm_SHIFTED_POTENTIAL:
1146 <      if (i_is_Charge) {        
1147 <        self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1146 >      if (i_is_Charge) {
1147 >        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1148          (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1149        }
1150        break;

Diff Legend

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