48 |
|
|
49 |
|
#include <algorithm> |
50 |
|
#include <set> |
51 |
+ |
#include <map> |
52 |
|
|
53 |
|
#include "brains/SimInfo.hpp" |
54 |
|
#include "math/Vector3.hpp" |
55 |
|
#include "primitives/Molecule.hpp" |
56 |
|
#include "UseTheForce/fCutoffPolicy.h" |
57 |
|
#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
58 |
< |
#include "UseTheForce/DarkSide/fScreeningMethod.h" |
58 |
> |
#include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" |
59 |
> |
#include "UseTheForce/DarkSide/fSwitchingFunctionType.h" |
60 |
|
#include "UseTheForce/doForces_interface.h" |
61 |
|
#include "UseTheForce/DarkSide/electrostatic_interface.h" |
62 |
|
#include "UseTheForce/notifyCutoffs_interface.h" |
63 |
+ |
#include "UseTheForce/DarkSide/switcheroo_interface.h" |
64 |
|
#include "utils/MemoryUtils.hpp" |
65 |
|
#include "utils/simError.h" |
66 |
|
#include "selection/SelectionManager.hpp" |
71 |
|
#endif |
72 |
|
|
73 |
|
namespace oopse { |
74 |
+ |
std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { |
75 |
+ |
std::map<int, std::set<int> >::iterator i = container.find(index); |
76 |
+ |
std::set<int> result; |
77 |
+ |
if (i != container.end()) { |
78 |
+ |
result = i->second; |
79 |
+ |
} |
80 |
|
|
81 |
+ |
return result; |
82 |
+ |
} |
83 |
+ |
|
84 |
|
SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, |
85 |
|
ForceField* ff, Globals* simParams) : |
86 |
|
stamps_(stamps), forceField_(ff), simParams_(simParams), |
355 |
|
int b; |
356 |
|
int c; |
357 |
|
int d; |
358 |
+ |
|
359 |
+ |
std::map<int, std::set<int> > atomGroups; |
360 |
+ |
|
361 |
+ |
Molecule::RigidBodyIterator rbIter; |
362 |
+ |
RigidBody* rb; |
363 |
+ |
Molecule::IntegrableObjectIterator ii; |
364 |
+ |
StuntDouble* integrableObject; |
365 |
+ |
|
366 |
+ |
for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
367 |
+ |
integrableObject = mol->nextIntegrableObject(ii)) { |
368 |
+ |
|
369 |
+ |
if (integrableObject->isRigidBody()) { |
370 |
+ |
rb = static_cast<RigidBody*>(integrableObject); |
371 |
+ |
std::vector<Atom*> atoms = rb->getAtoms(); |
372 |
+ |
std::set<int> rigidAtoms; |
373 |
+ |
for (int i = 0; i < atoms.size(); ++i) { |
374 |
+ |
rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
375 |
+ |
} |
376 |
+ |
for (int i = 0; i < atoms.size(); ++i) { |
377 |
+ |
atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
378 |
+ |
} |
379 |
+ |
} else { |
380 |
+ |
std::set<int> oneAtomSet; |
381 |
+ |
oneAtomSet.insert(integrableObject->getGlobalIndex()); |
382 |
+ |
atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
383 |
+ |
} |
384 |
+ |
} |
385 |
+ |
|
386 |
|
|
387 |
+ |
|
388 |
|
for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
389 |
|
a = bond->getAtomA()->getGlobalIndex(); |
390 |
|
b = bond->getAtomB()->getGlobalIndex(); |
395 |
|
a = bend->getAtomA()->getGlobalIndex(); |
396 |
|
b = bend->getAtomB()->getGlobalIndex(); |
397 |
|
c = bend->getAtomC()->getGlobalIndex(); |
398 |
+ |
std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
399 |
+ |
std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
400 |
+ |
std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
401 |
|
|
402 |
< |
exclude_.addPair(a, b); |
403 |
< |
exclude_.addPair(a, c); |
404 |
< |
exclude_.addPair(b, c); |
402 |
> |
exclude_.addPairs(rigidSetA, rigidSetB); |
403 |
> |
exclude_.addPairs(rigidSetA, rigidSetC); |
404 |
> |
exclude_.addPairs(rigidSetB, rigidSetC); |
405 |
> |
|
406 |
> |
//exclude_.addPair(a, b); |
407 |
> |
//exclude_.addPair(a, c); |
408 |
> |
//exclude_.addPair(b, c); |
409 |
|
} |
410 |
|
|
411 |
|
for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
413 |
|
b = torsion->getAtomB()->getGlobalIndex(); |
414 |
|
c = torsion->getAtomC()->getGlobalIndex(); |
415 |
|
d = torsion->getAtomD()->getGlobalIndex(); |
416 |
+ |
std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
417 |
+ |
std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
418 |
+ |
std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
419 |
+ |
std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
420 |
|
|
421 |
+ |
exclude_.addPairs(rigidSetA, rigidSetB); |
422 |
+ |
exclude_.addPairs(rigidSetA, rigidSetC); |
423 |
+ |
exclude_.addPairs(rigidSetA, rigidSetD); |
424 |
+ |
exclude_.addPairs(rigidSetB, rigidSetC); |
425 |
+ |
exclude_.addPairs(rigidSetB, rigidSetD); |
426 |
+ |
exclude_.addPairs(rigidSetC, rigidSetD); |
427 |
+ |
|
428 |
+ |
/* |
429 |
+ |
exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
430 |
+ |
exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
431 |
+ |
exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
432 |
+ |
exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
433 |
+ |
exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
434 |
+ |
exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
435 |
+ |
|
436 |
+ |
|
437 |
|
exclude_.addPair(a, b); |
438 |
|
exclude_.addPair(a, c); |
439 |
|
exclude_.addPair(a, d); |
440 |
|
exclude_.addPair(b, c); |
441 |
|
exclude_.addPair(b, d); |
442 |
|
exclude_.addPair(c, d); |
443 |
+ |
*/ |
444 |
|
} |
445 |
|
|
377 |
– |
Molecule::RigidBodyIterator rbIter; |
378 |
– |
RigidBody* rb; |
446 |
|
for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
447 |
|
std::vector<Atom*> atoms = rb->getAtoms(); |
448 |
|
for (int i = 0; i < atoms.size() -1 ; ++i) { |
467 |
|
int b; |
468 |
|
int c; |
469 |
|
int d; |
470 |
+ |
|
471 |
+ |
std::map<int, std::set<int> > atomGroups; |
472 |
+ |
|
473 |
+ |
Molecule::RigidBodyIterator rbIter; |
474 |
+ |
RigidBody* rb; |
475 |
+ |
Molecule::IntegrableObjectIterator ii; |
476 |
+ |
StuntDouble* integrableObject; |
477 |
|
|
478 |
+ |
for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
479 |
+ |
integrableObject = mol->nextIntegrableObject(ii)) { |
480 |
+ |
|
481 |
+ |
if (integrableObject->isRigidBody()) { |
482 |
+ |
rb = static_cast<RigidBody*>(integrableObject); |
483 |
+ |
std::vector<Atom*> atoms = rb->getAtoms(); |
484 |
+ |
std::set<int> rigidAtoms; |
485 |
+ |
for (int i = 0; i < atoms.size(); ++i) { |
486 |
+ |
rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
487 |
+ |
} |
488 |
+ |
for (int i = 0; i < atoms.size(); ++i) { |
489 |
+ |
atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
490 |
+ |
} |
491 |
+ |
} else { |
492 |
+ |
std::set<int> oneAtomSet; |
493 |
+ |
oneAtomSet.insert(integrableObject->getGlobalIndex()); |
494 |
+ |
atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
495 |
+ |
} |
496 |
+ |
} |
497 |
+ |
|
498 |
+ |
|
499 |
|
for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
500 |
|
a = bond->getAtomA()->getGlobalIndex(); |
501 |
|
b = bond->getAtomB()->getGlobalIndex(); |
507 |
|
b = bend->getAtomB()->getGlobalIndex(); |
508 |
|
c = bend->getAtomC()->getGlobalIndex(); |
509 |
|
|
510 |
< |
exclude_.removePair(a, b); |
511 |
< |
exclude_.removePair(a, c); |
512 |
< |
exclude_.removePair(b, c); |
510 |
> |
std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
511 |
> |
std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
512 |
> |
std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
513 |
> |
|
514 |
> |
exclude_.removePairs(rigidSetA, rigidSetB); |
515 |
> |
exclude_.removePairs(rigidSetA, rigidSetC); |
516 |
> |
exclude_.removePairs(rigidSetB, rigidSetC); |
517 |
> |
|
518 |
> |
//exclude_.removePair(a, b); |
519 |
> |
//exclude_.removePair(a, c); |
520 |
> |
//exclude_.removePair(b, c); |
521 |
|
} |
522 |
|
|
523 |
|
for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
526 |
|
c = torsion->getAtomC()->getGlobalIndex(); |
527 |
|
d = torsion->getAtomD()->getGlobalIndex(); |
528 |
|
|
529 |
+ |
std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
530 |
+ |
std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
531 |
+ |
std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
532 |
+ |
std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
533 |
+ |
|
534 |
+ |
exclude_.removePairs(rigidSetA, rigidSetB); |
535 |
+ |
exclude_.removePairs(rigidSetA, rigidSetC); |
536 |
+ |
exclude_.removePairs(rigidSetA, rigidSetD); |
537 |
+ |
exclude_.removePairs(rigidSetB, rigidSetC); |
538 |
+ |
exclude_.removePairs(rigidSetB, rigidSetD); |
539 |
+ |
exclude_.removePairs(rigidSetC, rigidSetD); |
540 |
+ |
|
541 |
+ |
/* |
542 |
+ |
exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
543 |
+ |
exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
544 |
+ |
exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
545 |
+ |
exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
546 |
+ |
exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
547 |
+ |
exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
548 |
+ |
|
549 |
+ |
|
550 |
|
exclude_.removePair(a, b); |
551 |
|
exclude_.removePair(a, c); |
552 |
|
exclude_.removePair(a, d); |
553 |
|
exclude_.removePair(b, c); |
554 |
|
exclude_.removePair(b, d); |
555 |
|
exclude_.removePair(c, d); |
556 |
+ |
*/ |
557 |
|
} |
558 |
|
|
434 |
– |
Molecule::RigidBodyIterator rbIter; |
435 |
– |
RigidBody* rb; |
559 |
|
for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
560 |
|
std::vector<Atom*> atoms = rb->getAtoms(); |
561 |
|
for (int i = 0; i < atoms.size() -1 ; ++i) { |
595 |
|
int isError = 0; |
596 |
|
|
597 |
|
setupElectrostaticSummationMethod( isError ); |
598 |
+ |
setupSwitchingFunction(); |
599 |
|
|
600 |
|
if(isError){ |
601 |
|
sprintf( painCave.errMsg, |
640 |
|
int useLennardJones = 0; |
641 |
|
int useElectrostatic = 0; |
642 |
|
int useEAM = 0; |
643 |
+ |
int useSC = 0; |
644 |
|
int useCharge = 0; |
645 |
|
int useDirectional = 0; |
646 |
|
int useDipole = 0; |
654 |
|
//usePBC and useRF are from simParams |
655 |
|
int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
656 |
|
int useRF; |
657 |
< |
int useDW; |
657 |
> |
int useSF; |
658 |
|
std::string myMethod; |
659 |
|
|
660 |
|
// set the useRF logical |
661 |
|
useRF = 0; |
662 |
< |
useDW = 0; |
662 |
> |
useSF = 0; |
663 |
|
|
664 |
|
|
665 |
|
if (simParams_->haveElectrostaticSummationMethod()) { |
668 |
|
if (myMethod == "REACTION_FIELD") { |
669 |
|
useRF=1; |
670 |
|
} else { |
671 |
< |
if (myMethod == "SHIFTED_POTENTIAL") { |
672 |
< |
useDW = 1; |
671 |
> |
if (myMethod == "SHIFTED_FORCE") { |
672 |
> |
useSF = 1; |
673 |
|
} |
674 |
|
} |
675 |
|
} |
679 |
|
useLennardJones |= (*i)->isLennardJones(); |
680 |
|
useElectrostatic |= (*i)->isElectrostatic(); |
681 |
|
useEAM |= (*i)->isEAM(); |
682 |
+ |
useSC |= (*i)->isSC(); |
683 |
|
useCharge |= (*i)->isCharge(); |
684 |
|
useDirectional |= (*i)->isDirectional(); |
685 |
|
useDipole |= (*i)->isDipole(); |
730 |
|
temp = useEAM; |
731 |
|
MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
732 |
|
|
733 |
+ |
temp = useSC; |
734 |
+ |
MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
735 |
+ |
|
736 |
|
temp = useShape; |
737 |
|
MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
738 |
|
|
742 |
|
temp = useRF; |
743 |
|
MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
744 |
|
|
745 |
< |
temp = useDW; |
746 |
< |
MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
745 |
> |
temp = useSF; |
746 |
> |
MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
747 |
|
|
748 |
|
#endif |
749 |
|
|
757 |
|
fInfo_.SIM_uses_StickyPower = useStickyPower; |
758 |
|
fInfo_.SIM_uses_GayBerne = useGayBerne; |
759 |
|
fInfo_.SIM_uses_EAM = useEAM; |
760 |
+ |
fInfo_.SIM_uses_SC = useSC; |
761 |
|
fInfo_.SIM_uses_Shapes = useShape; |
762 |
|
fInfo_.SIM_uses_FLARB = useFLARB; |
763 |
|
fInfo_.SIM_uses_RF = useRF; |
764 |
< |
fInfo_.SIM_uses_DampedWolf = useDW; |
764 |
> |
fInfo_.SIM_uses_SF = useSF; |
765 |
|
|
766 |
|
if( myMethod == "REACTION_FIELD") { |
767 |
|
|
1075 |
|
} |
1076 |
|
} |
1077 |
|
|
1078 |
< |
if (simParams_->haveScreeningMethod()) { |
1079 |
< |
std::string myScreen = simParams_->getScreeningMethod(); |
1078 |
> |
if (simParams_->haveElectrostaticScreeningMethod()) { |
1079 |
> |
std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
1080 |
|
toUpper(myScreen); |
1081 |
|
if (myScreen == "UNDAMPED") { |
1082 |
|
sm = UNDAMPED; |
1089 |
|
"SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); |
1090 |
|
painCave.isFatal = 0; |
1091 |
|
simError(); |
962 |
– |
} else { |
963 |
– |
// throw error |
964 |
– |
sprintf( painCave.errMsg, |
965 |
– |
"SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
966 |
– |
painCave.isFatal = 1; |
967 |
– |
simError(); |
1092 |
|
} |
1093 |
+ |
} else { |
1094 |
+ |
// throw error |
1095 |
+ |
sprintf( painCave.errMsg, |
1096 |
+ |
"SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
1097 |
+ |
painCave.isFatal = 1; |
1098 |
+ |
simError(); |
1099 |
|
} |
1100 |
|
} |
1101 |
|
} |
1102 |
+ |
|
1103 |
|
// let's pass some summation method variables to fortran |
1104 |
|
setElectrostaticSummationMethod( &esm ); |
1105 |
|
setScreeningMethod( &sm ); |
1108 |
|
initFortranFF( &esm, &errorOut ); |
1109 |
|
} |
1110 |
|
|
1111 |
+ |
void SimInfo::setupSwitchingFunction() { |
1112 |
+ |
int ft = CUBIC; |
1113 |
+ |
|
1114 |
+ |
if (simParams_->haveSwitchingFunctionType()) { |
1115 |
+ |
std::string funcType = simParams_->getSwitchingFunctionType(); |
1116 |
+ |
toUpper(funcType); |
1117 |
+ |
if (funcType == "CUBIC") { |
1118 |
+ |
ft = CUBIC; |
1119 |
+ |
} else { |
1120 |
+ |
if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1121 |
+ |
ft = FIFTH_ORDER_POLY; |
1122 |
+ |
} else { |
1123 |
+ |
// throw error |
1124 |
+ |
sprintf( painCave.errMsg, |
1125 |
+ |
"SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1126 |
+ |
painCave.isFatal = 1; |
1127 |
+ |
simError(); |
1128 |
+ |
} |
1129 |
+ |
} |
1130 |
+ |
} |
1131 |
+ |
|
1132 |
+ |
// send switching function notification to switcheroo |
1133 |
+ |
setFunctionType(&ft); |
1134 |
+ |
|
1135 |
+ |
} |
1136 |
+ |
|
1137 |
|
void SimInfo::addProperty(GenericData* genData) { |
1138 |
|
properties_.addProperty(genData); |
1139 |
|
} |