--- branches/development/src/nonbonded/Electrostatic.cpp 2012/12/18 16:23:13 1820 +++ 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. @@ -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 @@ -695,8 +700,8 @@ namespace OpenMD { RealType DadDb, trQaQb, DadQbr, DbdQar; // Cross-interaction scalars RealType rQaQbr; Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors - Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr, QbQar, rQbxQar; - Mat3x3d QaQb, QbQa; // Cross-interaction matrices + Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr; + Mat3x3d QaQb; // Cross-interaction matrices RealType U(0.0); // Potential Vector3d F(0.0); // Force @@ -730,7 +735,6 @@ namespace OpenMD { ri = 1.0 / *(idat.rij); Vector3d rhat = *(idat.d) * ri; - ri2 = ri * ri; // logicals @@ -793,7 +797,6 @@ namespace OpenMD { v46 = v46s->getValueAt( *(idat.rij) ); } - // calculate the single-site contributions (fields, etc). if (a_is_Charge) { @@ -805,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) { @@ -823,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) { @@ -834,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) { @@ -852,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) { @@ -1011,18 +1022,15 @@ namespace OpenMD { 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; - QbQa = Q_b * Q_a; trQaQb = QaQb.trace(); rQaQb = rhat * QaQb; - QbQar = QbQa * rhat; - QaQbr = QaQb * rhat; + QaQbr = QaQb * rhat; QaxQb = cross(Q_a, Q_b); rQaQbr = dot(rQa, Qbr); rQaxQbr = cross(rQa, Qbr); - rQbxQar = cross(rQb, Qar); - + U += pref * (trQa * trQb + 2.0 * trQaQb) * v41; U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42; U += pref * (rdQar * rdQbr) * v43; @@ -1044,17 +1052,14 @@ namespace OpenMD { Tb += pref * (+ 4.0 * QaxQb * v41); Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) - + 4.0 * cross(rhat, QbQar) - - 4.0 * rQbxQar) * v42; + - 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; - // Tb += pref * (+ 4.0 * QaxQb * v41); - // Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) - // - 4.0 * cross(rQaQb, rhat) - // + 4.0 * rQaxQbr) * v42; - // Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - - cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; + // cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; } } @@ -1095,7 +1100,7 @@ namespace OpenMD { *(idat.t2) += *(idat.sw) * indirect_Tb; } return; - } + } void Electrostatic::calcSelfCorrection(SelfData &sdat) { @@ -1138,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;