55 |
|
#include "primitives/Bend.hpp" |
56 |
|
#include "primitives/Torsion.hpp" |
57 |
|
#include "primitives/Inversion.hpp" |
58 |
– |
#include "parallel/ForceMatrixDecomposition.hpp" |
58 |
|
#include "nonbonded/NonBondedInteraction.hpp" |
59 |
+ |
#include "parallel/ForceMatrixDecomposition.hpp" |
60 |
|
|
61 |
|
using namespace std; |
62 |
|
namespace OpenMD { |
72 |
|
|
73 |
|
void ForceManager::calcForces() { |
74 |
|
|
75 |
< |
if (!info_->isFortranInitialized()) { |
75 |
> |
if (!info_->isTopologyDone()) { |
76 |
|
info_->update(); |
77 |
|
interactionMan_->setSimInfo(info_); |
78 |
|
interactionMan_->initialize(); |
79 |
|
swfun_ = interactionMan_->getSwitchingFunction(); |
80 |
|
fDecomp_->distributeInitialData(); |
81 |
< |
info_->setupFortran(); |
81 |
> |
info_->prepareTopology(); |
82 |
|
} |
83 |
|
|
84 |
|
preCalculation(); |
323 |
|
rCutSq = groupCutoffMap[gtypes].first; |
324 |
|
|
325 |
|
if (rgrpsq < rCutSq) { |
326 |
< |
idat.rcut = groupCutoffMap[gtypes].second; |
326 |
> |
*(idat.rcut) = groupCutoffMap[gtypes].second; |
327 |
|
if (iLoop == PAIR_LOOP) { |
328 |
|
vij *= 0.0; |
329 |
|
fij = V3Zero; |
330 |
|
} |
331 |
|
|
332 |
< |
in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, dswdr, rgrp); |
332 |
> |
in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr, |
333 |
> |
rgrp); |
334 |
|
atomListRow = fDecomp_->getAtomsInGroupRow(cg1); |
335 |
|
atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2); |
336 |
|
|
347 |
|
idat = fDecomp_->fillInteractionData(atom1, atom2); |
348 |
|
|
349 |
|
if (atomListRow.size() == 1 && atomListColumn.size() == 1) { |
350 |
< |
idat.d = d_grp; |
351 |
< |
idat.r2 = rgrpsq; |
350 |
> |
*(idat.d) = d_grp; |
351 |
> |
*(idat.r2) = rgrpsq; |
352 |
|
} else { |
353 |
< |
idat.d = fDecomp_->getInteratomicVector(atom1, atom2); |
354 |
< |
curSnapshot->wrapVector(idat.d); |
355 |
< |
idat.r2 = idat.d.lengthSquare(); |
353 |
> |
*(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2); |
354 |
> |
curSnapshot->wrapVector( *(idat.d) ); |
355 |
> |
*(idat.r2) = idat.d->lengthSquare(); |
356 |
|
} |
357 |
|
|
358 |
< |
idat.rij = sqrt(idat.r2); |
358 |
> |
*(idat.rij) = sqrt( *(idat.r2) ); |
359 |
|
|
360 |
|
if (iLoop == PREPAIR_LOOP) { |
361 |
|
interactionMan_->doPrePair(idat); |
362 |
|
} else { |
363 |
|
interactionMan_->doPair(idat); |
364 |
< |
vij += idat.vpair; |
365 |
< |
fij += idat.f1; |
366 |
< |
tau -= outProduct(idat.d, idat.f1); |
364 |
> |
vij += *(idat.vpair); |
365 |
> |
fij += *(idat.f1); |
366 |
> |
tau -= outProduct( *(idat.d), *(idat.f1)); |
367 |
|
} |
368 |
|
} |
369 |
|
} |
377 |
|
fij += fg; |
378 |
|
|
379 |
|
if (atomListRow.size() == 1 && atomListColumn.size() == 1) { |
380 |
< |
tau -= outProduct(idat.d, fg); |
380 |
> |
tau -= outProduct( *(idat.d), fg); |
381 |
|
} |
382 |
|
|
383 |
|
for (vector<int>::iterator ia = atomListRow.begin(); |
384 |
|
ia != atomListRow.end(); ++ia) { |
385 |
|
atom1 = (*ia); |
386 |
< |
mf = fDecomp_->getMfactRow(atom1); |
386 |
> |
mf = fDecomp_->getMassFactorRow(atom1); |
387 |
|
// fg is the force on atom ia due to cutoff group's |
388 |
|
// presence in switching region |
389 |
|
fg = swderiv * d_grp * mf; |
401 |
|
for (vector<int>::iterator jb = atomListColumn.begin(); |
402 |
|
jb != atomListColumn.end(); ++jb) { |
403 |
|
atom2 = (*jb); |
404 |
< |
mf = fDecomp_->getMfactColumn(atom2); |
404 |
> |
mf = fDecomp_->getMassFactorColumn(atom2); |
405 |
|
// fg is the force on atom jb due to cutoff group's |
406 |
|
// presence in switching region |
407 |
|
fg = -swderiv * d_grp * mf; |
427 |
|
if (iLoop == PREPAIR_LOOP) { |
428 |
|
if (info_->requiresPrepair()) { |
429 |
|
fDecomp_->collectIntermediateData(); |
430 |
< |
atomListLocal = fDecomp_->getAtomList(); |
431 |
< |
for (vector<int>::iterator ia = atomListLocal.begin(); |
431 |
< |
ia != atomListLocal.end(); ++ia) { |
432 |
< |
atom1 = (*ia); |
430 |
> |
|
431 |
> |
for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
432 |
|
sdat = fDecomp_->fillSelfData(atom1); |
433 |
|
interactionMan_->doPreForce(sdat); |
434 |
|
} |
435 |
+ |
|
436 |
|
fDecomp_->distributeIntermediateData(); |
437 |
|
} |
438 |
|
} |
441 |
|
|
442 |
|
fDecomp_->collectData(); |
443 |
|
|
444 |
< |
if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) { |
445 |
< |
atomListLocal = fDecomp_->getAtomList(); |
446 |
< |
for (vector<int>::iterator ia = atomListLocal.begin(); |
447 |
< |
ia != atomListLocal.end(); ++ia) { |
448 |
< |
atom1 = (*ia); |
444 |
> |
if ( info_->requiresSkipCorrection() ) { |
445 |
> |
|
446 |
> |
for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) { |
447 |
|
|
448 |
< |
if (info_->requiresSkipCorrection()) { |
449 |
< |
vector<int> skipList = fDecomp_->getSkipsForAtom(atom1); |
450 |
< |
for (vector<int>::iterator jb = skipList.begin(); |
451 |
< |
jb != skipList.end(); ++jb) { |
452 |
< |
atom2 = (*jb); |
453 |
< |
idat = fDecomp_->fillSkipData(atom1, atom2); |
454 |
< |
interactionMan_->doSkipCorrection(idat); |
455 |
< |
} |
448 |
> |
vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 ); |
449 |
> |
|
450 |
> |
for (vector<int>::iterator jb = skipList.begin(); |
451 |
> |
jb != skipList.end(); ++jb) { |
452 |
> |
|
453 |
> |
atom2 = (*jb); |
454 |
> |
idat = fDecomp_->fillSkipData(atom1, atom2); |
455 |
> |
interactionMan_->doSkipCorrection(idat); |
456 |
> |
|
457 |
|
} |
459 |
– |
|
460 |
– |
if (info_->requiresSelfCorrection()) { |
461 |
– |
sdat = fDecomp_->fillSelfData(atom1); |
462 |
– |
interactionMan_->doSelfCorrection(sdat); |
463 |
– |
} |
458 |
|
} |
459 |
|
} |
460 |
+ |
|
461 |
+ |
if (info_->requiresSelfCorrection()) { |
462 |
|
|
463 |
+ |
for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { |
464 |
+ |
sdat = fDecomp_->fillSelfData(atom1); |
465 |
+ |
interactionMan_->doSelfCorrection(sdat); |
466 |
+ |
} |
467 |
+ |
|
468 |
+ |
} |
469 |
+ |
|
470 |
|
// dangerous to iterate over enums, but we'll live on the edge: |
471 |
|
for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){ |
472 |
|
lrPot += longRangePotential[i]; //Quick hack |