| 40 |
|
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 42 |
|
|
| 43 |
+ |
#ifdef IS_MPI |
| 44 |
+ |
#include <mpi.h> |
| 45 |
+ |
#endif |
| 46 |
+ |
|
| 47 |
|
#include <stdio.h> |
| 48 |
|
#include <string.h> |
| 49 |
|
|
| 61 |
|
#include "math/erfc.hpp" |
| 62 |
|
#include "math/SquareMatrix.hpp" |
| 63 |
|
#include "primitives/Molecule.hpp" |
| 60 |
– |
#ifdef IS_MPI |
| 61 |
– |
#include <mpi.h> |
| 62 |
– |
#endif |
| 64 |
|
|
| 65 |
|
namespace OpenMD { |
| 66 |
|
|
| 766 |
|
|
| 767 |
|
// Excluded potential that is still computed for fluctuating charges |
| 768 |
|
excluded_Pot= 0.0; |
| 768 |
– |
|
| 769 |
|
|
| 770 |
|
// some variables we'll need independent of electrostatic type: |
| 771 |
|
|
| 887 |
|
Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or |
| 888 |
|
+ rdQbr * rhat * (dv22 - 2.0*v22or)); |
| 889 |
|
} |
| 890 |
< |
|
| 890 |
> |
|
| 891 |
> |
|
| 892 |
|
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
| 893 |
|
J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; |
| 894 |
|
} |
| 895 |
< |
|
| 895 |
> |
|
| 896 |
|
if (a_is_Charge) { |
| 897 |
|
|
| 898 |
|
if (b_is_Charge) { |
| 899 |
|
pref = pre11_ * *(idat.electroMult); |
| 900 |
|
U += C_a * C_b * pref * v01; |
| 901 |
|
F += C_a * C_b * pref * dv01 * rhat; |
| 902 |
< |
|
| 902 |
> |
|
| 903 |
|
// If this is an excluded pair, there are still indirect |
| 904 |
|
// interactions via the reaction field we must worry about: |
| 905 |
|
|
| 908 |
|
indirect_Pot += rfContrib; |
| 909 |
|
indirect_F += rfContrib * 2.0 * ri * rhat; |
| 910 |
|
} |
| 911 |
< |
|
| 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. |
| 915 |
|
|
| 916 |
< |
if (idat.excluded) { |
| 916 |
> |
if (idat.excluded) { |
| 917 |
|
if (a_is_Fluctuating || b_is_Fluctuating) { |
| 918 |
|
coulInt = J->getValueAt( *(idat.rij) ); |
| 919 |
< |
if (a_is_Fluctuating) dUdCa += coulInt * C_b; |
| 920 |
< |
if (b_is_Fluctuating) dUdCb += coulInt * C_a; |
| 921 |
< |
excluded_Pot += C_a * C_b * coulInt; |
| 921 |
< |
} |
| 919 |
> |
if (a_is_Fluctuating) dUdCa += C_b * coulInt; |
| 920 |
> |
if (b_is_Fluctuating) dUdCb += C_a * coulInt; |
| 921 |
> |
} |
| 922 |
|
} else { |
| 923 |
|
if (a_is_Fluctuating) dUdCa += C_b * pref * v01; |
| 924 |
< |
if (a_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 925 |
< |
} |
| 924 |
> |
if (b_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 925 |
> |
} |
| 926 |
|
} |
| 927 |
|
|
| 928 |
|
if (b_is_Dipole) { |
| 988 |
|
F -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat; |
| 989 |
|
Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa); |
| 990 |
|
Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb); |
| 991 |
– |
|
| 991 |
|
// Even if we excluded this pair from direct interactions, we |
| 992 |
|
// still have the reaction-field-mediated dipole-dipole |
| 993 |
|
// interaction: |
| 1047 |
|
trQaQb = QaQb.trace(); |
| 1048 |
|
rQaQb = rhat * QaQb; |
| 1049 |
|
QaQbr = QaQb * rhat; |
| 1050 |
< |
QaxQb = cross(Q_a, Q_b); |
| 1050 |
> |
QaxQb = mCross(Q_a, Q_b); |
| 1051 |
|
rQaQbr = dot(rQa, Qbr); |
| 1052 |
|
rQaxQbr = cross(rQa, Qbr); |
| 1053 |
|
|
| 1078 |
|
// + 4.0 * cross(rhat, QbQar) |
| 1079 |
|
|
| 1080 |
|
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1082 |
– |
|
| 1081 |
|
} |
| 1082 |
|
} |
| 1083 |
|
|
| 1140 |
|
|
| 1141 |
|
if (i_is_Fluctuating) { |
| 1142 |
|
C_a += *(sdat.flucQ); |
| 1143 |
< |
// dVdFQ is really a force, so this is negative the derivative |
| 1146 |
< |
*(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1143 |
> |
*(sdat.flucQfrc) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1144 |
|
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * |
| 1145 |
|
(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); |
| 1146 |
|
} |
| 1438 |
|
dks[i] = dk * skr[i]; |
| 1439 |
|
} |
| 1440 |
|
if (data.is_Quadrupole) { |
| 1441 |
< |
Q = atom->getQuadrupole(); |
| 1442 |
< |
Q *= mPoleConverter; |
| 1446 |
< |
Qk = Q * kVec; |
| 1441 |
> |
Q = atom->getQuadrupole() * mPoleConverter; |
| 1442 |
> |
Qk = Q * kVec; |
| 1443 |
|
qk = dot(kVec, Qk); |
| 1444 |
< |
qxk[i] = cross(kVec, Qk); |
| 1444 |
> |
qxk[i] = -cross(kVec, Qk); |
| 1445 |
|
qkc[i] = qk * ckr[i]; |
| 1446 |
|
qks[i] = qk * skr[i]; |
| 1447 |
|
} |
| 1500 |
|
RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) |
| 1501 |
|
-ckr[i]*(ckss+dkcs-qkss)); |
| 1502 |
|
RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) |
| 1503 |
< |
+skr[i]*(ckss+dkcs-qkss)); |
| 1503 |
> |
+skr[i]*(ckss+dkcs-qkss)); |
| 1504 |
|
|
| 1505 |
|
atom->addFrc( 4.0 * rvol * qfrc * kVec ); |
| 1506 |
|
|
