| 765 |
|
|
| 766 |
|
// Excluded potential that is still computed for fluctuating charges |
| 767 |
|
excluded_Pot= 0.0; |
| 768 |
– |
|
| 768 |
|
|
| 769 |
|
// some variables we'll need independent of electrostatic type: |
| 770 |
|
|
| 886 |
|
Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or |
| 887 |
|
+ rdQbr * rhat * (dv22 - 2.0*v22or)); |
| 888 |
|
} |
| 889 |
< |
|
| 889 |
> |
|
| 890 |
> |
|
| 891 |
|
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
| 892 |
|
J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; |
| 893 |
|
} |
| 894 |
< |
|
| 894 |
> |
|
| 895 |
|
if (a_is_Charge) { |
| 896 |
|
|
| 897 |
|
if (b_is_Charge) { |
| 898 |
|
pref = pre11_ * *(idat.electroMult); |
| 899 |
|
U += C_a * C_b * pref * v01; |
| 900 |
|
F += C_a * C_b * pref * dv01 * rhat; |
| 901 |
< |
|
| 901 |
> |
|
| 902 |
|
// If this is an excluded pair, there are still indirect |
| 903 |
|
// interactions via the reaction field we must worry about: |
| 904 |
|
|
| 907 |
|
indirect_Pot += rfContrib; |
| 908 |
|
indirect_F += rfContrib * 2.0 * ri * rhat; |
| 909 |
|
} |
| 910 |
< |
|
| 910 |
> |
|
| 911 |
|
// Fluctuating charge forces are handled via Coulomb integrals |
| 912 |
|
// for excluded pairs (i.e. those connected via bonds) and |
| 913 |
|
// with the standard charge-charge interaction otherwise. |
| 914 |
|
|
| 915 |
< |
if (idat.excluded) { |
| 915 |
> |
if (idat.excluded) { |
| 916 |
|
if (a_is_Fluctuating || b_is_Fluctuating) { |
| 917 |
|
coulInt = J->getValueAt( *(idat.rij) ); |
| 918 |
< |
if (a_is_Fluctuating) dUdCa += coulInt * C_b; |
| 919 |
< |
if (b_is_Fluctuating) dUdCb += coulInt * C_a; |
| 920 |
< |
excluded_Pot += C_a * C_b * coulInt; |
| 921 |
< |
} |
| 918 |
> |
if (a_is_Fluctuating) dUdCa += C_b * coulInt; |
| 919 |
> |
if (b_is_Fluctuating) dUdCb += C_a * coulInt; |
| 920 |
> |
} |
| 921 |
|
} else { |
| 922 |
|
if (a_is_Fluctuating) dUdCa += C_b * pref * v01; |
| 923 |
|
if (b_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 924 |
< |
} |
| 924 |
> |
} |
| 925 |
|
} |
| 926 |
|
|
| 927 |
|
if (b_is_Dipole) { |
| 1141 |
|
|
| 1142 |
|
if (i_is_Fluctuating) { |
| 1143 |
|
C_a += *(sdat.flucQ); |
| 1144 |
< |
// dVdFQ is really a force, so this is negative the derivative |
| 1146 |
< |
*(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1144 |
> |
*(sdat.flucQfrc) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1145 |
|
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * |
| 1146 |
|
(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); |
| 1147 |
|
} |
| 1439 |
|
dks[i] = dk * skr[i]; |
| 1440 |
|
} |
| 1441 |
|
if (data.is_Quadrupole) { |
| 1442 |
< |
Q = atom->getQuadrupole(); |
| 1443 |
< |
Q *= mPoleConverter; |
| 1444 |
< |
Qk = Q * kVec; |
| 1445 |
< |
qk = dot(kVec, Qk); |
| 1448 |
< |
qxk[i] = cross(kVec, Qk); |
| 1442 |
> |
Q = atom->getQuadrupole() * mPoleConverter; |
| 1443 |
> |
Qk = Q * kVec; |
| 1444 |
> |
qk = dot(Qk, kVec); |
| 1445 |
> |
qxk[i] = cross(Qk, kVec); |
| 1446 |
|
qkc[i] = qk * ckr[i]; |
| 1447 |
|
qks[i] = qk * skr[i]; |
| 1448 |
|
} |
