--- branches/development/src/nonbonded/Electrostatic.cpp 2012/10/22 20:42:10 1808 +++ branches/development/src/nonbonded/Electrostatic.cpp 2013/02/20 15:39:39 1850 @@ -35,7 +35,7 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ @@ -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) { @@ -804,15 +808,18 @@ namespace OpenMD { if (idat.excluded) { *(idat.skippedCharge2) += C_a; + } else { + // only do the field if we're not excluded: + Eb -= C_a * pre11_ * v02 * rhat; } - Eb -= C_a * pre11_ * v02 * rhat; } if (a_is_Dipole) { D_a = *(idat.dipole1); rdDa = dot(rhat, D_a); rxDa = cross(rhat, D_a); - Eb -= pre12_ * (v13 * rdDa * rhat + v12 * D_a); + if (!idat.excluded) + Eb -= pre12_ * (v13 * rdDa * rhat + v12 * D_a); } if (a_is_Quadrupole) { @@ -822,7 +829,8 @@ namespace OpenMD { rQa = rhat * Q_a; rdQar = dot(rhat, Qar); rxQar = cross(rhat, Qar); - Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24); + if (!idat.excluded) + Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24); } if (b_is_Charge) { @@ -833,15 +841,18 @@ namespace OpenMD { if (idat.excluded) { *(idat.skippedCharge1) += C_b; + } else { + // only do the field if we're not excluded: + Ea += C_b * pre11_ * v02 * rhat; } - Ea += C_b * pre11_ * v02 * rhat; } if (b_is_Dipole) { D_b = *(idat.dipole2); rdDb = dot(rhat, D_b); rxDb = cross(rhat, D_b); - Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b); + if (!idat.excluded) + Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b); } if (b_is_Quadrupole) { @@ -851,7 +862,8 @@ namespace OpenMD { rQb = rhat * Q_b; rdQbr = dot(rhat, Qbr); rxQbr = cross(rhat, Qbr); - Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24); + if (!idat.excluded) + Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24); } if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { @@ -953,7 +965,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 +1019,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; + 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; - 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; + 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; - 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 * 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 +1100,7 @@ namespace OpenMD { *(idat.t2) += *(idat.sw) * indirect_Tb; } return; - } + } void Electrostatic::calcSelfCorrection(SelfData &sdat) { @@ -1124,8 +1143,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;