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 1787 by gezelter, Wed Aug 29 18:13:11 2012 UTC vs.
Revision 1821 by gezelter, Mon Jan 7 20:05:43 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;
699 <    Mat3x3d  QaQb;                                // Cross-interaction matrices
698 >    Vector3d rQaQb, QaQbr, QaxQb, QbxQa, rQaxQbr, QbQar, rQbxQar;
699 >    Mat3x3d  QaQb, QbQa;                          // Cross-interaction matrices
700  
701      RealType U(0.0);  // Potential
702      Vector3d F(0.0);  // Force
# Line 746 | Line 747 | namespace OpenMD {
747      // Obtain all of the required radial function values from the
748      // spline structures:
749      
750 <    if (a_is_Charge && b_is_Charge) {
751 <      v01 = v01s->getValueAt( *(idat.rij) );
750 >    // needed for fields (and forces):
751 >    if (a_is_Charge || b_is_Charge) {
752        v02 = v02s->getValueAt( *(idat.rij) );
753      }
754 <    if ((a_is_Charge && b_is_Dipole) || (b_is_Charge && a_is_Dipole)) {
754 <      v11 = v11s->getValueAt( *(idat.rij) );
754 >    if (a_is_Dipole || b_is_Dipole) {
755        v12 = v12s->getValueAt( *(idat.rij) );
756        v13 = v13s->getValueAt( *(idat.rij) );
757      }
758 <    if ((a_is_Charge && b_is_Quadrupole) ||
759 <        (b_is_Charge && a_is_Quadrupole) ||
760 <        (a_is_Dipole && b_is_Dipole)) {
758 >    if (a_is_Quadrupole || b_is_Quadrupole) {
759 >      v23 = v23s->getValueAt( *(idat.rij) );
760 >      v24 = v24s->getValueAt( *(idat.rij) );
761 >    }
762 >
763 >    // needed for potentials (and torques):
764 >    if (a_is_Charge && b_is_Charge) {
765 >      v01 = v01s->getValueAt( *(idat.rij) );
766 >    }
767 >    if ((a_is_Charge && b_is_Dipole) || (b_is_Charge && a_is_Dipole)) {
768 >      v11 = v11s->getValueAt( *(idat.rij) );
769 >    }
770 >    if ((a_is_Charge && b_is_Quadrupole) || (b_is_Charge && a_is_Quadrupole)) {
771        v21 = v21s->getValueAt( *(idat.rij) );
772        v22 = v22s->getValueAt( *(idat.rij) );
773 +    } else if (a_is_Dipole && b_is_Dipole) {
774 +      v21 = v21s->getValueAt( *(idat.rij) );
775 +      v22 = v22s->getValueAt( *(idat.rij) );
776        v23 = v23s->getValueAt( *(idat.rij) );
777        v24 = v24s->getValueAt( *(idat.rij) );
778      }
# Line 941 | Line 954 | namespace OpenMD {
954          F  -= pref * (DadDb * rhat + rdDb * D_a + rdDa * D_b)*v23;
955          F  -= pref * (rdDa * rdDb) * v24 * rhat;
956          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
957 <        Tb += pref * (-v21 * DaxDb + v22 * rdDa * rxDb);
957 >        Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
958  
959          // Even if we excluded this pair from direct interactions, we
960          // still have the reaction-field-mediated dipole-dipole
# Line 995 | Line 1008 | namespace OpenMD {
1008          F  += pref * (rdDb * rdQar * rhat * v35);
1009          Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31
1010                        + 2.0*rdDb*rxQar*v32);
1011 <        Tb += pref * ((trQa*rxDb + 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1011 >        Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1012        }
1013        if (b_is_Quadrupole) {
1014          pref = pre44_ * *(idat.electroMult);
1015          QaQb = Q_a * Q_b;
1016 +        QbQa = Q_b * Q_a;
1017          trQaQb = QaQb.trace();
1018          rQaQb = rhat * QaQb;
1019 +        QbQar = QbQa * rhat;
1020          QaQbr = QaQb * rhat;        
1021          QaxQb = cross(Q_a, Q_b);
1022 +        QbxQa = cross(Q_b, Q_a);
1023 +        rQaQbr = dot(rQa, Qbr);
1024 +        rQaxQbr = cross(rQa, Qbr);
1025 +        //rQbxQar = cross(rQb, Qar);
1026  
1008 –        U  += pref * (trQa * trQb + 2.0*trQaQb) * v41;
1009 –        U  += pref * (trQa*rdQbr + trQb*rdQar  + 4.0*dot(rQa, Qbr)) * v42;
1010 –        U  += pref * (rdQar * rdQbr) * v43;
1027  
1028 <        F  += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44;
1029 <        F  += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44;
1030 <        F  += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44;
1031 <        F  += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat
1032 <                      + 4.0*dot(rQa, Qbr)*rhat)*v45;
1017 <        F  += pref * (2.0*rQa*rdQbr + 2.0*rdQar*rQb)*v45;
1018 <        F  += pref * (rdQar*rdQbr*rhat) * v46;
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 <        Ta += pref * (-4.0 * QaxQb  * v41);
1035 <        Ta += pref * (-2.0*trQb*cross(rQa, rhat)
1036 <                      + 4.0*cross(rhat, QaQbr)
1037 <                      - 4.0*cross(rQa, Qbr)) * v42;
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;
1064 >        
1065 >        U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1066 >        U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
1067 >        U  += pref * (rdQar * rdQbr) * v43;
1068 >
1069 >        F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44;
1070 >        F  += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45;
1071 >        F  += pref * rhat * (rdQar * rdQbr)*v46;
1072 >
1073 >        F  += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44;
1074 >        F  += pref * 4.0 * (rQaQb + QaQbr)*v44;
1075 >        F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45;
1076 >
1077 >        Ta += pref * (- 4.0 * QaxQb * v41);
1078 >        Ta += pref * (- 2.0 * trQb * cross(rQa, rhat)
1079 >                      + 4.0 * cross(rhat, QaQbr)
1080 >                      - 4.0 * rQaxQbr) * v42;
1081          Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
1082  
1026 –        Tb += pref * (4.0 * QaxQb * v41);
1027 –        Tb += pref * (-2.0*trQa*cross(rQb, rhat)
1028 –                      - 4.0*cross(rQaQb, rhat)
1029 –                      + 4.0*cross(rQa, Qbr))*v42;
1030 –        Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;        
1083  
1084 +        Tb += pref * (+ 4.0 * QaxQb * v41);
1085 +        Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1086 +                      + 4.0 * cross(rhat, QbQar)
1087 +                      + 4.0 * rQaxQbr) * v42;
1088 +        Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1089 +
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 +
1096 +        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1097        }
1098      }
1099  

Diff Legend

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