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 1825 by gezelter, Wed Jan 9 19:27:52 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 499 | Line 504 | namespace OpenMD {
504        v33v.push_back(v33);
505        v34v.push_back(v34);
506        v35v.push_back(v35);
507 <
507 >      
508        v41v.push_back(v41);
509        v42v.push_back(v42);
510        v43v.push_back(v43);
# 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 693 | Line 698 | namespace OpenMD {
698      RealType pref;
699  
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;
703 >    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr;
704      Mat3x3d  QaQb;                                // Cross-interaction matrices
705  
706      RealType U(0.0);  // Potential
# Line 729 | Line 735 | namespace OpenMD {
735  
736      ri = 1.0 /  *(idat.rij);
737      Vector3d rhat =  *(idat.d)  * ri;
732 –    ri2 = ri * ri;
738        
739      // logicals
740  
# Line 746 | Line 751 | namespace OpenMD {
751      // Obtain all of the required radial function values from the
752      // spline structures:
753      
754 <    if (a_is_Charge && b_is_Charge) {
755 <      v01 = v01s->getValueAt( *(idat.rij) );
754 >    // needed for fields (and forces):
755 >    if (a_is_Charge || b_is_Charge) {
756        v02 = v02s->getValueAt( *(idat.rij) );
757      }
758 <    if ((a_is_Charge && b_is_Dipole) || (b_is_Charge && a_is_Dipole)) {
754 <      v11 = v11s->getValueAt( *(idat.rij) );
758 >    if (a_is_Dipole || b_is_Dipole) {
759        v12 = v12s->getValueAt( *(idat.rij) );
760        v13 = v13s->getValueAt( *(idat.rij) );
761      }
762 <    if ((a_is_Charge && b_is_Quadrupole) ||
763 <        (b_is_Charge && a_is_Quadrupole) ||
764 <        (a_is_Dipole && b_is_Dipole)) {
762 >    if (a_is_Quadrupole || b_is_Quadrupole) {
763 >      v23 = v23s->getValueAt( *(idat.rij) );
764 >      v24 = v24s->getValueAt( *(idat.rij) );
765 >    }
766 >
767 >    // needed for potentials (and torques):
768 >    if (a_is_Charge && b_is_Charge) {
769 >      v01 = v01s->getValueAt( *(idat.rij) );
770 >    }
771 >    if ((a_is_Charge && b_is_Dipole) || (b_is_Charge && a_is_Dipole)) {
772 >      v11 = v11s->getValueAt( *(idat.rij) );
773 >    }
774 >    if ((a_is_Charge && b_is_Quadrupole) || (b_is_Charge && a_is_Quadrupole)) {
775        v21 = v21s->getValueAt( *(idat.rij) );
776        v22 = v22s->getValueAt( *(idat.rij) );
777 +    } else if (a_is_Dipole && b_is_Dipole) {
778 +      v21 = v21s->getValueAt( *(idat.rij) );
779 +      v22 = v22s->getValueAt( *(idat.rij) );
780        v23 = v23s->getValueAt( *(idat.rij) );
781        v24 = v24s->getValueAt( *(idat.rij) );
782      }
# Line 780 | Line 797 | namespace OpenMD {
797        v46 = v46s->getValueAt( *(idat.rij) );
798      }
799  
783 –
800      // calculate the single-site contributions (fields, etc).
801      
802      if (a_is_Charge) {
# Line 941 | Line 957 | namespace OpenMD {
957          F  -= pref * (DadDb * rhat + rdDb * D_a + rdDa * D_b)*v23;
958          F  -= pref * (rdDa * rdDb) * v24 * rhat;
959          Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
960 <        Tb += pref * (-v21 * DaxDb + v22 * rdDa * rxDb);
960 >        Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
961  
962          // Even if we excluded this pair from direct interactions, we
963          // still have the reaction-field-mediated dipole-dipole
# Line 995 | Line 1011 | namespace OpenMD {
1011          F  += pref * (rdDb * rdQar * rhat * v35);
1012          Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31
1013                        + 2.0*rdDb*rxQar*v32);
1014 <        Tb += pref * ((trQa*rxDb + 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1014 >        Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1015        }
1016        if (b_is_Quadrupole) {
1017 <        pref = pre44_ * *(idat.electroMult);
1017 >        pref = pre44_ * *(idat.electroMult);  // yes
1018          QaQb = Q_a * Q_b;
1019          trQaQb = QaQb.trace();
1020          rQaQb = rhat * QaQb;
1021 <        QaQbr = QaQb * rhat;        
1021 >        QaQbr = QaQb * rhat;
1022          QaxQb = cross(Q_a, Q_b);
1023 <
1024 <        U  += pref * (trQa * trQb + 2.0*trQaQb) * v41;
1025 <        U  += pref * (trQa*rdQbr + trQb*rdQar  + 4.0*dot(rQa, Qbr)) * v42;
1023 >        rQaQbr = dot(rQa, Qbr);
1024 >        rQaxQbr = cross(rQa, Qbr);
1025 >        
1026 >        U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1027 >        U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
1028          U  += pref * (rdQar * rdQbr) * v43;
1029  
1030 <        F  += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44;
1031 <        F  += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44;
1032 <        F  += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44;
1015 <        F  += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat
1016 <                      + 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;
1030 >        F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44;
1031 >        F  += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45;
1032 >        F  += pref * rhat * (rdQar * rdQbr)*v46;
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 >        F  += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44;
1035 >        F  += pref * 4.0 * (rQaQb + QaQbr)*v44;
1036 >        F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45;
1037 >
1038 >        Ta += pref * (- 4.0 * QaxQb * v41);
1039 >        Ta += pref * (- 2.0 * trQb * cross(rQa, rhat)
1040 >                      + 4.0 * cross(rhat, QaQbr)
1041 >                      - 4.0 * rQaxQbr) * v42;
1042          Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
1043  
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;        
1044  
1045 +        Tb += pref * (+ 4.0 * QaxQb * v41);
1046 +        Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1047 +                      - 4.0 * cross(rQaQb, rhat)
1048 +                      + 4.0 * rQaxQbr) * v42;
1049 +        // Possible replacement for line 2 above:
1050 +        //             + 4.0 * cross(rhat, QbQar)
1051 +
1052 +        Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1053 +
1054 +        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
1055        }
1056      }
1057  
# Line 1069 | Line 1092 | namespace OpenMD {
1092          *(idat.t2) += *(idat.sw) * indirect_Tb;
1093      }
1094      return;
1095 <  }  
1095 >  }
1096      
1097    void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1098  
# Line 1112 | Line 1135 | namespace OpenMD {
1135        
1136      case esm_SHIFTED_FORCE:
1137      case esm_SHIFTED_POTENTIAL:
1138 <      if (i_is_Charge) {        
1139 <        self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1138 >      if (i_is_Charge) {
1139 >        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1140          (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1141        }
1142        break;

Diff Legend

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