| 48 |  |  | 
| 49 |  | #include <algorithm> | 
| 50 |  | #include <set> | 
| 51 | + | #include <map> | 
| 52 |  |  | 
| 53 |  | #include "brains/SimInfo.hpp" | 
| 54 |  | #include "math/Vector3.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 |  |  | 
| 379 | – | Molecule::RigidBodyIterator rbIter; | 
| 380 | – | 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)) { | 
| 525 |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 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 |  |  | 
| 436 | – | Molecule::RigidBodyIterator rbIter; | 
| 437 | – | 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) { | 
| 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; | 
| 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 |  |  | 
| 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; |