| 57 |
|
#include "math/erfc.hpp" |
| 58 |
|
#include "math/SquareMatrix.hpp" |
| 59 |
|
#include "primitives/Molecule.hpp" |
| 60 |
+ |
#ifdef IS_MPI |
| 61 |
+ |
#include <mpi.h> |
| 62 |
+ |
#endif |
| 63 |
|
|
| 61 |
– |
|
| 64 |
|
namespace OpenMD { |
| 65 |
|
|
| 66 |
|
Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), |
| 679 |
|
FQtids[atid] = fqtid; |
| 680 |
|
Jij[fqtid].resize(nFlucq_); |
| 681 |
|
|
| 682 |
< |
// Now, iterate over all known fluctuating and add to the coulomb integral map: |
| 682 |
> |
// Now, iterate over all known fluctuating and add to the |
| 683 |
> |
// coulomb integral map: |
| 684 |
|
|
| 685 |
|
std::set<int>::iterator it; |
| 686 |
|
for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { |
| 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 |
| 1143 |
< |
*(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 |
|
} |
| 1170 |
|
case esm_SHIFTED_FORCE: |
| 1171 |
|
case esm_SHIFTED_POTENTIAL: |
| 1172 |
|
case esm_TAYLOR_SHIFTED: |
| 1173 |
+ |
case esm_EWALD_FULL: |
| 1174 |
|
if (i_is_Charge) |
| 1175 |
|
self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); |
| 1176 |
|
if (i_is_Dipole) |
| 1198 |
|
} |
| 1199 |
|
|
| 1200 |
|
|
| 1201 |
< |
void Electrostatic::ReciprocalSpaceSum () { |
| 1201 |
> |
void Electrostatic::ReciprocalSpaceSum(RealType& pot) { |
| 1202 |
|
|
| 1203 |
|
RealType kPot = 0.0; |
| 1204 |
|
RealType kVir = 0.0; |
| 1223 |
|
Vector3d box = hmat.diagonals(); |
| 1224 |
|
RealType boxMax = box.max(); |
| 1225 |
|
|
| 1226 |
< |
//int kMax = int(pow(dampingAlpha_,2)*cutoffRadius_ * boxMax / M_PI); |
| 1227 |
< |
const int kMax = 5; |
| 1226 |
> |
//int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); |
| 1227 |
> |
int kMax = 7; |
| 1228 |
|
int kSqMax = kMax*kMax + 2; |
| 1229 |
|
|
| 1230 |
|
int kLimit = kMax+1; |
| 1244 |
|
|
| 1245 |
|
// Calculate and store exponential factors |
| 1246 |
|
|
| 1247 |
< |
vector<vector<Vector3d> > eCos; |
| 1248 |
< |
vector<vector<Vector3d> > eSin; |
| 1247 |
> |
vector<vector<RealType> > elc; |
| 1248 |
> |
vector<vector<RealType> > emc; |
| 1249 |
> |
vector<vector<RealType> > enc; |
| 1250 |
> |
vector<vector<RealType> > els; |
| 1251 |
> |
vector<vector<RealType> > ems; |
| 1252 |
> |
vector<vector<RealType> > ens; |
| 1253 |
> |
|
| 1254 |
|
|
| 1255 |
|
int nMax = info_->getNAtoms(); |
| 1256 |
|
|
| 1257 |
< |
eCos.resize(kLimit+1); |
| 1258 |
< |
eSin.resize(kLimit+1); |
| 1257 |
> |
elc.resize(kLimit+1); |
| 1258 |
> |
emc.resize(kLimit+1); |
| 1259 |
> |
enc.resize(kLimit+1); |
| 1260 |
> |
els.resize(kLimit+1); |
| 1261 |
> |
ems.resize(kLimit+1); |
| 1262 |
> |
ens.resize(kLimit+1); |
| 1263 |
> |
|
| 1264 |
|
for (int j = 0; j < kLimit+1; j++) { |
| 1265 |
< |
eCos[j].resize(nMax); |
| 1266 |
< |
eSin[j].resize(nMax); |
| 1265 |
> |
elc[j].resize(nMax); |
| 1266 |
> |
emc[j].resize(nMax); |
| 1267 |
> |
enc[j].resize(nMax); |
| 1268 |
> |
els[j].resize(nMax); |
| 1269 |
> |
ems[j].resize(nMax); |
| 1270 |
> |
ens[j].resize(nMax); |
| 1271 |
|
} |
| 1272 |
|
|
| 1273 |
|
Vector3d t( 2.0 * M_PI ); |
| 1274 |
|
t.Vdiv(t, box); |
| 1275 |
+ |
|
| 1276 |
|
|
| 1277 |
|
SimInfo::MoleculeIterator mi; |
| 1278 |
|
Molecule::AtomIterator ai; |
| 1279 |
|
int i; |
| 1280 |
|
Vector3d r; |
| 1281 |
|
Vector3d tt; |
| 1261 |
– |
Vector3d w; |
| 1262 |
– |
Vector3d u; |
| 1282 |
|
|
| 1283 |
|
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1284 |
|
mol = info_->nextMolecule(mi)) { |
| 1290 |
|
info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); |
| 1291 |
|
|
| 1292 |
|
tt.Vmul(t, r); |
| 1293 |
+ |
|
| 1294 |
+ |
elc[1][i] = 1.0; |
| 1295 |
+ |
emc[1][i] = 1.0; |
| 1296 |
+ |
enc[1][i] = 1.0; |
| 1297 |
+ |
els[1][i] = 0.0; |
| 1298 |
+ |
ems[1][i] = 0.0; |
| 1299 |
+ |
ens[1][i] = 0.0; |
| 1300 |
+ |
|
| 1301 |
+ |
elc[2][i] = cos(tt.x()); |
| 1302 |
+ |
emc[2][i] = cos(tt.y()); |
| 1303 |
+ |
enc[2][i] = cos(tt.z()); |
| 1304 |
+ |
els[2][i] = sin(tt.x()); |
| 1305 |
+ |
ems[2][i] = sin(tt.y()); |
| 1306 |
+ |
ens[2][i] = sin(tt.z()); |
| 1307 |
|
|
| 1275 |
– |
eCos[1][i] = Vector3d(1.0, 1.0, 1.0); |
| 1276 |
– |
eSin[1][i] = Vector3d(0.0, 0.0, 0.0); |
| 1277 |
– |
eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z())); |
| 1278 |
– |
eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z())); |
| 1279 |
– |
u = 2.0 * eCos[1][i]; |
| 1280 |
– |
eCos[3][i].Vmul(u, eCos[2][i]); |
| 1281 |
– |
eSin[3][i].Vmul(u, eSin[2][i]); |
| 1282 |
– |
|
| 1308 |
|
for(int l = 3; l <= kLimit; l++) { |
| 1309 |
< |
w.Vmul(u, eCos[l-1][i]); |
| 1310 |
< |
eCos[l][i] = w - eCos[l-2][i]; |
| 1311 |
< |
w.Vmul(u, eSin[l-1][i]); |
| 1312 |
< |
eSin[l][i] = w - eSin[l-2][i]; |
| 1309 |
> |
elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i]; |
| 1310 |
> |
emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i]; |
| 1311 |
> |
enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i]; |
| 1312 |
> |
els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i]; |
| 1313 |
> |
ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i]; |
| 1314 |
> |
ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i]; |
| 1315 |
|
} |
| 1316 |
|
} |
| 1317 |
|
} |
| 1355 |
|
std::vector<RealType> qks(nMax, 0.0); |
| 1356 |
|
std::vector<Vector3d> dxk(nMax, V3Zero); |
| 1357 |
|
std::vector<Vector3d> qxk(nMax, V3Zero); |
| 1358 |
< |
|
| 1358 |
> |
RealType rl, rm, rn; |
| 1359 |
> |
Vector3d kVec; |
| 1360 |
> |
Vector3d Qk; |
| 1361 |
> |
Mat3x3d k2; |
| 1362 |
> |
RealType ckcs, ckss, dkcs, dkss, qkcs, qkss; |
| 1363 |
> |
int atid; |
| 1364 |
> |
ElectrostaticAtomData data; |
| 1365 |
> |
RealType C, dk, qk; |
| 1366 |
> |
Vector3d D; |
| 1367 |
> |
Mat3x3d Q; |
| 1368 |
> |
|
| 1369 |
|
int mMin = kLimit; |
| 1370 |
|
int nMin = kLimit + 1; |
| 1371 |
|
for (int l = 1; l <= kLimit; l++) { |
| 1372 |
< |
int ll =l - 1; |
| 1373 |
< |
RealType rl = xcl * float(ll); |
| 1372 |
> |
int ll = l - 1; |
| 1373 |
> |
rl = xcl * float(ll); |
| 1374 |
|
for (int mmm = mMin; mmm <= kLim2; mmm++) { |
| 1375 |
|
int mm = mmm - kLimit; |
| 1376 |
|
int m = abs(mm) + 1; |
| 1377 |
< |
RealType rm = ycl * float(mm); |
| 1377 |
> |
rm = ycl * float(mm); |
| 1378 |
|
// Set temporary products of exponential terms |
| 1379 |
|
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1380 |
|
mol = info_->nextMolecule(mi)) { |
| 1383 |
|
|
| 1384 |
|
i = atom->getLocalIndex(); |
| 1385 |
|
if(mm < 0) { |
| 1386 |
< |
clm[i] = eCos[l][i].x()*eCos[m][i].y() |
| 1387 |
< |
+ eSin[l][i].x()*eSin[m][i].y(); |
| 1351 |
< |
slm[i] = eCos[l][i].x()*eCos[m][i].y() |
| 1352 |
< |
- eSin[m][i].y()*eCos[l][i].x(); |
| 1386 |
> |
clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i]; |
| 1387 |
> |
slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i]; |
| 1388 |
|
} else { |
| 1389 |
< |
clm[i] = eCos[l][i].x()*eCos[m][i].y() |
| 1390 |
< |
- eSin[l][i].x()*eSin[m][i].y(); |
| 1356 |
< |
slm[i] = eSin[l][i].x()*eCos[m][i].y() |
| 1357 |
< |
+ eSin[m][i].y()*eCos[l][i].x(); |
| 1389 |
> |
clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i]; |
| 1390 |
> |
slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i]; |
| 1391 |
|
} |
| 1392 |
|
} |
| 1393 |
|
} |
| 1394 |
|
for (int nnn = nMin; nnn <= kLim2; nnn++) { |
| 1395 |
|
int nn = nnn - kLimit; |
| 1396 |
|
int n = abs(nn) + 1; |
| 1397 |
< |
RealType rn = zcl * float(nn); |
| 1397 |
> |
rn = zcl * float(nn); |
| 1398 |
|
// Test on magnitude of k vector: |
| 1399 |
|
int kk=ll*ll + mm*mm + nn*nn; |
| 1400 |
|
if(kk <= kSqLim) { |
| 1401 |
< |
Vector3d kVec = Vector3d(rl, rm, rn); |
| 1402 |
< |
Mat3x3d k2 = outProduct(kVec, kVec); |
| 1401 |
> |
kVec = Vector3d(rl, rm, rn); |
| 1402 |
> |
k2 = outProduct(kVec, kVec); |
| 1403 |
|
// Calculate exp(ikr) terms |
| 1404 |
|
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1405 |
|
mol = info_->nextMolecule(mi)) { |
| 1408 |
|
i = atom->getLocalIndex(); |
| 1409 |
|
|
| 1410 |
|
if (nn < 0) { |
| 1411 |
< |
ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z(); |
| 1412 |
< |
skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z(); |
| 1411 |
> |
ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i]; |
| 1412 |
> |
skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i]; |
| 1413 |
> |
|
| 1414 |
|
} else { |
| 1415 |
< |
ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z(); |
| 1416 |
< |
skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z(); |
| 1415 |
> |
ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i]; |
| 1416 |
> |
skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i]; |
| 1417 |
|
} |
| 1418 |
|
} |
| 1419 |
|
} |
| 1424 |
|
mol = info_->nextMolecule(mi)) { |
| 1425 |
|
for(Atom* atom = mol->beginAtom(ai); atom != NULL; |
| 1426 |
|
atom = mol->nextAtom(ai)) { |
| 1427 |
< |
i = atom->getGlobalIndex(); |
| 1427 |
> |
i = atom->getLocalIndex(); |
| 1428 |
|
int atid = atom->getAtomType()->getIdent(); |
| 1429 |
< |
ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; |
| 1429 |
> |
data = ElectrostaticMap[Etids[atid]]; |
| 1430 |
|
|
| 1431 |
|
if (data.is_Charge) { |
| 1432 |
< |
RealType C = data.fixedCharge; |
| 1432 |
> |
C = data.fixedCharge; |
| 1433 |
|
if (atom->isFluctuatingCharge()) C += atom->getFlucQPos(); |
| 1434 |
|
ckc[i] = C * ckr[i]; |
| 1435 |
< |
cks[i] = C * cks[i]; |
| 1435 |
> |
cks[i] = C * skr[i]; |
| 1436 |
|
} |
| 1437 |
|
|
| 1438 |
|
if (data.is_Dipole) { |
| 1439 |
< |
Vector3d D = atom->getDipole() * mPoleConverter; |
| 1440 |
< |
RealType dk = dot(kVec, D); |
| 1441 |
< |
dxk[i] = cross(kVec, D); |
| 1439 |
> |
D = atom->getDipole() * mPoleConverter; |
| 1440 |
> |
dk = dot(D, kVec); |
| 1441 |
> |
dxk[i] = cross(D, kVec); |
| 1442 |
|
dkc[i] = dk * ckr[i]; |
| 1443 |
|
dks[i] = dk * skr[i]; |
| 1444 |
|
} |
| 1445 |
|
if (data.is_Quadrupole) { |
| 1446 |
< |
Mat3x3d Q = atom->getQuadrupole(); |
| 1446 |
> |
Q = atom->getQuadrupole(); |
| 1447 |
|
Q *= mPoleConverter; |
| 1448 |
< |
RealType qk = -( Q * k2 ).trace(); |
| 1449 |
< |
qxk[i] = -2.0 * cross(k2, Q); |
| 1448 |
> |
Qk = Q * kVec; |
| 1449 |
> |
qk = dot(kVec, Qk); |
| 1450 |
> |
qxk[i] = cross(kVec, Qk); |
| 1451 |
|
qkc[i] = qk * ckr[i]; |
| 1452 |
|
qks[i] = qk * skr[i]; |
| 1453 |
|
} |
| 1454 |
|
} |
| 1455 |
|
} |
| 1456 |
< |
|
| 1456 |
> |
|
| 1457 |
|
// calculate vector sums |
| 1458 |
|
|
| 1459 |
< |
RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); |
| 1460 |
< |
RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0); |
| 1461 |
< |
RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0); |
| 1462 |
< |
RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0); |
| 1463 |
< |
RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); |
| 1464 |
< |
RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0); |
| 1459 |
> |
ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); |
| 1460 |
> |
ckss = std::accumulate(cks.begin(),cks.end(),0.0); |
| 1461 |
> |
dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0); |
| 1462 |
> |
dkss = std::accumulate(dks.begin(),dks.end(),0.0); |
| 1463 |
> |
qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); |
| 1464 |
> |
qkss = std::accumulate(qks.begin(),qks.end(),0.0); |
| 1465 |
|
|
| 1466 |
|
#ifdef IS_MPI |
| 1467 |
|
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, |
| 1479 |
|
#endif |
| 1480 |
|
|
| 1481 |
|
// Accumulate potential energy and virial contribution: |
| 1482 |
+ |
|
| 1483 |
+ |
kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) |
| 1484 |
+ |
+ (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs)); |
| 1485 |
+ |
|
| 1486 |
+ |
kVir += 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss |
| 1487 |
+ |
+4.0*(ckss*dkcs-ckcs*dkss) |
| 1488 |
+ |
+3.0*(dkcs*dkcs+dkss*dkss) |
| 1489 |
+ |
-6.0*(ckss*qkss+ckcs*qkcs) |
| 1490 |
+ |
+8.0*(dkss*qkcs-dkcs*qkss) |
| 1491 |
+ |
+5.0*(qkss*qkss+qkcs*qkcs)); |
| 1492 |
|
|
| 1448 |
– |
//cerr << "l, m, n = " << l << " " << m << " " << n << "\n"; |
| 1449 |
– |
cerr << "kVec = " << kVec << "\n"; |
| 1450 |
– |
cerr << "ckss = " << ckss << " ckcs = " << ckcs << "\n"; |
| 1451 |
– |
kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) |
| 1452 |
– |
+ (ckcs-dkss-qkcs)*(ckcs-dkss-qkss)); |
| 1453 |
– |
//cerr << "kspace pot = " << kPot << "\n"; |
| 1454 |
– |
kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss |
| 1455 |
– |
+4.0*(ckss*dkcs-ckcs*dkss) |
| 1456 |
– |
+3.0*(dkcs*dkcs+dkss*dkss) |
| 1457 |
– |
-6.0*(ckss*qkss+ckcs*qkcs) |
| 1458 |
– |
+8.0*(dkss*qkcs-dkcs*qkss) |
| 1459 |
– |
+5.0*(qkss*qkss+qkcs*qkcs)); |
| 1460 |
– |
|
| 1493 |
|
// Calculate force and torque for each site: |
| 1494 |
|
|
| 1495 |
|
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1498 |
|
atom = mol->nextAtom(ai)) { |
| 1499 |
|
|
| 1500 |
|
i = atom->getLocalIndex(); |
| 1501 |
< |
int atid = atom->getAtomType()->getIdent(); |
| 1502 |
< |
ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; |
| 1503 |
< |
|
| 1501 |
> |
atid = atom->getAtomType()->getIdent(); |
| 1502 |
> |
data = ElectrostaticMap[Etids[atid]]; |
| 1503 |
> |
|
| 1504 |
|
RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs) |
| 1505 |
< |
- (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); |
| 1505 |
> |
- (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); |
| 1506 |
|
RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) |
| 1507 |
|
-ckr[i]*(ckss+dkcs-qkss)); |
| 1508 |
< |
RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+ |
| 1509 |
< |
skr[i]*(ckss+dkcs-qkss)); |
| 1510 |
< |
|
| 1479 |
< |
|
| 1508 |
> |
RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) |
| 1509 |
> |
+skr[i]*(ckss+dkcs-qkss)); |
| 1510 |
> |
|
| 1511 |
|
atom->addFrc( 4.0 * rvol * qfrc * kVec ); |
| 1512 |
|
|
| 1513 |
|
if (data.is_Dipole) { |
| 1520 |
|
} |
| 1521 |
|
} |
| 1522 |
|
} |
| 1523 |
+ |
nMin = 1; |
| 1524 |
|
} |
| 1525 |
+ |
mMin = 1; |
| 1526 |
|
} |
| 1527 |
+ |
pot += kPot; |
| 1528 |
|
} |
| 1529 |
|
} |
