47 |
|
* @version 1.0 |
48 |
|
*/ |
49 |
|
|
50 |
+ |
|
51 |
|
#include "brains/ForceManager.hpp" |
52 |
|
#include "primitives/Molecule.hpp" |
53 |
|
#define __OPENMD_C |
59 |
|
#include "nonbonded/NonBondedInteraction.hpp" |
60 |
|
#include "parallel/ForceMatrixDecomposition.hpp" |
61 |
|
|
62 |
+ |
#include <cstdio> |
63 |
+ |
#include <iostream> |
64 |
+ |
#include <iomanip> |
65 |
+ |
|
66 |
|
using namespace std; |
67 |
|
namespace OpenMD { |
68 |
|
|
69 |
|
ForceManager::ForceManager(SimInfo * info) : info_(info) { |
70 |
< |
|
71 |
< |
#ifdef IS_MPI |
72 |
< |
fDecomp_ = new ForceMatrixDecomposition(info_); |
68 |
< |
#else |
69 |
< |
// fDecomp_ = new ForceSerialDecomposition(info); |
70 |
< |
#endif |
70 |
> |
forceField_ = info_->getForceField(); |
71 |
> |
interactionMan_ = new InteractionManager(); |
72 |
> |
fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_); |
73 |
|
} |
74 |
+ |
|
75 |
+ |
/** |
76 |
+ |
* setupCutoffs |
77 |
+ |
* |
78 |
+ |
* Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy |
79 |
+ |
* |
80 |
+ |
* cutoffRadius : realType |
81 |
+ |
* If the cutoffRadius was explicitly set, use that value. |
82 |
+ |
* If the cutoffRadius was not explicitly set: |
83 |
+ |
* Are there electrostatic atoms? Use 12.0 Angstroms. |
84 |
+ |
* No electrostatic atoms? Poll the atom types present in the |
85 |
+ |
* simulation for suggested cutoff values (e.g. 2.5 * sigma). |
86 |
+ |
* Use the maximum suggested value that was found. |
87 |
+ |
* |
88 |
+ |
* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) |
89 |
+ |
* If cutoffMethod was explicitly set, use that choice. |
90 |
+ |
* If cutoffMethod was not explicitly set, use SHIFTED_FORCE |
91 |
+ |
* |
92 |
+ |
* cutoffPolicy : (one of MIX, MAX, TRADITIONAL) |
93 |
+ |
* If cutoffPolicy was explicitly set, use that choice. |
94 |
+ |
* If cutoffPolicy was not explicitly set, use TRADITIONAL |
95 |
+ |
*/ |
96 |
+ |
void ForceManager::setupCutoffs() { |
97 |
+ |
|
98 |
+ |
Globals* simParams_ = info_->getSimParams(); |
99 |
+ |
ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
100 |
+ |
|
101 |
+ |
if (simParams_->haveCutoffRadius()) { |
102 |
+ |
rCut_ = simParams_->getCutoffRadius(); |
103 |
+ |
} else { |
104 |
+ |
if (info_->usesElectrostaticAtoms()) { |
105 |
+ |
sprintf(painCave.errMsg, |
106 |
+ |
"ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n" |
107 |
+ |
"\tOpenMD will use a default value of 12.0 angstroms" |
108 |
+ |
"\tfor the cutoffRadius.\n"); |
109 |
+ |
painCave.isFatal = 0; |
110 |
+ |
painCave.severity = OPENMD_INFO; |
111 |
+ |
simError(); |
112 |
+ |
rCut_ = 12.0; |
113 |
+ |
} else { |
114 |
+ |
RealType thisCut; |
115 |
+ |
set<AtomType*>::iterator i; |
116 |
+ |
set<AtomType*> atomTypes; |
117 |
+ |
atomTypes = info_->getSimulatedAtomTypes(); |
118 |
+ |
for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
119 |
+ |
thisCut = interactionMan_->getSuggestedCutoffRadius((*i)); |
120 |
+ |
rCut_ = max(thisCut, rCut_); |
121 |
+ |
} |
122 |
+ |
sprintf(painCave.errMsg, |
123 |
+ |
"ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n" |
124 |
+ |
"\tOpenMD will use %lf angstroms.\n", |
125 |
+ |
rCut_); |
126 |
+ |
painCave.isFatal = 0; |
127 |
+ |
painCave.severity = OPENMD_INFO; |
128 |
+ |
simError(); |
129 |
+ |
} |
130 |
+ |
} |
131 |
+ |
|
132 |
+ |
fDecomp_->setUserCutoff(rCut_); |
133 |
+ |
|
134 |
+ |
map<string, CutoffMethod> stringToCutoffMethod; |
135 |
+ |
stringToCutoffMethod["HARD"] = HARD; |
136 |
+ |
stringToCutoffMethod["SWITCHED"] = SWITCHED; |
137 |
+ |
stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; |
138 |
+ |
stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; |
139 |
|
|
140 |
< |
void ForceManager::calcForces() { |
140 |
> |
if (simParams_->haveCutoffMethod()) { |
141 |
> |
string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); |
142 |
> |
map<string, CutoffMethod>::iterator i; |
143 |
> |
i = stringToCutoffMethod.find(cutMeth); |
144 |
> |
if (i == stringToCutoffMethod.end()) { |
145 |
> |
sprintf(painCave.errMsg, |
146 |
> |
"ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" |
147 |
> |
"\tShould be one of: " |
148 |
> |
"HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", |
149 |
> |
cutMeth.c_str()); |
150 |
> |
painCave.isFatal = 1; |
151 |
> |
painCave.severity = OPENMD_ERROR; |
152 |
> |
simError(); |
153 |
> |
} else { |
154 |
> |
cutoffMethod_ = i->second; |
155 |
> |
} |
156 |
> |
} else { |
157 |
> |
sprintf(painCave.errMsg, |
158 |
> |
"ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n" |
159 |
> |
"\tOpenMD will use SHIFTED_FORCE.\n"); |
160 |
> |
painCave.isFatal = 0; |
161 |
> |
painCave.severity = OPENMD_INFO; |
162 |
> |
simError(); |
163 |
> |
cutoffMethod_ = SHIFTED_FORCE; |
164 |
> |
} |
165 |
> |
|
166 |
> |
map<string, CutoffPolicy> stringToCutoffPolicy; |
167 |
> |
stringToCutoffPolicy["MIX"] = MIX; |
168 |
> |
stringToCutoffPolicy["MAX"] = MAX; |
169 |
> |
stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL; |
170 |
> |
|
171 |
> |
std::string cutPolicy; |
172 |
> |
if (forceFieldOptions_.haveCutoffPolicy()){ |
173 |
> |
cutPolicy = forceFieldOptions_.getCutoffPolicy(); |
174 |
> |
}else if (simParams_->haveCutoffPolicy()) { |
175 |
> |
cutPolicy = simParams_->getCutoffPolicy(); |
176 |
> |
} |
177 |
> |
|
178 |
> |
if (!cutPolicy.empty()){ |
179 |
> |
toUpper(cutPolicy); |
180 |
> |
map<string, CutoffPolicy>::iterator i; |
181 |
> |
i = stringToCutoffPolicy.find(cutPolicy); |
182 |
> |
|
183 |
> |
if (i == stringToCutoffPolicy.end()) { |
184 |
> |
sprintf(painCave.errMsg, |
185 |
> |
"ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n" |
186 |
> |
"\tShould be one of: " |
187 |
> |
"MIX, MAX, or TRADITIONAL\n", |
188 |
> |
cutPolicy.c_str()); |
189 |
> |
painCave.isFatal = 1; |
190 |
> |
painCave.severity = OPENMD_ERROR; |
191 |
> |
simError(); |
192 |
> |
} else { |
193 |
> |
cutoffPolicy_ = i->second; |
194 |
> |
} |
195 |
> |
} else { |
196 |
> |
sprintf(painCave.errMsg, |
197 |
> |
"ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n" |
198 |
> |
"\tOpenMD will use TRADITIONAL.\n"); |
199 |
> |
painCave.isFatal = 0; |
200 |
> |
painCave.severity = OPENMD_INFO; |
201 |
> |
simError(); |
202 |
> |
cutoffPolicy_ = TRADITIONAL; |
203 |
> |
} |
204 |
> |
fDecomp_->setCutoffPolicy(cutoffPolicy_); |
205 |
> |
} |
206 |
> |
|
207 |
> |
/** |
208 |
> |
* setupSwitching |
209 |
> |
* |
210 |
> |
* Sets the values of switchingRadius and |
211 |
> |
* If the switchingRadius was explicitly set, use that value (but check it) |
212 |
> |
* If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ |
213 |
> |
*/ |
214 |
> |
void ForceManager::setupSwitching() { |
215 |
> |
Globals* simParams_ = info_->getSimParams(); |
216 |
> |
|
217 |
> |
// create the switching function object: |
218 |
> |
switcher_ = new SwitchingFunction(); |
219 |
|
|
220 |
+ |
if (simParams_->haveSwitchingRadius()) { |
221 |
+ |
rSwitch_ = simParams_->getSwitchingRadius(); |
222 |
+ |
if (rSwitch_ > rCut_) { |
223 |
+ |
sprintf(painCave.errMsg, |
224 |
+ |
"ForceManager::setupSwitching: switchingRadius (%f) is larger " |
225 |
+ |
"than the cutoffRadius(%f)\n", rSwitch_, rCut_); |
226 |
+ |
painCave.isFatal = 1; |
227 |
+ |
painCave.severity = OPENMD_ERROR; |
228 |
+ |
simError(); |
229 |
+ |
} |
230 |
+ |
} else { |
231 |
+ |
rSwitch_ = 0.85 * rCut_; |
232 |
+ |
sprintf(painCave.errMsg, |
233 |
+ |
"ForceManager::setupSwitching: No value was set for the switchingRadius.\n" |
234 |
+ |
"\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
235 |
+ |
"\tswitchingRadius = %f. for this simulation\n", rSwitch_); |
236 |
+ |
painCave.isFatal = 0; |
237 |
+ |
painCave.severity = OPENMD_WARNING; |
238 |
+ |
simError(); |
239 |
+ |
} |
240 |
+ |
|
241 |
+ |
// Default to cubic switching function. |
242 |
+ |
sft_ = cubic; |
243 |
+ |
if (simParams_->haveSwitchingFunctionType()) { |
244 |
+ |
string funcType = simParams_->getSwitchingFunctionType(); |
245 |
+ |
toUpper(funcType); |
246 |
+ |
if (funcType == "CUBIC") { |
247 |
+ |
sft_ = cubic; |
248 |
+ |
} else { |
249 |
+ |
if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
250 |
+ |
sft_ = fifth_order_poly; |
251 |
+ |
} else { |
252 |
+ |
// throw error |
253 |
+ |
sprintf( painCave.errMsg, |
254 |
+ |
"ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n" |
255 |
+ |
"\tswitchingFunctionType must be one of: " |
256 |
+ |
"\"cubic\" or \"fifth_order_polynomial\".", |
257 |
+ |
funcType.c_str() ); |
258 |
+ |
painCave.isFatal = 1; |
259 |
+ |
painCave.severity = OPENMD_ERROR; |
260 |
+ |
simError(); |
261 |
+ |
} |
262 |
+ |
} |
263 |
+ |
} |
264 |
+ |
switcher_->setSwitchType(sft_); |
265 |
+ |
switcher_->setSwitch(rSwitch_, rCut_); |
266 |
+ |
} |
267 |
+ |
|
268 |
+ |
void ForceManager::initialize() { |
269 |
+ |
|
270 |
|
if (!info_->isTopologyDone()) { |
271 |
|
info_->update(); |
272 |
|
interactionMan_->setSimInfo(info_); |
273 |
|
interactionMan_->initialize(); |
274 |
< |
swfun_ = interactionMan_->getSwitchingFunction(); |
275 |
< |
fDecomp_->distributeInitialData(); |
276 |
< |
info_->prepareTopology(); |
274 |
> |
|
275 |
> |
// We want to delay the cutoffs until after the interaction |
276 |
> |
// manager has set up the atom-atom interactions so that we can |
277 |
> |
// query them for suggested cutoff values |
278 |
> |
|
279 |
> |
setupCutoffs(); |
280 |
> |
setupSwitching(); |
281 |
> |
|
282 |
> |
info_->prepareTopology(); |
283 |
|
} |
284 |
+ |
|
285 |
+ |
ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); |
286 |
|
|
287 |
+ |
// Force fields can set options on how to scale van der Waals and electrostatic |
288 |
+ |
// interactions for atoms connected via bonds, bends and torsions |
289 |
+ |
// in this case the topological distance between atoms is: |
290 |
+ |
// 0 = topologically unconnected |
291 |
+ |
// 1 = bonded together |
292 |
+ |
// 2 = connected via a bend |
293 |
+ |
// 3 = connected via a torsion |
294 |
+ |
|
295 |
+ |
vdwScale_.reserve(4); |
296 |
+ |
fill(vdwScale_.begin(), vdwScale_.end(), 0.0); |
297 |
+ |
|
298 |
+ |
electrostaticScale_.reserve(4); |
299 |
+ |
fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0); |
300 |
+ |
|
301 |
+ |
vdwScale_[0] = 1.0; |
302 |
+ |
vdwScale_[1] = fopts.getvdw12scale(); |
303 |
+ |
vdwScale_[2] = fopts.getvdw13scale(); |
304 |
+ |
vdwScale_[3] = fopts.getvdw14scale(); |
305 |
+ |
|
306 |
+ |
electrostaticScale_[0] = 1.0; |
307 |
+ |
electrostaticScale_[1] = fopts.getelectrostatic12scale(); |
308 |
+ |
electrostaticScale_[2] = fopts.getelectrostatic13scale(); |
309 |
+ |
electrostaticScale_[3] = fopts.getelectrostatic14scale(); |
310 |
+ |
|
311 |
+ |
fDecomp_->distributeInitialData(); |
312 |
+ |
|
313 |
+ |
initialized_ = true; |
314 |
+ |
|
315 |
+ |
} |
316 |
+ |
|
317 |
+ |
void ForceManager::calcForces() { |
318 |
+ |
|
319 |
+ |
if (!initialized_) initialize(); |
320 |
+ |
|
321 |
|
preCalculation(); |
322 |
|
shortRangeInteractions(); |
323 |
|
longRangeInteractions(); |
324 |
< |
postCalculation(); |
88 |
< |
|
324 |
> |
postCalculation(); |
325 |
|
} |
326 |
|
|
327 |
|
void ForceManager::preCalculation() { |
481 |
|
|
482 |
|
void ForceManager::longRangeInteractions() { |
483 |
|
|
248 |
– |
// some of this initial stuff will go away: |
484 |
|
Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
485 |
|
DataStorage* config = &(curSnapshot->atomData); |
486 |
|
DataStorage* cgConfig = &(curSnapshot->cgData); |
252 |
– |
RealType* frc = config->getArrayPointer(DataStorage::dslForce); |
253 |
– |
RealType* pos = config->getArrayPointer(DataStorage::dslPosition); |
254 |
– |
RealType* trq = config->getArrayPointer(DataStorage::dslTorque); |
255 |
– |
RealType* A = config->getArrayPointer(DataStorage::dslAmat); |
256 |
– |
RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame); |
257 |
– |
RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot); |
258 |
– |
RealType* rc; |
487 |
|
|
488 |
< |
if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){ |
489 |
< |
rc = cgConfig->getArrayPointer(DataStorage::dslPosition); |
488 |
> |
//calculate the center of mass of cutoff group |
489 |
> |
|
490 |
> |
SimInfo::MoleculeIterator mi; |
491 |
> |
Molecule* mol; |
492 |
> |
Molecule::CutoffGroupIterator ci; |
493 |
> |
CutoffGroup* cg; |
494 |
> |
|
495 |
> |
if(info_->getNCutoffGroups() > 0){ |
496 |
> |
for (mol = info_->beginMolecule(mi); mol != NULL; |
497 |
> |
mol = info_->nextMolecule(mi)) { |
498 |
> |
for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
499 |
> |
cg = mol->nextCutoffGroup(ci)) { |
500 |
> |
cg->updateCOM(); |
501 |
> |
} |
502 |
> |
} |
503 |
|
} else { |
504 |
|
// center of mass of the group is the same as position of the atom |
505 |
|
// if cutoff group does not exist |
506 |
< |
rc = pos; |
266 |
< |
} |
267 |
< |
|
268 |
< |
//initialize data before passing to fortran |
269 |
< |
RealType longRangePotential[N_INTERACTION_FAMILIES]; |
270 |
< |
RealType lrPot = 0.0; |
271 |
< |
int isError = 0; |
272 |
< |
|
273 |
< |
// dangerous to iterate over enums, but we'll live on the edge: |
274 |
< |
for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){ |
275 |
< |
longRangePotential[i]=0.0; //Initialize array |
506 |
> |
cgConfig->position = config->position; |
507 |
|
} |
508 |
|
|
509 |
< |
// new stuff starts here: |
279 |
< |
|
509 |
> |
fDecomp_->zeroWorkArrays(); |
510 |
|
fDecomp_->distributeData(); |
511 |
< |
|
512 |
< |
int cg1, cg2, atom1, atom2; |
513 |
< |
Vector3d d_grp, dag; |
514 |
< |
RealType rgrpsq, rgrp; |
511 |
> |
|
512 |
> |
int cg1, cg2, atom1, atom2, topoDist; |
513 |
> |
Vector3d d_grp, dag, d; |
514 |
> |
RealType rgrpsq, rgrp, r2, r; |
515 |
> |
RealType electroMult, vdwMult; |
516 |
|
RealType vij; |
517 |
< |
Vector3d fij, fg; |
518 |
< |
pair<int, int> gtypes; |
517 |
> |
Vector3d fij, fg, f1; |
518 |
> |
tuple3<RealType, RealType, RealType> cuts; |
519 |
|
RealType rCutSq; |
520 |
|
bool in_switching_region; |
521 |
|
RealType sw, dswdr, swderiv; |
523 |
|
InteractionData idat; |
524 |
|
SelfData sdat; |
525 |
|
RealType mf; |
526 |
+ |
RealType lrPot; |
527 |
+ |
RealType vpair; |
528 |
+ |
potVec longRangePotential(0.0); |
529 |
+ |
potVec workPot(0.0); |
530 |
|
|
531 |
|
int loopStart, loopEnd; |
532 |
|
|
533 |
+ |
idat.vdwMult = &vdwMult; |
534 |
+ |
idat.electroMult = &electroMult; |
535 |
+ |
idat.pot = &workPot; |
536 |
+ |
sdat.pot = fDecomp_->getEmbeddingPotential(); |
537 |
+ |
idat.vpair = &vpair; |
538 |
+ |
idat.f1 = &f1; |
539 |
+ |
idat.sw = &sw; |
540 |
+ |
idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false; |
541 |
+ |
idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false; |
542 |
+ |
|
543 |
|
loopEnd = PAIR_LOOP; |
544 |
|
if (info_->requiresPrepair() ) { |
545 |
|
loopStart = PREPAIR_LOOP; |
546 |
|
} else { |
547 |
|
loopStart = PAIR_LOOP; |
548 |
|
} |
549 |
< |
|
550 |
< |
for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) { |
551 |
< |
|
549 |
> |
|
550 |
> |
for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) { |
551 |
> |
|
552 |
|
if (iLoop == loopStart) { |
553 |
|
bool update_nlist = fDecomp_->checkNeighborList(); |
554 |
|
if (update_nlist) |
555 |
|
neighborList = fDecomp_->buildNeighborList(); |
556 |
< |
} |
557 |
< |
|
556 |
> |
} |
557 |
> |
|
558 |
|
for (vector<pair<int, int> >::iterator it = neighborList.begin(); |
559 |
|
it != neighborList.end(); ++it) { |
560 |
< |
|
560 |
> |
|
561 |
|
cg1 = (*it).first; |
562 |
|
cg2 = (*it).second; |
563 |
+ |
|
564 |
+ |
cuts = fDecomp_->getGroupCutoffs(cg1, cg2); |
565 |
|
|
319 |
– |
gtypes = fDecomp_->getGroupTypes(cg1, cg2); |
566 |
|
d_grp = fDecomp_->getIntergroupVector(cg1, cg2); |
567 |
|
curSnapshot->wrapVector(d_grp); |
568 |
|
rgrpsq = d_grp.lengthSquare(); |
323 |
– |
rCutSq = groupCutoffMap[gtypes].first; |
569 |
|
|
570 |
+ |
rCutSq = cuts.second; |
571 |
+ |
|
572 |
|
if (rgrpsq < rCutSq) { |
573 |
< |
*(idat.rcut) = groupCutoffMap[gtypes].second; |
573 |
> |
idat.rcut = &cuts.first; |
574 |
|
if (iLoop == PAIR_LOOP) { |
575 |
|
vij *= 0.0; |
576 |
|
fij = V3Zero; |
577 |
|
} |
578 |
|
|
579 |
< |
in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr, |
580 |
< |
rgrp); |
579 |
> |
in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, |
580 |
> |
rgrp); |
581 |
> |
|
582 |
|
atomListRow = fDecomp_->getAtomsInGroupRow(cg1); |
583 |
|
atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2); |
584 |
|
|
589 |
|
for (vector<int>::iterator jb = atomListColumn.begin(); |
590 |
|
jb != atomListColumn.end(); ++jb) { |
591 |
|
atom2 = (*jb); |
592 |
< |
|
592 |
> |
|
593 |
|
if (!fDecomp_->skipAtomPair(atom1, atom2)) { |
594 |
+ |
vpair = 0.0; |
595 |
+ |
workPot = 0.0; |
596 |
+ |
f1 = V3Zero; |
597 |
+ |
|
598 |
+ |
fDecomp_->fillInteractionData(idat, atom1, atom2); |
599 |
|
|
600 |
< |
idat = fDecomp_->fillInteractionData(atom1, atom2); |
600 |
> |
topoDist = fDecomp_->getTopologicalDistance(atom1, atom2); |
601 |
> |
vdwMult = vdwScale_[topoDist]; |
602 |
> |
electroMult = electrostaticScale_[topoDist]; |
603 |
|
|
604 |
|
if (atomListRow.size() == 1 && atomListColumn.size() == 1) { |
605 |
< |
*(idat.d) = d_grp; |
606 |
< |
*(idat.r2) = rgrpsq; |
605 |
> |
idat.d = &d_grp; |
606 |
> |
idat.r2 = &rgrpsq; |
607 |
|
} else { |
608 |
< |
*(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2); |
609 |
< |
curSnapshot->wrapVector( *(idat.d) ); |
610 |
< |
*(idat.r2) = idat.d->lengthSquare(); |
608 |
> |
d = fDecomp_->getInteratomicVector(atom1, atom2); |
609 |
> |
curSnapshot->wrapVector( d ); |
610 |
> |
r2 = d.lengthSquare(); |
611 |
> |
idat.d = &d; |
612 |
> |
idat.r2 = &r2; |
613 |
|
} |
614 |
|
|
615 |
< |
*(idat.rij) = sqrt( *(idat.r2) ); |
615 |
> |
r = sqrt( *(idat.r2) ); |
616 |
> |
idat.rij = &r; |
617 |
|
|
618 |
|
if (iLoop == PREPAIR_LOOP) { |
619 |
|
interactionMan_->doPrePair(idat); |
620 |
|
} else { |
621 |
|
interactionMan_->doPair(idat); |
622 |
< |
vij += *(idat.vpair); |
623 |
< |
fij += *(idat.f1); |
624 |
< |
tau -= outProduct( *(idat.d), *(idat.f1)); |
622 |
> |
fDecomp_->unpackInteractionData(idat, atom1, atom2); |
623 |
> |
vij += vpair; |
624 |
> |
fij += f1; |
625 |
> |
tau -= outProduct( *(idat.d), f1); |
626 |
|
} |
627 |
|
} |
628 |
|
} |
686 |
|
if (iLoop == PREPAIR_LOOP) { |
687 |
|
if (info_->requiresPrepair()) { |
688 |
|
fDecomp_->collectIntermediateData(); |
689 |
< |
atomListLocal = fDecomp_->getAtomList(); |
690 |
< |
for (vector<int>::iterator ia = atomListLocal.begin(); |
691 |
< |
ia != atomListLocal.end(); ++ia) { |
433 |
< |
atom1 = (*ia); |
434 |
< |
sdat = fDecomp_->fillSelfData(atom1); |
689 |
> |
|
690 |
> |
for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
691 |
> |
fDecomp_->fillSelfData(sdat, atom1); |
692 |
|
interactionMan_->doPreForce(sdat); |
693 |
|
} |
694 |
+ |
|
695 |
+ |
|
696 |
|
fDecomp_->distributeIntermediateData(); |
697 |
|
} |
698 |
|
} |
701 |
|
|
702 |
|
fDecomp_->collectData(); |
703 |
|
|
704 |
< |
if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) { |
705 |
< |
atomListLocal = fDecomp_->getAtomList(); |
706 |
< |
for (vector<int>::iterator ia = atomListLocal.begin(); |
448 |
< |
ia != atomListLocal.end(); ++ia) { |
449 |
< |
atom1 = (*ia); |
704 |
> |
if ( info_->requiresSkipCorrection() ) { |
705 |
> |
|
706 |
> |
for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) { |
707 |
|
|
708 |
< |
if (info_->requiresSkipCorrection()) { |
709 |
< |
vector<int> skipList = fDecomp_->getSkipsForAtom(atom1); |
710 |
< |
for (vector<int>::iterator jb = skipList.begin(); |
711 |
< |
jb != skipList.end(); ++jb) { |
712 |
< |
atom2 = (*jb); |
713 |
< |
idat = fDecomp_->fillSkipData(atom1, atom2); |
714 |
< |
interactionMan_->doSkipCorrection(idat); |
715 |
< |
} |
708 |
> |
vector<int> skipList = fDecomp_->getSkipsForAtom( atom1 ); |
709 |
> |
|
710 |
> |
for (vector<int>::iterator jb = skipList.begin(); |
711 |
> |
jb != skipList.end(); ++jb) { |
712 |
> |
|
713 |
> |
atom2 = (*jb); |
714 |
> |
fDecomp_->fillSkipData(idat, atom1, atom2); |
715 |
> |
interactionMan_->doSkipCorrection(idat); |
716 |
> |
fDecomp_->unpackSkipData(idat, atom1, atom2); |
717 |
> |
|
718 |
|
} |
460 |
– |
|
461 |
– |
if (info_->requiresSelfCorrection()) { |
462 |
– |
sdat = fDecomp_->fillSelfData(atom1); |
463 |
– |
interactionMan_->doSelfCorrection(sdat); |
464 |
– |
} |
719 |
|
} |
720 |
|
} |
721 |
+ |
|
722 |
+ |
if (info_->requiresSelfCorrection()) { |
723 |
|
|
724 |
< |
// dangerous to iterate over enums, but we'll live on the edge: |
725 |
< |
for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){ |
726 |
< |
lrPot += longRangePotential[i]; //Quick hack |
724 |
> |
for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
725 |
> |
fDecomp_->fillSelfData(sdat, atom1); |
726 |
> |
interactionMan_->doSelfCorrection(sdat); |
727 |
> |
} |
728 |
> |
|
729 |
|
} |
730 |
< |
|
730 |
> |
|
731 |
> |
longRangePotential = *(fDecomp_->getEmbeddingPotential()) + |
732 |
> |
*(fDecomp_->getPairwisePotential()); |
733 |
> |
|
734 |
> |
lrPot = longRangePotential.sum(); |
735 |
> |
|
736 |
|
//store the tau and long range potential |
737 |
|
curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; |
738 |
|
curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY]; |