--- trunk/src/nonbonded/Electrostatic.cpp 2013/07/29 17:55:17 1915 +++ trunk/src/nonbonded/Electrostatic.cpp 2013/08/05 16:13:46 1923 @@ -677,7 +677,8 @@ namespace OpenMD { FQtids[atid] = fqtid; Jij[fqtid].resize(nFlucq_); - // Now, iterate over all known fluctuating and add to the coulomb integral map: + // Now, iterate over all known fluctuating and add to the + // coulomb integral map: std::set::iterator it; for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { @@ -1165,6 +1166,7 @@ namespace OpenMD { case esm_SHIFTED_FORCE: case esm_SHIFTED_POTENTIAL: case esm_TAYLOR_SHIFTED: + case esm_EWALD_FULL: if (i_is_Charge) self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); if (i_is_Dipole) @@ -1192,7 +1194,7 @@ namespace OpenMD { } - void Electrostatic::ReciprocalSpaceSum () { + void Electrostatic::ReciprocalSpaceSum(potVec& pot) { RealType kPot = 0.0; RealType kVir = 0.0; @@ -1217,8 +1219,8 @@ namespace OpenMD { Vector3d box = hmat.diagonals(); RealType boxMax = box.max(); - //int kMax = int(pow(dampingAlpha_,2)*cutoffRadius_ * boxMax / M_PI); - const int kMax = 5; + //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); + int kMax = 7; int kSqMax = kMax*kMax + 2; int kLimit = kMax+1; @@ -1252,6 +1254,7 @@ namespace OpenMD { Vector3d t( 2.0 * M_PI ); t.Vdiv(t, box); + SimInfo::MoleculeIterator mi; Molecule::AtomIterator ai; @@ -1260,6 +1263,8 @@ namespace OpenMD { Vector3d tt; Vector3d w; Vector3d u; + Vector3d a; + Vector3d b; for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { @@ -1271,20 +1276,33 @@ namespace OpenMD { info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); tt.Vmul(t, r); + eCos[1][i] = Vector3d(1.0, 1.0, 1.0); eSin[1][i] = Vector3d(0.0, 0.0, 0.0); eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z())); eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z())); - u = 2.0 * eCos[1][i]; - eCos[3][i].Vmul(u, eCos[2][i]); - eSin[3][i].Vmul(u, eSin[2][i]); + + u = eCos[2][i]; + w = eSin[2][i]; for(int l = 3; l <= kLimit; l++) { - w.Vmul(u, eCos[l-1][i]); - eCos[l][i] = w - eCos[l-2][i]; - w.Vmul(u, eSin[l-1][i]); - eSin[l][i] = w - eSin[l-2][i]; + eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x(); + eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y(); + eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z(); + + eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x(); + eSin[l][i].y() = eSin[l-1][i].y()*eCos[2][i].y() + eCos[l-1][i].y()*eSin[2][i].y(); + eSin[l][i].z() = eSin[l-1][i].z()*eCos[2][i].z() + eCos[l-1][i].z()*eSin[2][i].z(); + + + // a.Vmul(eCos[l-1][i], u); + // b.Vmul(eSin[l-1][i], w); + // eCos[l][i] = a - b; + // a.Vmul(eSin[l-1][i], u); + // b.Vmul(eCos[l-1][i], w); + // eSin[l][i] = a + b; + } } } @@ -1332,7 +1350,7 @@ namespace OpenMD { int mMin = kLimit; int nMin = kLimit + 1; for (int l = 1; l <= kLimit; l++) { - int ll =l - 1; + int ll = l - 1; RealType rl = xcl * float(ll); for (int mmm = mMin; mmm <= kLim2; mmm++) { int mm = mmm - kLimit; @@ -1348,13 +1366,13 @@ namespace OpenMD { if(mm < 0) { clm[i] = eCos[l][i].x()*eCos[m][i].y() + eSin[l][i].x()*eSin[m][i].y(); - slm[i] = eCos[l][i].x()*eCos[m][i].y() + slm[i] = eSin[l][i].x()*eCos[m][i].y() - eSin[m][i].y()*eCos[l][i].x(); } else { clm[i] = eCos[l][i].x()*eCos[m][i].y() - eSin[l][i].x()*eSin[m][i].y(); slm[i] = eSin[l][i].x()*eCos[m][i].y() - + eSin[m][i].y()*eCos[l][i].x(); + + eSin[m][i].y()*eCos[l][i].x(); } } } @@ -1390,7 +1408,7 @@ namespace OpenMD { mol = info_->nextMolecule(mi)) { for(Atom* atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - i = atom->getGlobalIndex(); + i = atom->getLocalIndex(); int atid = atom->getAtomType()->getIdent(); ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; @@ -1398,7 +1416,7 @@ namespace OpenMD { RealType C = data.fixedCharge; if (atom->isFluctuatingCharge()) C += atom->getFlucQPos(); ckc[i] = C * ckr[i]; - cks[i] = C * cks[i]; + cks[i] = C * skr[i]; } if (data.is_Dipole) { @@ -1418,7 +1436,7 @@ namespace OpenMD { } } } - + // calculate vector sums RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); @@ -1427,6 +1445,7 @@ namespace OpenMD { RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0); RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0); + #ifdef IS_MPI MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, @@ -1444,13 +1463,10 @@ namespace OpenMD { #endif // Accumulate potential energy and virial contribution: - - //cerr << "l, m, n = " << l << " " << m << " " << n << "\n"; - cerr << "kVec = " << kVec << "\n"; - cerr << "ckss = " << ckss << " ckcs = " << ckcs << "\n"; + kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss)); - //cerr << "kspace pot = " << kPot << "\n"; + kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss +4.0*(ckss*dkcs-ckcs*dkss) +3.0*(dkcs*dkcs+dkss*dkss) @@ -1470,13 +1486,12 @@ namespace OpenMD { ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs) - - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); + - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) -ckr[i]*(ckss+dkcs-qkss)); RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+ skr[i]*(ckss+dkcs-qkss)); - - + atom->addFrc( 4.0 * rvol * qfrc * kVec ); if (data.is_Dipole) { @@ -1489,7 +1504,11 @@ namespace OpenMD { } } } + nMin = 1; } + mMin = 1; } + cerr << "kPot = " << kPot << "\n"; + pot[ELECTROSTATIC_FAMILY] += kPot; } }