| 99 |
|
* Use the maximum suggested value that was found. |
| 100 |
|
* |
| 101 |
|
* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED, |
| 102 |
< |
* or SHIFTED_POTENTIAL) |
| 102 |
> |
* SHIFTED_POTENTIAL, or EWALD_FULL) |
| 103 |
|
* If cutoffMethod was explicitly set, use that choice. |
| 104 |
|
* If cutoffMethod was not explicitly set, use SHIFTED_FORCE |
| 105 |
|
* |
| 172 |
|
stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; |
| 173 |
|
stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; |
| 174 |
|
stringToCutoffMethod["TAYLOR_SHIFTED"] = TAYLOR_SHIFTED; |
| 175 |
+ |
stringToCutoffMethod["EWALD_FULL"] = EWALD_FULL; |
| 176 |
|
|
| 177 |
|
if (simParams_->haveCutoffMethod()) { |
| 178 |
|
string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); |
| 183 |
|
"ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" |
| 184 |
|
"\tShould be one of: " |
| 185 |
|
"HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n" |
| 186 |
< |
"\tor SHIFTED_FORCE\n", |
| 186 |
> |
"\tSHIFTED_FORCE, or EWALD_FULL\n", |
| 187 |
|
cutMeth.c_str()); |
| 188 |
|
painCave.isFatal = 1; |
| 189 |
|
painCave.severity = OPENMD_ERROR; |
| 229 |
|
cutoffMethod_ = SHIFTED_FORCE; |
| 230 |
|
} else if (myMethod == "TAYLOR_SHIFTED") { |
| 231 |
|
cutoffMethod_ = TAYLOR_SHIFTED; |
| 232 |
+ |
} else if (myMethod == "EWALD_FULL") { |
| 233 |
+ |
cutoffMethod_ = EWALD_FULL; |
| 234 |
|
} |
| 235 |
|
|
| 236 |
|
if (simParams_->haveSwitchingRadius()) |
| 237 |
|
rSwitch_ = simParams_->getSwitchingRadius(); |
| 238 |
|
|
| 239 |
|
if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE" || |
| 240 |
< |
myMethod == "TAYLOR_SHIFTED") { |
| 240 |
> |
myMethod == "TAYLOR_SHIFTED" || myMethod == "EWALD_FULL") { |
| 241 |
|
if (simParams_->haveSwitchingRadius()){ |
| 242 |
|
sprintf(painCave.errMsg, |
| 243 |
|
"ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n" |
| 698 |
|
RealType vij; |
| 699 |
|
Vector3d fij, fg, f1; |
| 700 |
|
tuple3<RealType, RealType, RealType> cuts; |
| 701 |
< |
RealType rCutSq; |
| 701 |
> |
RealType rCut, rCutSq, rListSq; |
| 702 |
|
bool in_switching_region; |
| 703 |
|
RealType sw, dswdr, swderiv; |
| 704 |
|
vector<int> atomListColumn, atomListRow; |
| 709 |
|
RealType dVdFQ1(0.0); |
| 710 |
|
RealType dVdFQ2(0.0); |
| 711 |
|
potVec longRangePotential(0.0); |
| 712 |
+ |
potVec reciprocalPotential(0.0); |
| 713 |
|
potVec workPot(0.0); |
| 714 |
|
potVec exPot(0.0); |
| 715 |
|
Vector3d eField1(0.0); |
| 717 |
|
vector<int>::iterator ia, jb; |
| 718 |
|
|
| 719 |
|
int loopStart, loopEnd; |
| 720 |
< |
|
| 720 |
> |
|
| 721 |
> |
idat.rcut = &rCut; |
| 722 |
|
idat.vdwMult = &vdwMult; |
| 723 |
|
idat.electroMult = &electroMult; |
| 724 |
|
idat.pot = &workPot; |
| 751 |
|
if (update_nlist) { |
| 752 |
|
if (!usePeriodicBoundaryConditions_) |
| 753 |
|
Mat3x3d bbox = thermo->getBoundingBox(); |
| 754 |
< |
neighborList = fDecomp_->buildNeighborList(); |
| 754 |
> |
fDecomp_->buildNeighborList(neighborList_); |
| 755 |
|
} |
| 756 |
|
} |
| 757 |
|
|
| 758 |
< |
for (vector<pair<int, int> >::iterator it = neighborList.begin(); |
| 759 |
< |
it != neighborList.end(); ++it) { |
| 758 |
> |
for (vector<pair<int, int> >::iterator it = neighborList_.begin(); |
| 759 |
> |
it != neighborList_.end(); ++it) { |
| 760 |
|
|
| 761 |
|
cg1 = (*it).first; |
| 762 |
|
cg2 = (*it).second; |
| 763 |
|
|
| 764 |
< |
cuts = fDecomp_->getGroupCutoffs(cg1, cg2); |
| 764 |
> |
fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq); |
| 765 |
|
|
| 766 |
|
d_grp = fDecomp_->getIntergroupVector(cg1, cg2); |
| 767 |
|
|
| 768 |
< |
curSnapshot->wrapVector(d_grp); |
| 768 |
> |
// already wrapped in the getIntergroupVector call: |
| 769 |
> |
// curSnapshot->wrapVector(d_grp); |
| 770 |
|
rgrpsq = d_grp.lengthSquare(); |
| 765 |
– |
rCutSq = cuts.second; |
| 771 |
|
|
| 772 |
|
if (rgrpsq < rCutSq) { |
| 773 |
< |
idat.rcut = &cuts.first; |
| 773 |
> |
|
| 774 |
|
if (iLoop == PAIR_LOOP) { |
| 775 |
|
vij = 0.0; |
| 776 |
|
fij.zero(); |
| 926 |
|
|
| 927 |
|
// collects pairwise information |
| 928 |
|
fDecomp_->collectData(); |
| 929 |
+ |
if (cutoffMethod_ == EWALD_FULL) { |
| 930 |
+ |
interactionMan_->doReciprocalSpaceSum(reciprocalPotential); |
| 931 |
+ |
} |
| 932 |
|
|
| 933 |
|
if (info_->requiresSelfCorrection()) { |
| 934 |
|
for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
| 941 |
|
fDecomp_->collectSelfData(); |
| 942 |
|
|
| 943 |
|
longRangePotential = *(fDecomp_->getEmbeddingPotential()) + |
| 944 |
< |
*(fDecomp_->getPairwisePotential()); |
| 944 |
> |
*(fDecomp_->getPairwisePotential()) + reciprocalPotential; |
| 945 |
|
|
| 946 |
|
curSnapshot->setLongRangePotential(longRangePotential); |
| 947 |
|
|
| 950 |
|
|
| 951 |
|
} |
| 952 |
|
|
| 945 |
– |
|
| 953 |
|
void ForceManager::postCalculation() { |
| 954 |
|
|
| 955 |
|
vector<Perturbation*>::iterator pi; |
