| 262 |
|
b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; |
| 263 |
|
b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; |
| 264 |
|
b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; |
| 265 |
< |
selfMult_ = b0c + a2 * invArootPi; |
| 265 |
> |
//selfMult1_ = - 2.0 * a2 * invArootPi; |
| 266 |
> |
//selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0; |
| 267 |
> |
//selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0; |
| 268 |
> |
// Half the Smith self piece: |
| 269 |
> |
selfMult1_ = - a2 * invArootPi; |
| 270 |
> |
selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0; |
| 271 |
> |
selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0; |
| 272 |
|
} else { |
| 273 |
|
a2 = 0.0; |
| 274 |
|
b0c = 1.0 / r; |
| 277 |
|
b3c = (5.0 * b2c) / r2; |
| 278 |
|
b4c = (7.0 * b3c) / r2; |
| 279 |
|
b5c = (9.0 * b4c) / r2; |
| 280 |
< |
selfMult_ = b0c; |
| 280 |
> |
selfMult1_ = 0.0; |
| 281 |
> |
selfMult2_ = 0.0; |
| 282 |
> |
selfMult4_ = 0.0; |
| 283 |
|
} |
| 284 |
|
|
| 285 |
|
// higher derivatives of B_0 at the cutoff radius: |
| 287 |
|
db0c_2 = -b1c + r2 * b2c; |
| 288 |
|
db0c_3 = 3.0*r*b2c - r2*r*b3c; |
| 289 |
|
db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c; |
| 290 |
< |
db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; |
| 290 |
> |
db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; |
| 291 |
|
|
| 292 |
+ |
selfMult1_ -= b0c; |
| 293 |
+ |
selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; |
| 294 |
+ |
selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; |
| 295 |
|
|
| 296 |
|
// working variables for the splines: |
| 297 |
|
RealType ri, ri2; |
| 1115 |
|
|
| 1116 |
|
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1117 |
|
|
| 1107 |
– |
// cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; |
| 1118 |
|
} |
| 1119 |
|
} |
| 1120 |
|
|
| 1166 |
|
// logicals |
| 1167 |
|
bool i_is_Charge = data.is_Charge; |
| 1168 |
|
bool i_is_Dipole = data.is_Dipole; |
| 1169 |
+ |
bool i_is_Quadrupole = data.is_Quadrupole; |
| 1170 |
|
bool i_is_Fluctuating = data.is_Fluctuating; |
| 1171 |
|
RealType C_a = data.fixedCharge; |
| 1172 |
< |
RealType self, preVal, DadDa; |
| 1173 |
< |
|
| 1172 |
> |
RealType self(0.0), preVal, DdD, trQ, trQQ; |
| 1173 |
> |
|
| 1174 |
> |
if (i_is_Dipole) { |
| 1175 |
> |
DdD = data.dipole.lengthSquare(); |
| 1176 |
> |
} |
| 1177 |
> |
|
| 1178 |
|
if (i_is_Fluctuating) { |
| 1179 |
|
C_a += *(sdat.flucQ); |
| 1180 |
|
// dVdFQ is really a force, so this is negative the derivative |
| 1195 |
|
} |
| 1196 |
|
|
| 1197 |
|
if (i_is_Dipole) { |
| 1198 |
< |
DadDa = data.dipole.lengthSquare(); |
| 1184 |
< |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa; |
| 1198 |
> |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD; |
| 1199 |
|
} |
| 1200 |
|
|
| 1201 |
|
break; |
| 1202 |
|
|
| 1203 |
|
case esm_SHIFTED_FORCE: |
| 1204 |
|
case esm_SHIFTED_POTENTIAL: |
| 1205 |
< |
if (i_is_Charge) { |
| 1206 |
< |
self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; |
| 1207 |
< |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; |
| 1205 |
> |
case esm_TAYLOR_SHIFTED: |
| 1206 |
> |
if (i_is_Charge) |
| 1207 |
> |
self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); |
| 1208 |
> |
if (i_is_Dipole) |
| 1209 |
> |
self += selfMult2_ * pre22_ * DdD; |
| 1210 |
> |
if (i_is_Quadrupole) { |
| 1211 |
> |
trQ = data.quadrupole.trace(); |
| 1212 |
> |
trQQ = (data.quadrupole * data.quadrupole).trace(); |
| 1213 |
> |
self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ); |
| 1214 |
> |
if (i_is_Charge) |
| 1215 |
> |
self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ; |
| 1216 |
|
} |
| 1217 |
+ |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; |
| 1218 |
|
break; |
| 1219 |
|
default: |
| 1220 |
|
break; |
