| 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" |
| 637 |
|
// Collect from all nodes. This should eventually be moved into a |
| 638 |
|
// SystemDecomposition, but this is a better place than in |
| 639 |
|
// Thermo to do the collection. |
| 640 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE, |
| 641 |
< |
MPI::SUM); |
| 642 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE, |
| 643 |
< |
MPI::SUM); |
| 644 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1, |
| 645 |
< |
MPI::REALTYPE, MPI::SUM); |
| 646 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1, |
| 647 |
< |
MPI::REALTYPE, MPI::SUM); |
| 640 |
> |
|
| 641 |
> |
MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE, |
| 642 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 643 |
> |
MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE, |
| 644 |
> |
MPI_SUM, MPI_COMM_WORLD); |
| 645 |
> |
MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1, |
| 646 |
> |
MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
| 647 |
> |
MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1, |
| 648 |
> |
MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
| 649 |
> |
// MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE, |
| 650 |
> |
// MPI::SUM); |
| 651 |
> |
// MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE, |
| 652 |
> |
// MPI::SUM); |
| 653 |
> |
// MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1, |
| 654 |
> |
// MPI::REALTYPE, MPI::SUM); |
| 655 |
> |
// MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1, |
| 656 |
> |
// MPI::REALTYPE, MPI::SUM); |
| 657 |
|
#endif |
| 658 |
|
|
| 659 |
|
Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 675 |
|
DataStorage* config = &(curSnapshot->atomData); |
| 676 |
|
DataStorage* cgConfig = &(curSnapshot->cgData); |
| 677 |
|
|
| 678 |
+ |
|
| 679 |
|
//calculate the center of mass of cutoff group |
| 680 |
|
|
| 681 |
|
SimInfo::MoleculeIterator mi; |
| 708 |
|
RealType vij; |
| 709 |
|
Vector3d fij, fg, f1; |
| 710 |
|
tuple3<RealType, RealType, RealType> cuts; |
| 711 |
< |
RealType rCutSq; |
| 711 |
> |
RealType rCut, rCutSq, rListSq; |
| 712 |
|
bool in_switching_region; |
| 713 |
|
RealType sw, dswdr, swderiv; |
| 714 |
|
vector<int> atomListColumn, atomListRow; |
| 718 |
|
RealType vpair; |
| 719 |
|
RealType dVdFQ1(0.0); |
| 720 |
|
RealType dVdFQ2(0.0); |
| 721 |
< |
Vector3d eField1(0.0); |
| 722 |
< |
Vector3d eField2(0.0); |
| 710 |
< |
potVec longRangePotential(0.0); |
| 721 |
> |
potVec longRangePotential(0.0); |
| 722 |
> |
RealType reciprocalPotential(0.0); |
| 723 |
|
potVec workPot(0.0); |
| 724 |
|
potVec exPot(0.0); |
| 725 |
+ |
Vector3d eField1(0.0); |
| 726 |
+ |
Vector3d eField2(0.0); |
| 727 |
|
vector<int>::iterator ia, jb; |
| 728 |
|
|
| 729 |
|
int loopStart, loopEnd; |
| 730 |
< |
|
| 730 |
> |
|
| 731 |
> |
idat.rcut = &rCut; |
| 732 |
|
idat.vdwMult = &vdwMult; |
| 733 |
|
idat.electroMult = &electroMult; |
| 734 |
|
idat.pot = &workPot; |
| 761 |
|
if (update_nlist) { |
| 762 |
|
if (!usePeriodicBoundaryConditions_) |
| 763 |
|
Mat3x3d bbox = thermo->getBoundingBox(); |
| 764 |
< |
neighborList = fDecomp_->buildNeighborList(); |
| 764 |
> |
fDecomp_->buildNeighborList(neighborList_); |
| 765 |
|
} |
| 766 |
|
} |
| 767 |
|
|
| 768 |
< |
for (vector<pair<int, int> >::iterator it = neighborList.begin(); |
| 769 |
< |
it != neighborList.end(); ++it) { |
| 768 |
> |
for (vector<pair<int, int> >::iterator it = neighborList_.begin(); |
| 769 |
> |
it != neighborList_.end(); ++it) { |
| 770 |
|
|
| 771 |
|
cg1 = (*it).first; |
| 772 |
|
cg2 = (*it).second; |
| 773 |
|
|
| 774 |
< |
cuts = fDecomp_->getGroupCutoffs(cg1, cg2); |
| 774 |
> |
fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq); |
| 775 |
|
|
| 776 |
|
d_grp = fDecomp_->getIntergroupVector(cg1, cg2); |
| 777 |
|
|
| 778 |
< |
curSnapshot->wrapVector(d_grp); |
| 778 |
> |
// already wrapped in the getIntergroupVector call: |
| 779 |
> |
// curSnapshot->wrapVector(d_grp); |
| 780 |
|
rgrpsq = d_grp.lengthSquare(); |
| 765 |
– |
rCutSq = cuts.second; |
| 781 |
|
|
| 782 |
|
if (rgrpsq < rCutSq) { |
| 768 |
– |
idat.rcut = &cuts.first; |
| 783 |
|
if (iLoop == PAIR_LOOP) { |
| 784 |
|
vij = 0.0; |
| 785 |
|
fij.zero(); |
| 836 |
|
|
| 837 |
|
r = sqrt( *(idat.r2) ); |
| 838 |
|
idat.rij = &r; |
| 839 |
< |
|
| 839 |
> |
|
| 840 |
|
if (iLoop == PREPAIR_LOOP) { |
| 841 |
|
interactionMan_->doPrePair(idat); |
| 842 |
|
} else { |
| 932 |
|
} |
| 933 |
|
} |
| 934 |
|
} |
| 935 |
< |
|
| 935 |
> |
|
| 936 |
|
// collects pairwise information |
| 937 |
|
fDecomp_->collectData(); |
| 938 |
+ |
if (cutoffMethod_ == EWALD_FULL) { |
| 939 |
+ |
interactionMan_->doReciprocalSpaceSum(reciprocalPotential); |
| 940 |
+ |
|
| 941 |
+ |
curSnapshot->setReciprocalPotential(reciprocalPotential); |
| 942 |
+ |
} |
| 943 |
|
|
| 944 |
|
if (info_->requiresSelfCorrection()) { |
| 945 |
|
for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
| 961 |
|
|
| 962 |
|
} |
| 963 |
|
|
| 945 |
– |
|
| 964 |
|
void ForceManager::postCalculation() { |
| 965 |
|
|
| 966 |
|
vector<Perturbation*>::iterator pi; |
| 986 |
|
} |
| 987 |
|
|
| 988 |
|
#ifdef IS_MPI |
| 989 |
< |
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9, |
| 990 |
< |
MPI::REALTYPE, MPI::SUM); |
| 989 |
> |
MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9, |
| 990 |
> |
MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
| 991 |
> |
// MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9, |
| 992 |
> |
// MPI::REALTYPE, MPI::SUM); |
| 993 |
|
#endif |
| 994 |
|
curSnapshot->setStressTensor(stressTensor); |
| 995 |
|
|
