--- branches/development/src/nonbonded/Electrostatic.cpp 2012/10/22 20:42:10 1808 +++ branches/development/src/nonbonded/Electrostatic.cpp 2013/01/09 19:27:52 1825 @@ -106,8 +106,8 @@ namespace OpenMD { angstromToM_ = 1.0e-10; debyeToCm_ = 3.33564095198e-30; - // number of points for electrostatic splines - np_ = 100; + // Default number of points for electrostatic splines + np_ = 140; // variables to handle different summation methods for long-range // electrostatics: @@ -246,7 +246,7 @@ namespace OpenMD { b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; - selfMult_ = b0c + 2.0 * a2 * invArootPi; + selfMult_ = b0c + a2 * invArootPi; } else { a2 = 0.0; b0c = 1.0 / r; @@ -279,6 +279,11 @@ namespace OpenMD { // working variables for Taylor expansion: RealType rmRc, rmRc2, rmRc3, rmRc4; + // Approximate using splines using a maximum of 0.1 Angstroms + // between points. + int nptest = int((cutoffRadius_ + 2.0) / 0.1); + np_ = (np_ > nptest) ? np_ : nptest; + // Add a 2 angstrom safety window to deal with cutoffGroups that // have charged atoms longer than the cutoffRadius away from each // other. Splining is almost certainly the best choice here. @@ -499,7 +504,7 @@ namespace OpenMD { v33v.push_back(v33); v34v.push_back(v34); v35v.push_back(v35); - + v41v.push_back(v41); v42v.push_back(v42); v43v.push_back(v43); @@ -685,7 +690,7 @@ namespace OpenMD { Vector3d D_a, D_b; // Dipoles (space-fixed) Mat3x3d Q_a, Q_b; // Quadrupoles (space-fixed) - RealType ri, ri2, ri3, ri4; // Distance utility scalars + RealType ri; // Distance utility scalar RealType rdDa, rdDb; // Dipole utility scalars Vector3d rxDa, rxDb; // Dipole utility vectors RealType rdQar, rdQbr, trQa, trQb; // Quadrupole utility scalars @@ -693,8 +698,9 @@ namespace OpenMD { RealType pref; RealType DadDb, trQaQb, DadQbr, DbdQar; // Cross-interaction scalars + RealType rQaQbr; Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors - Vector3d rQaQb, QaQbr, QaxQb; + Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr; Mat3x3d QaQb; // Cross-interaction matrices RealType U(0.0); // Potential @@ -729,7 +735,6 @@ namespace OpenMD { ri = 1.0 / *(idat.rij); Vector3d rhat = *(idat.d) * ri; - ri2 = ri * ri; // logicals @@ -792,7 +797,6 @@ namespace OpenMD { v46 = v46s->getValueAt( *(idat.rij) ); } - // calculate the single-site contributions (fields, etc). if (a_is_Charge) { @@ -953,7 +957,7 @@ namespace OpenMD { F -= pref * (DadDb * rhat + rdDb * D_a + rdDa * D_b)*v23; F -= pref * (rdDa * rdDb) * v24 * rhat; Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa); - Tb += pref * (-v21 * DaxDb + v22 * rdDa * rxDb); + Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb); // Even if we excluded this pair from direct interactions, we // still have the reaction-field-mediated dipole-dipole @@ -1007,40 +1011,47 @@ namespace OpenMD { F += pref * (rdDb * rdQar * rhat * v35); Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31 + 2.0*rdDb*rxQar*v32); - Tb += pref * ((trQa*rxDb + 2.0 * DbxQar)*v31 + rxDb*rdQar*v32); + Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32); } if (b_is_Quadrupole) { - pref = pre44_ * *(idat.electroMult); + pref = pre44_ * *(idat.electroMult); // yes QaQb = Q_a * Q_b; trQaQb = QaQb.trace(); rQaQb = rhat * QaQb; - QaQbr = QaQb * rhat; + QaQbr = QaQb * rhat; QaxQb = cross(Q_a, Q_b); - - U += pref * (trQa * trQb + 2.0*trQaQb) * v41; - U += pref * (trQa*rdQbr + trQb*rdQar + 4.0*dot(rQa, Qbr)) * v42; - U += pref * (rdQar * rdQbr) * v43; - - F += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44; - F += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44; - F += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44; - F += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat - + 4.0*dot(rQa, Qbr)*rhat)*v45; - F += pref * (2.0*rQa*rdQbr + 2.0*rdQar*rQb)*v45; - F += pref * (rdQar*rdQbr*rhat) * v46; + rQaQbr = dot(rQa, Qbr); + rQaxQbr = cross(rQa, Qbr); + + U += pref * (trQa * trQb + 2.0 * trQaQb) * v41; + U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42; + U += pref * (rdQar * rdQbr) * v43; - Ta += pref * (-4.0 * QaxQb * v41); - Ta += pref * (-2.0*trQb*cross(rQa, rhat) - + 4.0*cross(rhat, QaQbr) - - 4.0*cross(rQa, Qbr)) * v42; + F += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44; + F += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45; + F += pref * rhat * (rdQar * rdQbr)*v46; + + F += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44; + F += pref * 4.0 * (rQaQb + QaQbr)*v44; + F += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45; + + Ta += pref * (- 4.0 * QaxQb * v41); + Ta += pref * (- 2.0 * trQb * cross(rQa, rhat) + + 4.0 * cross(rhat, QaQbr) + - 4.0 * rQaxQbr) * v42; Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43; - Tb += pref * (4.0 * QaxQb * v41); - Tb += pref * (-2.0*trQa*cross(rQb, rhat) - - 4.0*cross(rQaQb, rhat) - + 4.0*cross(rQa, Qbr))*v42; - Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; + Tb += pref * (+ 4.0 * QaxQb * v41); + Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) + - 4.0 * cross(rQaQb, rhat) + + 4.0 * rQaxQbr) * v42; + // Possible replacement for line 2 above: + // + 4.0 * cross(rhat, QbQar) + + Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; + + // cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; } } @@ -1081,7 +1092,7 @@ namespace OpenMD { *(idat.t2) += *(idat.sw) * indirect_Tb; } return; - } + } void Electrostatic::calcSelfCorrection(SelfData &sdat) { @@ -1124,8 +1135,8 @@ namespace OpenMD { case esm_SHIFTED_FORCE: case esm_SHIFTED_POTENTIAL: - if (i_is_Charge) { - self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; + if (i_is_Charge) { + self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; } break;