106 |
|
angstromToM_ = 1.0e-10; |
107 |
|
debyeToCm_ = 3.33564095198e-30; |
108 |
|
|
109 |
< |
// number of points for electrostatic splines |
110 |
< |
np_ = 1000; |
109 |
> |
// Default number of points for electrostatic splines |
110 |
> |
np_ = 140; |
111 |
|
|
112 |
|
// variables to handle different summation methods for long-range |
113 |
|
// electrostatics: |
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; |
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. |
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 |
735 |
|
|
736 |
|
ri = 1.0 / *(idat.rij); |
737 |
|
Vector3d rhat = *(idat.d) * ri; |
733 |
– |
ri2 = ri * ri; |
738 |
|
|
739 |
|
// logicals |
740 |
|
|
808 |
|
|
809 |
|
if (idat.excluded) { |
810 |
|
*(idat.skippedCharge2) += C_a; |
811 |
+ |
} else { |
812 |
+ |
// only do the field if we're not excluded: |
813 |
+ |
Eb -= C_a * pre11_ * v02 * rhat; |
814 |
|
} |
808 |
– |
Eb -= C_a * pre11_ * v02 * rhat; |
815 |
|
} |
816 |
|
|
817 |
|
if (a_is_Dipole) { |
818 |
|
D_a = *(idat.dipole1); |
819 |
|
rdDa = dot(rhat, D_a); |
820 |
|
rxDa = cross(rhat, D_a); |
821 |
< |
Eb -= pre12_ * (v13 * rdDa * rhat + v12 * D_a); |
821 |
> |
if (!idat.excluded) |
822 |
> |
Eb -= pre12_ * (v13 * rdDa * rhat + v12 * D_a); |
823 |
|
} |
824 |
|
|
825 |
|
if (a_is_Quadrupole) { |
829 |
|
rQa = rhat * Q_a; |
830 |
|
rdQar = dot(rhat, Qar); |
831 |
|
rxQar = cross(rhat, Qar); |
832 |
< |
Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24); |
832 |
> |
if (!idat.excluded) |
833 |
> |
Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24); |
834 |
|
} |
835 |
|
|
836 |
|
if (b_is_Charge) { |
841 |
|
|
842 |
|
if (idat.excluded) { |
843 |
|
*(idat.skippedCharge1) += C_b; |
844 |
+ |
} else { |
845 |
+ |
// only do the field if we're not excluded: |
846 |
+ |
Ea += C_b * pre11_ * v02 * rhat; |
847 |
|
} |
837 |
– |
Ea += C_b * pre11_ * v02 * rhat; |
848 |
|
} |
849 |
|
|
850 |
|
if (b_is_Dipole) { |
851 |
|
D_b = *(idat.dipole2); |
852 |
|
rdDb = dot(rhat, D_b); |
853 |
|
rxDb = cross(rhat, D_b); |
854 |
< |
Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b); |
854 |
> |
if (!idat.excluded) |
855 |
> |
Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b); |
856 |
|
} |
857 |
|
|
858 |
|
if (b_is_Quadrupole) { |
862 |
|
rQb = rhat * Q_b; |
863 |
|
rdQbr = dot(rhat, Qbr); |
864 |
|
rxQbr = cross(rhat, Qbr); |
865 |
< |
Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24); |
865 |
> |
if (!idat.excluded) |
866 |
> |
Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24); |
867 |
|
} |
868 |
|
|
869 |
|
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
1100 |
|
*(idat.t2) += *(idat.sw) * indirect_Tb; |
1101 |
|
} |
1102 |
|
return; |
1103 |
< |
} |
1103 |
> |
} |
1104 |
|
|
1105 |
|
void Electrostatic::calcSelfCorrection(SelfData &sdat) { |
1106 |
|
|
1143 |
|
|
1144 |
|
case esm_SHIFTED_FORCE: |
1145 |
|
case esm_SHIFTED_POTENTIAL: |
1146 |
< |
if (i_is_Charge) { |
1147 |
< |
self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; |
1146 |
> |
if (i_is_Charge) { |
1147 |
> |
self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; |
1148 |
|
(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; |
1149 |
|
} |
1150 |
|
break; |