| 888 |
|
+ rdQbr * rhat * (dv22 - 2.0*v22or)); |
| 889 |
|
} |
| 890 |
|
|
| 891 |
< |
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
| 891 |
> |
if ((a_is_Fluctuating || b_is_Fluctuating) |
| 892 |
> |
&& (idat.excluded || idat.sameRegion)) { |
| 893 |
|
J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; |
| 894 |
|
} |
| 895 |
|
|
| 910 |
|
} |
| 911 |
|
|
| 912 |
|
// Fluctuating charge forces are handled via Coulomb integrals |
| 913 |
< |
// for excluded pairs (i.e. those connected via bonds) and |
| 914 |
< |
// with the standard charge-charge interaction otherwise. |
| 913 |
> |
// for excluded pairs (i.e. those connected via bonds), atoms |
| 914 |
> |
// within the same region (for metals) and with the standard |
| 915 |
> |
// charge-charge interaction otherwise. |
| 916 |
|
|
| 917 |
< |
if (idat.excluded) { |
| 917 |
> |
if (idat.excluded || idat.sameRegion) { |
| 918 |
|
if (a_is_Fluctuating || b_is_Fluctuating) { |
| 919 |
|
coulInt = J->getValueAt( *(idat.rij) ); |
| 920 |
|
if (a_is_Fluctuating) dUdCa += coulInt * C_b; |
| 921 |
|
if (b_is_Fluctuating) dUdCb += coulInt * C_a; |
| 922 |
< |
excluded_Pot += C_a * C_b * coulInt; |
| 922 |
> |
if (idat.excluded) |
| 923 |
> |
excluded_Pot += C_a * C_b * coulInt; |
| 924 |
|
} |
| 925 |
|
} else { |
| 926 |
|
if (a_is_Fluctuating) dUdCa += C_b * pref * v01; |
| 927 |
< |
if (a_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 927 |
> |
if (b_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 928 |
|
} |
| 929 |
|
} |
| 930 |
|
|
| 1145 |
|
|
| 1146 |
|
if (i_is_Fluctuating) { |
| 1147 |
|
C_a += *(sdat.flucQ); |
| 1148 |
< |
// dVdFQ is really a force, so this is negative the derivative |
| 1146 |
< |
*(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1148 |
> |
*(sdat.flucQfrc) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1149 |
|
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * |
| 1150 |
|
(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); |
| 1151 |
|
} |
| 1443 |
|
dks[i] = dk * skr[i]; |
| 1444 |
|
} |
| 1445 |
|
if (data.is_Quadrupole) { |
| 1446 |
< |
Q = atom->getQuadrupole(); |
| 1447 |
< |
Q *= mPoleConverter; |
| 1446 |
< |
Qk = Q * kVec; |
| 1446 |
> |
Q = atom->getQuadrupole() * mPoleConverter; |
| 1447 |
> |
Qk = Q * kVec; |
| 1448 |
|
qk = dot(kVec, Qk); |
| 1449 |
< |
qxk[i] = cross(kVec, Qk); |
| 1449 |
> |
qxk[i] = cross(Qk, kVec); |
| 1450 |
|
qkc[i] = qk * ckr[i]; |
| 1451 |
|
qks[i] = qk * skr[i]; |
| 1452 |
|
} |
| 1507 |
|
RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) |
| 1508 |
|
+skr[i]*(ckss+dkcs-qkss)); |
| 1509 |
|
|
| 1510 |
+ |
|
| 1511 |
+ |
cerr << "qkcs = " << qkcs << "\n"; |
| 1512 |
+ |
cerr << "qkss = " << qkss << "\n"; |
| 1513 |
+ |
cerr << "qkc[i] = " << qkc[i] << "\n"; |
| 1514 |
+ |
cerr << "qks[i] = " << qks[i] << "\n"; |
| 1515 |
+ |
cerr << "pot = " << 2.0 * rvol * AK[kk]*(qkss*qkss + qkcs*qkcs) << "\n"; |
| 1516 |
+ |
cerr << "qfrc = " << 4.0 * rvol * qfrc << "\n"; |
| 1517 |
+ |
cerr << "qtrq2 = " << 4.0 * rvol * qtrq2 << "\n"; |
| 1518 |
|
atom->addFrc( 4.0 * rvol * qfrc * kVec ); |
| 1519 |
|
|
| 1520 |
|
if (data.is_Dipole) { |
