| 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" |
| 64 |
< |
#ifdef IS_MPI |
| 61 |
< |
#include <mpi.h> |
| 62 |
< |
#endif |
| 64 |
> |
#include "flucq/FluctuatingChargeForces.hpp" |
| 65 |
|
|
| 66 |
|
namespace OpenMD { |
| 67 |
|
|
| 71 |
|
haveDampingAlpha_(false), |
| 72 |
|
haveDielectric_(false), |
| 73 |
|
haveElectroSplines_(false) |
| 74 |
< |
{} |
| 74 |
> |
{ |
| 75 |
> |
flucQ_ = new FluctuatingChargeForces(info_); |
| 76 |
> |
} |
| 77 |
|
|
| 78 |
+ |
void Electrostatic::setForceField(ForceField *ff) { |
| 79 |
+ |
forceField_ = ff; |
| 80 |
+ |
flucQ_->setForceField(forceField_); |
| 81 |
+ |
} |
| 82 |
+ |
|
| 83 |
+ |
void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) { |
| 84 |
+ |
simTypes_ = simtypes; |
| 85 |
+ |
flucQ_->setSimulatedAtomTypes(simTypes_); |
| 86 |
+ |
} |
| 87 |
+ |
|
| 88 |
|
void Electrostatic::initialize() { |
| 89 |
|
|
| 90 |
|
Globals* simParams_ = info_->getSimParams(); |
| 780 |
|
// Excluded potential that is still computed for fluctuating charges |
| 781 |
|
excluded_Pot= 0.0; |
| 782 |
|
|
| 769 |
– |
|
| 783 |
|
// some variables we'll need independent of electrostatic type: |
| 784 |
|
|
| 785 |
|
ri = 1.0 / *(idat.rij); |
| 900 |
|
Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or |
| 901 |
|
+ rdQbr * rhat * (dv22 - 2.0*v22or)); |
| 902 |
|
} |
| 903 |
< |
|
| 903 |
> |
|
| 904 |
> |
|
| 905 |
|
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
| 906 |
|
J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; |
| 907 |
|
} |
| 908 |
< |
|
| 908 |
> |
|
| 909 |
|
if (a_is_Charge) { |
| 910 |
|
|
| 911 |
|
if (b_is_Charge) { |
| 912 |
|
pref = pre11_ * *(idat.electroMult); |
| 913 |
|
U += C_a * C_b * pref * v01; |
| 914 |
|
F += C_a * C_b * pref * dv01 * rhat; |
| 915 |
< |
|
| 915 |
> |
|
| 916 |
|
// If this is an excluded pair, there are still indirect |
| 917 |
|
// interactions via the reaction field we must worry about: |
| 918 |
|
|
| 921 |
|
indirect_Pot += rfContrib; |
| 922 |
|
indirect_F += rfContrib * 2.0 * ri * rhat; |
| 923 |
|
} |
| 924 |
< |
|
| 924 |
> |
|
| 925 |
|
// Fluctuating charge forces are handled via Coulomb integrals |
| 926 |
|
// for excluded pairs (i.e. those connected via bonds) and |
| 927 |
|
// with the standard charge-charge interaction otherwise. |
| 928 |
|
|
| 929 |
< |
if (idat.excluded) { |
| 929 |
> |
if (idat.excluded) { |
| 930 |
|
if (a_is_Fluctuating || b_is_Fluctuating) { |
| 931 |
|
coulInt = J->getValueAt( *(idat.rij) ); |
| 932 |
< |
if (a_is_Fluctuating) dUdCa += coulInt * C_b; |
| 933 |
< |
if (b_is_Fluctuating) dUdCb += coulInt * C_a; |
| 934 |
< |
excluded_Pot += C_a * C_b * coulInt; |
| 921 |
< |
} |
| 932 |
> |
if (a_is_Fluctuating) dUdCa += C_b * coulInt; |
| 933 |
> |
if (b_is_Fluctuating) dUdCb += C_a * coulInt; |
| 934 |
> |
} |
| 935 |
|
} else { |
| 936 |
|
if (a_is_Fluctuating) dUdCa += C_b * pref * v01; |
| 937 |
< |
if (a_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 938 |
< |
} |
| 937 |
> |
if (b_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 938 |
> |
} |
| 939 |
|
} |
| 940 |
|
|
| 941 |
|
if (b_is_Dipole) { |
| 1001 |
|
F -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat; |
| 1002 |
|
Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa); |
| 1003 |
|
Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb); |
| 991 |
– |
|
| 1004 |
|
// Even if we excluded this pair from direct interactions, we |
| 1005 |
|
// still have the reaction-field-mediated dipole-dipole |
| 1006 |
|
// interaction: |
| 1060 |
|
trQaQb = QaQb.trace(); |
| 1061 |
|
rQaQb = rhat * QaQb; |
| 1062 |
|
QaQbr = QaQb * rhat; |
| 1063 |
< |
QaxQb = cross(Q_a, Q_b); |
| 1063 |
> |
QaxQb = mCross(Q_a, Q_b); |
| 1064 |
|
rQaQbr = dot(rQa, Qbr); |
| 1065 |
|
rQaxQbr = cross(rQa, Qbr); |
| 1066 |
|
|
| 1091 |
|
// + 4.0 * cross(rhat, QbQar) |
| 1092 |
|
|
| 1093 |
|
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1082 |
– |
|
| 1094 |
|
} |
| 1095 |
|
} |
| 1096 |
|
|
| 1153 |
|
|
| 1154 |
|
if (i_is_Fluctuating) { |
| 1155 |
|
C_a += *(sdat.flucQ); |
| 1156 |
< |
// dVdFQ is really a force, so this is negative the derivative |
| 1157 |
< |
*(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1158 |
< |
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * |
| 1159 |
< |
(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); |
| 1156 |
> |
|
| 1157 |
> |
flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ), |
| 1158 |
> |
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY], |
| 1159 |
> |
*(sdat.flucQfrc) ); |
| 1160 |
> |
|
| 1161 |
|
} |
| 1162 |
|
|
| 1163 |
|
switch (summationMethod_) { |
| 1260 |
|
vector<vector<RealType> > els; |
| 1261 |
|
vector<vector<RealType> > ems; |
| 1262 |
|
vector<vector<RealType> > ens; |
| 1251 |
– |
|
| 1263 |
|
|
| 1264 |
|
int nMax = info_->getNAtoms(); |
| 1265 |
|
|
| 1282 |
|
Vector3d t( 2.0 * M_PI ); |
| 1283 |
|
t.Vdiv(t, box); |
| 1284 |
|
|
| 1274 |
– |
|
| 1285 |
|
SimInfo::MoleculeIterator mi; |
| 1286 |
|
Molecule::AtomIterator ai; |
| 1287 |
|
int i; |
| 1451 |
|
dks[i] = dk * skr[i]; |
| 1452 |
|
} |
| 1453 |
|
if (data.is_Quadrupole) { |
| 1454 |
< |
Q = atom->getQuadrupole(); |
| 1455 |
< |
Q *= mPoleConverter; |
| 1446 |
< |
Qk = Q * kVec; |
| 1454 |
> |
Q = atom->getQuadrupole() * mPoleConverter; |
| 1455 |
> |
Qk = Q * kVec; |
| 1456 |
|
qk = dot(kVec, Qk); |
| 1457 |
< |
qxk[i] = cross(kVec, Qk); |
| 1457 |
> |
qxk[i] = -cross(kVec, Qk); |
| 1458 |
|
qkc[i] = qk * ckr[i]; |
| 1459 |
|
qks[i] = qk * skr[i]; |
| 1460 |
|
} |
| 1471 |
|
qkss = std::accumulate(qks.begin(),qks.end(),0.0); |
| 1472 |
|
|
| 1473 |
|
#ifdef IS_MPI |
| 1474 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, |
| 1475 |
< |
MPI::SUM); |
| 1476 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE, |
| 1477 |
< |
MPI::SUM); |
| 1478 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE, |
| 1479 |
< |
MPI::SUM); |
| 1480 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE, |
| 1481 |
< |
MPI::SUM); |
| 1482 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE, |
| 1483 |
< |
MPI::SUM); |
| 1484 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE, |
| 1485 |
< |
MPI::SUM); |
| 1474 |
> |
MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE, |
| 1475 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 1476 |
> |
MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE, |
| 1477 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 1478 |
> |
MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE, |
| 1479 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 1480 |
> |
MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE, |
| 1481 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 1482 |
> |
MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE, |
| 1483 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 1484 |
> |
MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE, |
| 1485 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 1486 |
|
#endif |
| 1487 |
|
|
| 1488 |
|
// Accumulate potential energy and virial contribution: |
| 1513 |
|
RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) |
| 1514 |
|
-ckr[i]*(ckss+dkcs-qkss)); |
| 1515 |
|
RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) |
| 1516 |
< |
+skr[i]*(ckss+dkcs-qkss)); |
| 1516 |
> |
+skr[i]*(ckss+dkcs-qkss)); |
| 1517 |
|
|
| 1518 |
|
atom->addFrc( 4.0 * rvol * qfrc * kVec ); |
| 1519 |
< |
|
| 1519 |
> |
|
| 1520 |
> |
if (atom->isFluctuatingCharge()) { |
| 1521 |
> |
atom->addFlucQFrc( - 2.0 * rvol * qtrq2 ); |
| 1522 |
> |
} |
| 1523 |
> |
|
| 1524 |
|
if (data.is_Dipole) { |
| 1525 |
|
atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] ); |
| 1526 |
|
} |
