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 1822 by gezelter, Tue Jan 8 15:29:03 2013 UTC

# Line 107 | Line 107 | namespace OpenMD {
107      debyeToCm_ = 3.33564095198e-30;
108      
109      // number of points for electrostatic splines
110 <    np_ = 100;
110 >    np_ = 1000;
111      
112      // variables to handle different summation methods for long-range
113      // electrostatics:
# Line 695 | Line 695 | namespace OpenMD {
695      RealType DadDb, trQaQb, DadQbr, DbdQar;       // Cross-interaction scalars
696      RealType rQaQbr;
697      Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors
698 <    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr, QbQar, rQbxQar;
699 <    Mat3x3d  QaQb, QbQa;                          // Cross-interaction matrices
698 >    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr;
699 >    Mat3x3d  QaQb;                                // Cross-interaction matrices
700  
701      RealType U(0.0);  // Potential
702      Vector3d F(0.0);  // Force
# Line 793 | Line 793 | namespace OpenMD {
793        v46 = v46s->getValueAt( *(idat.rij) );
794      }
795  
796 –
796      // calculate the single-site contributions (fields, etc).
797      
798      if (a_is_Charge) {
# Line 1011 | Line 1010 | namespace OpenMD {
1010          Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1011        }
1012        if (b_is_Quadrupole) {
1013 <        pref = pre44_ * *(idat.electroMult);
1013 >        pref = pre44_ * *(idat.electroMult);  // yes
1014          QaQb = Q_a * Q_b;
1016 –        QbQa = Q_b * Q_a;
1015          trQaQb = QaQb.trace();
1016          rQaQb = rhat * QaQb;
1017 <        QbQar = QbQa * rhat;
1020 <        QaQbr = QaQb * rhat;        
1017 >        QaQbr = QaQb * rhat;
1018          QaxQb = cross(Q_a, Q_b);
1019          rQaQbr = dot(rQa, Qbr);
1020          rQaxQbr = cross(rQa, Qbr);
1021 <        rQbxQar = cross(rQb, Qar);
1025 <
1021 >        
1022          U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1023          U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
1024          U  += pref * (rdQar * rdQbr) * v43;
# Line 1044 | Line 1040 | namespace OpenMD {
1040  
1041          Tb += pref * (+ 4.0 * QaxQb * v41);
1042          Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1043 <                      + 4.0 * cross(rhat, QbQar)
1044 <                      - 4.0 * rQbxQar) * v42;
1043 >                      - 4.0 * cross(rQaQb, rhat)
1044 >                      + 4.0 * rQaxQbr) * v42;
1045 >        // Possible replacement for line 2 above:
1046 >        //             + 4.0 * cross(rhat, QbQar)
1047 >
1048          Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1049  
1050 <        // 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";
1050 >        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1051        }
1052      }
1053  

Diff Legend

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