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 1808 by gezelter, Mon Oct 22 20:42:10 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 499 | Line 499 | namespace OpenMD {
499        v33v.push_back(v33);
500        v34v.push_back(v34);
501        v35v.push_back(v35);
502 <
502 >      
503        v41v.push_back(v41);
504        v42v.push_back(v42);
505        v43v.push_back(v43);
# Line 693 | Line 693 | namespace OpenMD {
693      RealType pref;
694  
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;
698 >    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr;
699      Mat3x3d  QaQb;                                // Cross-interaction matrices
700  
701      RealType U(0.0);  // Potential
# Line 791 | Line 792 | namespace OpenMD {
792        v45 = v45s->getValueAt( *(idat.rij) );
793        v46 = v46s->getValueAt( *(idat.rij) );
794      }
794 –
795  
796      // calculate the single-site contributions (fields, etc).
797      
# Line 953 | Line 953 | namespace OpenMD {
953          F  -= pref * (DadDb * rhat + rdDb * D_a + rdDa * D_b)*v23;
954          F  -= pref * (rdDa * rdDb) * v24 * rhat;
955          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
956 <        Tb += pref * (-v21 * DaxDb + v22 * rdDa * rxDb);
956 >        Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
957  
958          // Even if we excluded this pair from direct interactions, we
959          // still have the reaction-field-mediated dipole-dipole
# Line 1007 | Line 1007 | namespace OpenMD {
1007          F  += pref * (rdDb * rdQar * rhat * v35);
1008          Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31
1009                        + 2.0*rdDb*rxQar*v32);
1010 <        Tb += pref * ((trQa*rxDb + 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
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;
1015          trQaQb = QaQb.trace();
1016          rQaQb = rhat * QaQb;
1017 <        QaQbr = QaQb * rhat;        
1017 >        QaQbr = QaQb * rhat;
1018          QaxQb = cross(Q_a, Q_b);
1019 <
1020 <        U  += pref * (trQa * trQb + 2.0*trQaQb) * v41;
1021 <        U  += pref * (trQa*rdQbr + trQb*rdQar  + 4.0*dot(rQa, Qbr)) * v42;
1019 >        rQaQbr = dot(rQa, Qbr);
1020 >        rQaxQbr = cross(rQa, Qbr);
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;
1025  
1026 <        F  += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44;
1027 <        F  += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44;
1028 <        F  += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44;
1027 <        F  += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat
1028 <                      + 4.0*dot(rQa, Qbr)*rhat)*v45;
1029 <        F  += pref * (2.0*rQa*rdQbr + 2.0*rdQar*rQb)*v45;
1030 <        F  += pref * (rdQar*rdQbr*rhat) * v46;
1026 >        F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44;
1027 >        F  += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45;
1028 >        F  += pref * rhat * (rdQar * rdQbr)*v46;
1029  
1030 <        Ta += pref * (-4.0 * QaxQb  * v41);
1031 <        Ta += pref * (-2.0*trQb*cross(rQa, rhat)
1032 <                      + 4.0*cross(rhat, QaQbr)
1033 <                      - 4.0*cross(rQa, Qbr)) * v42;
1030 >        F  += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44;
1031 >        F  += pref * 4.0 * (rQaQb + QaQbr)*v44;
1032 >        F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45;
1033 >
1034 >        Ta += pref * (- 4.0 * QaxQb * v41);
1035 >        Ta += pref * (- 2.0 * trQb * cross(rQa, rhat)
1036 >                      + 4.0 * cross(rhat, QaQbr)
1037 >                      - 4.0 * rQaxQbr) * v42;
1038          Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
1039  
1038 –        Tb += pref * (4.0 * QaxQb * v41);
1039 –        Tb += pref * (-2.0*trQa*cross(rQb, rhat)
1040 –                      - 4.0*cross(rQaQb, rhat)
1041 –                      + 4.0*cross(rQa, Qbr))*v42;
1042 –        Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;        
1040  
1041 +        Tb += pref * (+ 4.0 * QaxQb * v41);
1042 +        Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
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 +        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1051        }
1052      }
1053  

Diff Legend

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