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 1820 by gezelter, Tue Dec 18 16:23:13 2012 UTC vs.
Revision 1823 by gezelter, Tue Jan 8 21:02:27 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_ = 100;
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 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, 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 793 | Line 798 | namespace OpenMD {
798        v46 = v46s->getValueAt( *(idat.rij) );
799      }
800  
796 –
801      // calculate the single-site contributions (fields, etc).
802      
803      if (a_is_Charge) {
# Line 1011 | Line 1015 | namespace OpenMD {
1015          Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1016        }
1017        if (b_is_Quadrupole) {
1018 <        pref = pre44_ * *(idat.electroMult);
1018 >        pref = pre44_ * *(idat.electroMult);  // yes
1019          QaQb = Q_a * Q_b;
1016 –        QbQa = Q_b * Q_a;
1020          trQaQb = QaQb.trace();
1021          rQaQb = rhat * QaQb;
1022 <        QbQar = QbQa * rhat;
1020 <        QaQbr = QaQb * rhat;        
1022 >        QaQbr = QaQb * rhat;
1023          QaxQb = cross(Q_a, Q_b);
1024          rQaQbr = dot(rQa, Qbr);
1025          rQaxQbr = cross(rQa, Qbr);
1026 <        rQbxQar = cross(rQb, Qar);
1025 <
1026 >        
1027          U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1028          U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
1029          U  += pref * (rdQar * rdQbr) * v43;
# Line 1044 | Line 1045 | namespace OpenMD {
1045  
1046          Tb += pref * (+ 4.0 * QaxQb * v41);
1047          Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1048 <                      + 4.0 * cross(rhat, QbQar)
1049 <                      - 4.0 * rQbxQar) * v42;
1048 >                      - 4.0 * cross(rQaQb, rhat)
1049 >                      + 4.0 * rQaxQbr) * v42;
1050 >        // Possible replacement for line 2 above:
1051 >        //             + 4.0 * cross(rhat, QbQar)
1052 >
1053          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1054  
1055 <        // Tb += pref * (+ 4.0 * QaxQb * v41);
1052 <        // Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1053 <        //               - 4.0 * cross(rQaQb, rhat)
1054 <        //               + 4.0 * rQaxQbr) * v42;
1055 <        // Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1056 <
1057 <        cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1055 >        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1056        }
1057      }
1058  
# Line 1095 | Line 1093 | namespace OpenMD {
1093          *(idat.t2) += *(idat.sw) * indirect_Tb;
1094      }
1095      return;
1096 <  }  
1096 >  }
1097      
1098    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1099  
# Line 1138 | Line 1136 | namespace OpenMD {
1136        
1137      case esm_SHIFTED_FORCE:
1138      case esm_SHIFTED_POTENTIAL:
1139 <      if (i_is_Charge) {        
1140 <        self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1139 >      if (i_is_Charge) {
1140 >        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1141          (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1142        }
1143        break;

Diff Legend

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