| 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 "primitives/StuntDouble.hpp" | 
| 57 |  | #include "UseTheForce/fCutoffPolicy.h" | 
| 58 |  | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" | 
| 59 | + | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" | 
| 60 | + | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" | 
| 61 |  | #include "UseTheForce/doForces_interface.h" | 
| 62 | + | #include "UseTheForce/DarkSide/neighborLists_interface.h" | 
| 63 |  | #include "UseTheForce/DarkSide/electrostatic_interface.h" | 
| 64 | < | #include "UseTheForce/notifyCutoffs_interface.h" | 
| 64 | > | #include "UseTheForce/DarkSide/switcheroo_interface.h" | 
| 65 |  | #include "utils/MemoryUtils.hpp" | 
| 66 |  | #include "utils/simError.h" | 
| 67 |  | #include "selection/SelectionManager.hpp" | 
| 68 | + | #include "io/ForceFieldOptions.hpp" | 
| 69 | + | #include "UseTheForce/ForceField.hpp" | 
| 70 |  |  | 
| 71 | + |  | 
| 72 |  | #ifdef IS_MPI | 
| 73 |  | #include "UseTheForce/mpiComponentPlan.h" | 
| 74 |  | #include "UseTheForce/DarkSide/simParallel_interface.h" | 
| 75 |  | #endif | 
| 76 |  |  | 
| 77 |  | namespace oopse { | 
| 78 | + | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { | 
| 79 | + | std::map<int, std::set<int> >::iterator i = container.find(index); | 
| 80 | + | std::set<int> result; | 
| 81 | + | if (i != container.end()) { | 
| 82 | + | result = i->second; | 
| 83 | + | } | 
| 84 |  |  | 
| 85 | < | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, | 
| 86 | < | ForceField* ff, Globals* simParams) : | 
| 87 | < | stamps_(stamps), forceField_(ff), simParams_(simParams), | 
| 88 | < | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 85 | > | return result; | 
| 86 | > | } | 
| 87 | > |  | 
| 88 | > | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | 
| 89 | > | forceField_(ff), simParams_(simParams), | 
| 90 | > | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 91 |  | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | 
| 92 |  | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 93 |  | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0), | 
| 94 |  | nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0), | 
| 95 | < | sman_(NULL), fortranInitialized_(false) { | 
| 95 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { | 
| 96 |  |  | 
| 81 | – |  | 
| 82 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; | 
| 97 |  | MoleculeStamp* molStamp; | 
| 98 |  | int nMolWithSameStamp; | 
| 99 |  | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 101 |  | CutoffGroupStamp* cgStamp; | 
| 102 |  | RigidBodyStamp* rbStamp; | 
| 103 |  | int nRigidAtoms = 0; | 
| 104 | < |  | 
| 105 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { | 
| 106 | < | molStamp = i->first; | 
| 107 | < | nMolWithSameStamp = i->second; | 
| 104 | > | std::vector<Component*> components = simParams->getComponents(); | 
| 105 | > |  | 
| 106 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 107 | > | molStamp = (*i)->getMoleculeStamp(); | 
| 108 | > | nMolWithSameStamp = (*i)->getNMol(); | 
| 109 |  |  | 
| 110 |  | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 111 |  |  | 
| 112 |  | //calculate atoms in molecules | 
| 113 |  | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 114 |  |  | 
| 100 | – |  | 
| 115 |  | //calculate atoms in cutoff groups | 
| 116 |  | int nAtomsInGroups = 0; | 
| 117 |  | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 118 |  |  | 
| 119 |  | for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 120 | < | cgStamp = molStamp->getCutoffGroup(j); | 
| 120 | > | cgStamp = molStamp->getCutoffGroupStamp(j); | 
| 121 |  | nAtomsInGroups += cgStamp->getNMembers(); | 
| 122 |  | } | 
| 123 |  |  | 
| 130 |  | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 131 |  |  | 
| 132 |  | for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 133 | < | rbStamp = molStamp->getRigidBody(j); | 
| 133 | > | rbStamp = molStamp->getRigidBodyStamp(j); | 
| 134 |  | nAtomsInRigidBodies += rbStamp->getNMembers(); | 
| 135 |  | } | 
| 136 |  |  | 
| 169 |  | } | 
| 170 |  | molecules_.clear(); | 
| 171 |  |  | 
| 158 | – | delete stamps_; | 
| 172 |  | delete sman_; | 
| 173 |  | delete simParams_; | 
| 174 |  | delete forceField_; | 
| 275 |  | } | 
| 276 |  | } | 
| 277 |  |  | 
| 278 | < | }//end for (integrableObject) | 
| 279 | < | }// end for (mol) | 
| 278 | > | } | 
| 279 | > | } | 
| 280 |  |  | 
| 281 |  | // n_constraints is local, so subtract them on each processor | 
| 282 |  | ndf_local -= nConstraints_; | 
| 293 |  |  | 
| 294 |  | } | 
| 295 |  |  | 
| 296 | + | int SimInfo::getFdf() { | 
| 297 | + | #ifdef IS_MPI | 
| 298 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 299 | + | #else | 
| 300 | + | fdf_ = fdf_local; | 
| 301 | + | #endif | 
| 302 | + | return fdf_; | 
| 303 | + | } | 
| 304 | + |  | 
| 305 |  | void SimInfo::calcNdfRaw() { | 
| 306 |  | int ndfRaw_local; | 
| 307 |  |  | 
| 364 |  | int b; | 
| 365 |  | int c; | 
| 366 |  | int d; | 
| 367 | + |  | 
| 368 | + | std::map<int, std::set<int> > atomGroups; | 
| 369 | + |  | 
| 370 | + | Molecule::RigidBodyIterator rbIter; | 
| 371 | + | RigidBody* rb; | 
| 372 | + | Molecule::IntegrableObjectIterator ii; | 
| 373 | + | StuntDouble* integrableObject; | 
| 374 |  |  | 
| 375 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 376 | + | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 377 | + |  | 
| 378 | + | if (integrableObject->isRigidBody()) { | 
| 379 | + | rb = static_cast<RigidBody*>(integrableObject); | 
| 380 | + | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 381 | + | std::set<int> rigidAtoms; | 
| 382 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 383 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 384 | + | } | 
| 385 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 386 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 387 | + | } | 
| 388 | + | } else { | 
| 389 | + | std::set<int> oneAtomSet; | 
| 390 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 391 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 392 | + | } | 
| 393 | + | } | 
| 394 | + |  | 
| 395 | + |  | 
| 396 | + |  | 
| 397 |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 398 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 399 |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 404 |  | a = bend->getAtomA()->getGlobalIndex(); | 
| 405 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 406 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 407 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 408 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 409 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 410 |  |  | 
| 411 | < | exclude_.addPair(a, b); | 
| 412 | < | exclude_.addPair(a, c); | 
| 413 | < | exclude_.addPair(b, c); | 
| 411 | > | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 412 | > | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 413 | > | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 414 | > |  | 
| 415 | > | //exclude_.addPair(a, b); | 
| 416 | > | //exclude_.addPair(a, c); | 
| 417 | > | //exclude_.addPair(b, c); | 
| 418 |  | } | 
| 419 |  |  | 
| 420 |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 422 |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 423 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 424 |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 425 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 426 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 427 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 428 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 429 |  |  | 
| 430 | + | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 431 | + | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 432 | + | exclude_.addPairs(rigidSetA, rigidSetD); | 
| 433 | + | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 434 | + | exclude_.addPairs(rigidSetB, rigidSetD); | 
| 435 | + | exclude_.addPairs(rigidSetC, rigidSetD); | 
| 436 | + |  | 
| 437 | + | /* | 
| 438 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 439 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 440 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 441 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 442 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 443 | + | exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 444 | + |  | 
| 445 | + |  | 
| 446 |  | exclude_.addPair(a, b); | 
| 447 |  | exclude_.addPair(a, c); | 
| 448 |  | exclude_.addPair(a, d); | 
| 449 |  | exclude_.addPair(b, c); | 
| 450 |  | exclude_.addPair(b, d); | 
| 451 |  | exclude_.addPair(c, d); | 
| 452 | + | */ | 
| 453 |  | } | 
| 454 |  |  | 
| 376 | – | Molecule::RigidBodyIterator rbIter; | 
| 377 | – | RigidBody* rb; | 
| 455 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 456 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 457 |  | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 476 |  | int b; | 
| 477 |  | int c; | 
| 478 |  | int d; | 
| 479 | + |  | 
| 480 | + | std::map<int, std::set<int> > atomGroups; | 
| 481 | + |  | 
| 482 | + | Molecule::RigidBodyIterator rbIter; | 
| 483 | + | RigidBody* rb; | 
| 484 | + | Molecule::IntegrableObjectIterator ii; | 
| 485 | + | StuntDouble* integrableObject; | 
| 486 |  |  | 
| 487 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 488 | + | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 489 | + |  | 
| 490 | + | if (integrableObject->isRigidBody()) { | 
| 491 | + | rb = static_cast<RigidBody*>(integrableObject); | 
| 492 | + | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 493 | + | std::set<int> rigidAtoms; | 
| 494 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 495 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 496 | + | } | 
| 497 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 498 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 499 | + | } | 
| 500 | + | } else { | 
| 501 | + | std::set<int> oneAtomSet; | 
| 502 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 503 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 504 | + | } | 
| 505 | + | } | 
| 506 | + |  | 
| 507 | + |  | 
| 508 |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 509 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 510 |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 516 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 517 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 518 |  |  | 
| 519 | < | exclude_.removePair(a, b); | 
| 520 | < | exclude_.removePair(a, c); | 
| 521 | < | exclude_.removePair(b, c); | 
| 519 | > | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 520 | > | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 521 | > | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 522 | > |  | 
| 523 | > | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 524 | > | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 525 | > | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 526 | > |  | 
| 527 | > | //exclude_.removePair(a, b); | 
| 528 | > | //exclude_.removePair(a, c); | 
| 529 | > | //exclude_.removePair(b, c); | 
| 530 |  | } | 
| 531 |  |  | 
| 532 |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 535 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 536 |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 537 |  |  | 
| 538 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 539 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 540 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 541 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 542 | + |  | 
| 543 | + | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 544 | + | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 545 | + | exclude_.removePairs(rigidSetA, rigidSetD); | 
| 546 | + | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 547 | + | exclude_.removePairs(rigidSetB, rigidSetD); | 
| 548 | + | exclude_.removePairs(rigidSetC, rigidSetD); | 
| 549 | + |  | 
| 550 | + | /* | 
| 551 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 552 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 553 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 554 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 555 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 556 | + | exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 557 | + |  | 
| 558 | + |  | 
| 559 |  | exclude_.removePair(a, b); | 
| 560 |  | exclude_.removePair(a, c); | 
| 561 |  | exclude_.removePair(a, d); | 
| 562 |  | exclude_.removePair(b, c); | 
| 563 |  | exclude_.removePair(b, d); | 
| 564 |  | exclude_.removePair(c, d); | 
| 565 | + | */ | 
| 566 |  | } | 
| 567 |  |  | 
| 433 | – | Molecule::RigidBodyIterator rbIter; | 
| 434 | – | RigidBody* rb; | 
| 568 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 569 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 570 |  | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 603 |  | /** @deprecate */ | 
| 604 |  | int isError = 0; | 
| 605 |  |  | 
| 606 | + | setupCutoff(); | 
| 607 | + |  | 
| 608 |  | setupElectrostaticSummationMethod( isError ); | 
| 609 | + | setupSwitchingFunction(); | 
| 610 | + | setupAccumulateBoxDipole(); | 
| 611 |  |  | 
| 612 |  | if(isError){ | 
| 613 |  | sprintf( painCave.errMsg, | 
| 615 |  | painCave.isFatal = 1; | 
| 616 |  | simError(); | 
| 617 |  | } | 
| 481 | – |  | 
| 482 | – |  | 
| 483 | – | setupCutoff(); | 
| 618 |  |  | 
| 619 |  | calcNdf(); | 
| 620 |  | calcNdfRaw(); | 
| 649 |  | int useLennardJones = 0; | 
| 650 |  | int useElectrostatic = 0; | 
| 651 |  | int useEAM = 0; | 
| 652 | + | int useSC = 0; | 
| 653 |  | int useCharge = 0; | 
| 654 |  | int useDirectional = 0; | 
| 655 |  | int useDipole = 0; | 
| 663 |  | //usePBC and useRF are from simParams | 
| 664 |  | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | 
| 665 |  | int useRF; | 
| 666 | + | int useSF; | 
| 667 | + | int useSP; | 
| 668 | + | int useBoxDipole; | 
| 669 |  | std::string myMethod; | 
| 670 |  |  | 
| 671 |  | // set the useRF logical | 
| 672 |  | useRF = 0; | 
| 673 | + | useSF = 0; | 
| 674 | + | useSP = 0; | 
| 675 | + |  | 
| 676 | + |  | 
| 677 |  | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 678 | < | myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 679 | < | if (myMethod == "REACTION_FIELD") | 
| 680 | < | useRF = 1; | 
| 678 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 679 | > | toUpper(myMethod); | 
| 680 | > | if (myMethod == "REACTION_FIELD"){ | 
| 681 | > | useRF = 1; | 
| 682 | > | } else if (myMethod == "SHIFTED_FORCE"){ | 
| 683 | > | useSF = 1; | 
| 684 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ | 
| 685 | > | useSP = 1; | 
| 686 | > | } | 
| 687 |  | } | 
| 688 | + |  | 
| 689 | + | if (simParams_->haveAccumulateBoxDipole()) | 
| 690 | + | if (simParams_->getAccumulateBoxDipole()) | 
| 691 | + | useBoxDipole = 1; | 
| 692 |  |  | 
| 693 |  | //loop over all of the atom types | 
| 694 |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 695 |  | useLennardJones |= (*i)->isLennardJones(); | 
| 696 |  | useElectrostatic |= (*i)->isElectrostatic(); | 
| 697 |  | useEAM |= (*i)->isEAM(); | 
| 698 | + | useSC |= (*i)->isSC(); | 
| 699 |  | useCharge |= (*i)->isCharge(); | 
| 700 |  | useDirectional |= (*i)->isDirectional(); | 
| 701 |  | useDipole |= (*i)->isDipole(); | 
| 746 |  | temp = useEAM; | 
| 747 |  | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 748 |  |  | 
| 749 | + | temp = useSC; | 
| 750 | + | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 751 | + |  | 
| 752 |  | temp = useShape; | 
| 753 |  | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 754 |  |  | 
| 758 |  | temp = useRF; | 
| 759 |  | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 760 |  |  | 
| 761 | + | temp = useSF; | 
| 762 | + | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 763 | + |  | 
| 764 | + | temp = useSP; | 
| 765 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 766 | + |  | 
| 767 | + | temp = useBoxDipole; | 
| 768 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 769 | + |  | 
| 770 |  | #endif | 
| 771 |  |  | 
| 772 |  | fInfo_.SIM_uses_PBC = usePBC; | 
| 779 |  | fInfo_.SIM_uses_StickyPower = useStickyPower; | 
| 780 |  | fInfo_.SIM_uses_GayBerne = useGayBerne; | 
| 781 |  | fInfo_.SIM_uses_EAM = useEAM; | 
| 782 | + | fInfo_.SIM_uses_SC = useSC; | 
| 783 |  | fInfo_.SIM_uses_Shapes = useShape; | 
| 784 |  | fInfo_.SIM_uses_FLARB = useFLARB; | 
| 785 |  | fInfo_.SIM_uses_RF = useRF; | 
| 786 | < |  | 
| 787 | < | if( fInfo_.SIM_uses_Dipoles && myMethod == "REACTION_FIELD") { | 
| 788 | < |  | 
| 623 | < | if (simParams_->haveDielectric()) { | 
| 624 | < | fInfo_.dielect = simParams_->getDielectric(); | 
| 625 | < | } else { | 
| 626 | < | sprintf(painCave.errMsg, | 
| 627 | < | "SimSetup Error: No Dielectric constant was set.\n" | 
| 628 | < | "\tYou are trying to use Reaction Field without" | 
| 629 | < | "\tsetting a dielectric constant!\n"); | 
| 630 | < | painCave.isFatal = 1; | 
| 631 | < | simError(); | 
| 632 | < | } | 
| 633 | < |  | 
| 634 | < | } else { | 
| 635 | < | fInfo_.dielect = 0.0; | 
| 636 | < | } | 
| 637 | < |  | 
| 786 | > | fInfo_.SIM_uses_SF = useSF; | 
| 787 | > | fInfo_.SIM_uses_SP = useSP; | 
| 788 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; | 
| 789 |  | } | 
| 790 |  |  | 
| 791 |  | void SimInfo::setupFortranSim() { | 
| 802 |  | } | 
| 803 |  |  | 
| 804 |  | //calculate mass ratio of cutoff group | 
| 805 | < | std::vector<double> mfact; | 
| 805 | > | std::vector<RealType> mfact; | 
| 806 |  | SimInfo::MoleculeIterator mi; | 
| 807 |  | Molecule* mol; | 
| 808 |  | Molecule::CutoffGroupIterator ci; | 
| 809 |  | CutoffGroup* cg; | 
| 810 |  | Molecule::AtomIterator ai; | 
| 811 |  | Atom* atom; | 
| 812 | < | double totalMass; | 
| 812 | > | RealType totalMass; | 
| 813 |  |  | 
| 814 |  | //to avoid memory reallocation, reserve enough space for mfact | 
| 815 |  | mfact.reserve(getNCutoffGroups()); | 
| 870 |  | "succesfully sent the simulation information to fortran.\n"); | 
| 871 |  | MPIcheckPoint(); | 
| 872 |  | #endif // is_mpi | 
| 873 | + |  | 
| 874 | + | // Setup number of neighbors in neighbor list if present | 
| 875 | + | if (simParams_->haveNeighborListNeighbors()) { | 
| 876 | + | setNeighbors(simParams_->getNeighborListNeighbors()); | 
| 877 | + | } | 
| 878 | + |  | 
| 879 | + |  | 
| 880 |  | } | 
| 881 |  |  | 
| 882 |  |  | 
| 939 |  |  | 
| 940 |  | #endif | 
| 941 |  |  | 
| 942 | < | double SimInfo::calcMaxCutoffRadius() { | 
| 942 | > | void SimInfo::setupCutoff() { | 
| 943 | > |  | 
| 944 | > | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | 
| 945 |  |  | 
| 946 | + | // Check the cutoff policy | 
| 947 | + | int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | 
| 948 |  |  | 
| 949 | < | std::set<AtomType*> atomTypes; | 
| 950 | < | std::set<AtomType*>::iterator i; | 
| 951 | < | std::vector<double> cutoffRadius; | 
| 952 | < |  | 
| 953 | < | //get the unique atom types | 
| 792 | < | atomTypes = getUniqueAtomTypes(); | 
| 793 | < |  | 
| 794 | < | //query the max cutoff radius among these atom types | 
| 795 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 796 | < | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); | 
| 949 | > | std::string myPolicy; | 
| 950 | > | if (forceFieldOptions_.haveCutoffPolicy()){ | 
| 951 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); | 
| 952 | > | }else if (simParams_->haveCutoffPolicy()) { | 
| 953 | > | myPolicy = simParams_->getCutoffPolicy(); | 
| 954 |  | } | 
| 955 |  |  | 
| 956 | < | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); | 
| 800 | < | #ifdef IS_MPI | 
| 801 | < | //pick the max cutoff radius among the processors | 
| 802 | < | #endif | 
| 803 | < |  | 
| 804 | < | return maxCutoffRadius; | 
| 805 | < | } | 
| 806 | < |  | 
| 807 | < | void SimInfo::getCutoff(double& rcut, double& rsw) { | 
| 808 | < |  | 
| 809 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 810 | < |  | 
| 811 | < | if (!simParams_->haveCutoffRadius()){ | 
| 812 | < | sprintf(painCave.errMsg, | 
| 813 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 814 | < | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 815 | < | "\tfor the cutoffRadius.\n"); | 
| 816 | < | painCave.isFatal = 0; | 
| 817 | < | simError(); | 
| 818 | < | rcut = 15.0; | 
| 819 | < | } else{ | 
| 820 | < | rcut = simParams_->getCutoffRadius(); | 
| 821 | < | } | 
| 822 | < |  | 
| 823 | < | if (!simParams_->haveSwitchingRadius()){ | 
| 824 | < | sprintf(painCave.errMsg, | 
| 825 | < | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 826 | < | "\tOOPSE will use a default value of\n" | 
| 827 | < | "\t0.95 * cutoffRadius for the switchingRadius\n"); | 
| 828 | < | painCave.isFatal = 0; | 
| 829 | < | simError(); | 
| 830 | < | rsw = 0.95 * rcut; | 
| 831 | < | } else{ | 
| 832 | < | rsw = simParams_->getSwitchingRadius(); | 
| 833 | < | } | 
| 834 | < |  | 
| 835 | < | } else { | 
| 836 | < | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | 
| 837 | < | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | 
| 838 | < |  | 
| 839 | < | if (simParams_->haveCutoffRadius()) { | 
| 840 | < | rcut = simParams_->getCutoffRadius(); | 
| 841 | < | } else { | 
| 842 | < | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | 
| 843 | < | rcut = calcMaxCutoffRadius(); | 
| 844 | < | } | 
| 845 | < |  | 
| 846 | < | if (simParams_->haveSwitchingRadius()) { | 
| 847 | < | rsw  = simParams_->getSwitchingRadius(); | 
| 848 | < | } else { | 
| 849 | < | rsw = rcut; | 
| 850 | < | } | 
| 851 | < |  | 
| 852 | < | } | 
| 853 | < | } | 
| 854 | < |  | 
| 855 | < | void SimInfo::setupCutoff() { | 
| 856 | < | getCutoff(rcut_, rsw_); | 
| 857 | < | double rnblist = rcut_ + 1; // skin of neighbor list | 
| 858 | < |  | 
| 859 | < | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | 
| 860 | < |  | 
| 861 | < | int cp =  TRADITIONAL_CUTOFF_POLICY; | 
| 862 | < | if (simParams_->haveCutoffPolicy()) { | 
| 863 | < | std::string myPolicy = simParams_->getCutoffPolicy(); | 
| 956 | > | if (!myPolicy.empty()){ | 
| 957 |  | toUpper(myPolicy); | 
| 958 |  | if (myPolicy == "MIX") { | 
| 959 |  | cp = MIX_CUTOFF_POLICY; | 
| 972 |  | } | 
| 973 |  | } | 
| 974 |  | } | 
| 975 | < | } | 
| 976 | < |  | 
| 975 | > | } | 
| 976 | > | notifyFortranCutoffPolicy(&cp); | 
| 977 |  |  | 
| 978 | + | // Check the Skin Thickness for neighborlists | 
| 979 | + | RealType skin; | 
| 980 |  | if (simParams_->haveSkinThickness()) { | 
| 981 | < | double skinThickness = simParams_->getSkinThickness(); | 
| 982 | < | } | 
| 981 | > | skin = simParams_->getSkinThickness(); | 
| 982 | > | notifyFortranSkinThickness(&skin); | 
| 983 | > | } | 
| 984 | > |  | 
| 985 | > | // Check if the cutoff was set explicitly: | 
| 986 | > | if (simParams_->haveCutoffRadius()) { | 
| 987 | > | rcut_ = simParams_->getCutoffRadius(); | 
| 988 | > | if (simParams_->haveSwitchingRadius()) { | 
| 989 | > | rsw_  = simParams_->getSwitchingRadius(); | 
| 990 | > | } else { | 
| 991 | > | if (fInfo_.SIM_uses_Charges | | 
| 992 | > | fInfo_.SIM_uses_Dipoles | | 
| 993 | > | fInfo_.SIM_uses_RF) { | 
| 994 | > |  | 
| 995 | > | rsw_ = 0.85 * rcut_; | 
| 996 | > | sprintf(painCave.errMsg, | 
| 997 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 998 | > | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" | 
| 999 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1000 | > | painCave.isFatal = 0; | 
| 1001 | > | simError(); | 
| 1002 | > | } else { | 
| 1003 | > | rsw_ = rcut_; | 
| 1004 | > | sprintf(painCave.errMsg, | 
| 1005 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 1006 | > | "\tOOPSE will use the same value as the cutoffRadius.\n" | 
| 1007 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1008 | > | painCave.isFatal = 0; | 
| 1009 | > | simError(); | 
| 1010 | > | } | 
| 1011 | > | } | 
| 1012 | > |  | 
| 1013 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1014 | > |  | 
| 1015 | > | } else { | 
| 1016 | > |  | 
| 1017 | > | // For electrostatic atoms, we'll assume a large safe value: | 
| 1018 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 1019 | > | sprintf(painCave.errMsg, | 
| 1020 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 1021 | > | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 1022 | > | "\tfor the cutoffRadius.\n"); | 
| 1023 | > | painCave.isFatal = 0; | 
| 1024 | > | simError(); | 
| 1025 | > | rcut_ = 15.0; | 
| 1026 | > |  | 
| 1027 | > | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1028 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1029 | > | toUpper(myMethod); | 
| 1030 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 1031 | > | if (simParams_->haveSwitchingRadius()){ | 
| 1032 | > | sprintf(painCave.errMsg, | 
| 1033 | > | "SimInfo Warning: A value was set for the switchingRadius\n" | 
| 1034 | > | "\teven though the electrostaticSummationMethod was\n" | 
| 1035 | > | "\tset to %s\n", myMethod.c_str()); | 
| 1036 | > | painCave.isFatal = 1; | 
| 1037 | > | simError(); | 
| 1038 | > | } | 
| 1039 | > | } | 
| 1040 | > | } | 
| 1041 | > |  | 
| 1042 | > | if (simParams_->haveSwitchingRadius()){ | 
| 1043 | > | rsw_ = simParams_->getSwitchingRadius(); | 
| 1044 | > | } else { | 
| 1045 | > | sprintf(painCave.errMsg, | 
| 1046 | > | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 1047 | > | "\tOOPSE will use a default value of\n" | 
| 1048 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 1049 | > | painCave.isFatal = 0; | 
| 1050 | > | simError(); | 
| 1051 | > | rsw_ = 0.85 * rcut_; | 
| 1052 | > | } | 
| 1053 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1054 | > | } else { | 
| 1055 | > | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | 
| 1056 | > | // We'll punt and let fortran figure out the cutoffs later. | 
| 1057 | > |  | 
| 1058 | > | notifyFortranYouAreOnYourOwn(); | 
| 1059 |  |  | 
| 1060 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | 
| 1061 | < | // also send cutoff notification to electrostatics | 
| 891 | < | setElectrostaticCutoffRadius(&rcut_); | 
| 1060 | > | } | 
| 1061 | > | } | 
| 1062 |  | } | 
| 1063 |  |  | 
| 1064 |  | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | 
| 1065 |  |  | 
| 1066 |  | int errorOut; | 
| 1067 |  | int esm =  NONE; | 
| 1068 | < | double alphaVal; | 
| 1069 | < | double dielectric; | 
| 1070 | < |  | 
| 1068 | > | int sm = UNDAMPED; | 
| 1069 | > | RealType alphaVal; | 
| 1070 | > | RealType dielectric; | 
| 1071 | > |  | 
| 1072 |  | errorOut = isError; | 
| 902 | – | alphaVal = simParams_->getDampingAlpha(); | 
| 903 | – | dielectric = simParams_->getDielectric(); | 
| 1073 |  |  | 
| 1074 |  | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1075 |  | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1077 |  | if (myMethod == "NONE") { | 
| 1078 |  | esm = NONE; | 
| 1079 |  | } else { | 
| 1080 | < | if (myMethod == "UNDAMPED_WOLF") { | 
| 1081 | < | esm = UNDAMPED_WOLF; | 
| 1080 | > | if (myMethod == "SWITCHING_FUNCTION") { | 
| 1081 | > | esm = SWITCHING_FUNCTION; | 
| 1082 |  | } else { | 
| 1083 | < | if (myMethod == "DAMPED_WOLF") { | 
| 1084 | < | esm = DAMPED_WOLF; | 
| 1085 | < | if (!simParams_->haveDampingAlpha()) { | 
| 1086 | < | //throw error | 
| 1087 | < | sprintf( painCave.errMsg, | 
| 919 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal); | 
| 920 | < | painCave.isFatal = 0; | 
| 921 | < | simError(); | 
| 922 | < | } | 
| 923 | < | } else { | 
| 924 | < | if (myMethod == "REACTION_FIELD") { | 
| 925 | < | esm = REACTION_FIELD; | 
| 1083 | > | if (myMethod == "SHIFTED_POTENTIAL") { | 
| 1084 | > | esm = SHIFTED_POTENTIAL; | 
| 1085 | > | } else { | 
| 1086 | > | if (myMethod == "SHIFTED_FORCE") { | 
| 1087 | > | esm = SHIFTED_FORCE; | 
| 1088 |  | } else { | 
| 1089 | < | // throw error | 
| 1090 | < | sprintf( painCave.errMsg, | 
| 1091 | < | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() ); | 
| 1092 | < | painCave.isFatal = 1; | 
| 1093 | < | simError(); | 
| 1094 | < | } | 
| 1095 | < | } | 
| 1089 | > | if (myMethod == "REACTION_FIELD") { | 
| 1090 | > | esm = REACTION_FIELD; | 
| 1091 | > | dielectric = simParams_->getDielectric(); | 
| 1092 | > | if (!simParams_->haveDielectric()) { | 
| 1093 | > | // throw warning | 
| 1094 | > | sprintf( painCave.errMsg, | 
| 1095 | > | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" | 
| 1096 | > | "\tA default value of %f will be used for the dielectric.\n", dielectric); | 
| 1097 | > | painCave.isFatal = 0; | 
| 1098 | > | simError(); | 
| 1099 | > | } | 
| 1100 | > | } else { | 
| 1101 | > | // throw error | 
| 1102 | > | sprintf( painCave.errMsg, | 
| 1103 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" | 
| 1104 | > | "\t(Input file specified %s .)\n" | 
| 1105 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" | 
| 1106 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" | 
| 1107 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); | 
| 1108 | > | painCave.isFatal = 1; | 
| 1109 | > | simError(); | 
| 1110 | > | } | 
| 1111 | > | } | 
| 1112 | > | } | 
| 1113 | > | } | 
| 1114 | > | } | 
| 1115 | > | } | 
| 1116 | > |  | 
| 1117 | > | if (simParams_->haveElectrostaticScreeningMethod()) { | 
| 1118 | > | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); | 
| 1119 | > | toUpper(myScreen); | 
| 1120 | > | if (myScreen == "UNDAMPED") { | 
| 1121 | > | sm = UNDAMPED; | 
| 1122 | > | } else { | 
| 1123 | > | if (myScreen == "DAMPED") { | 
| 1124 | > | sm = DAMPED; | 
| 1125 | > | if (!simParams_->haveDampingAlpha()) { | 
| 1126 | > | // first set a cutoff dependent alpha value | 
| 1127 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang | 
| 1128 | > | alphaVal = 0.5125 - rcut_* 0.025; | 
| 1129 | > | // for values rcut > 20.5, alpha is zero | 
| 1130 | > | if (alphaVal < 0) alphaVal = 0; | 
| 1131 | > |  | 
| 1132 | > | // throw warning | 
| 1133 | > | sprintf( painCave.errMsg, | 
| 1134 | > | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1135 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); | 
| 1136 | > | painCave.isFatal = 0; | 
| 1137 | > | simError(); | 
| 1138 | > | } else { | 
| 1139 | > | alphaVal = simParams_->getDampingAlpha(); | 
| 1140 | > | } | 
| 1141 | > |  | 
| 1142 | > | } else { | 
| 1143 | > | // throw error | 
| 1144 | > | sprintf( painCave.errMsg, | 
| 1145 | > | "SimInfo error: Unknown electrostaticScreeningMethod.\n" | 
| 1146 | > | "\t(Input file specified %s .)\n" | 
| 1147 | > | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" | 
| 1148 | > | "or \"damped\".\n", myScreen.c_str() ); | 
| 1149 | > | painCave.isFatal = 1; | 
| 1150 | > | simError(); | 
| 1151 |  | } | 
| 1152 |  | } | 
| 1153 |  | } | 
| 1154 | + |  | 
| 1155 |  | // let's pass some summation method variables to fortran | 
| 1156 |  | setElectrostaticSummationMethod( &esm ); | 
| 1157 | < | setDampedWolfAlpha( &alphaVal ); | 
| 1157 | > | setFortranElectrostaticMethod( &esm ); | 
| 1158 | > | setScreeningMethod( &sm ); | 
| 1159 | > | setDampingAlpha( &alphaVal ); | 
| 1160 |  | setReactionFieldDielectric( &dielectric ); | 
| 1161 | < | initFortranFF( &esm, &errorOut ); | 
| 1161 | > | initFortranFF( &errorOut ); | 
| 1162 |  | } | 
| 1163 |  |  | 
| 1164 | + | void SimInfo::setupSwitchingFunction() { | 
| 1165 | + | int ft = CUBIC; | 
| 1166 | + |  | 
| 1167 | + | if (simParams_->haveSwitchingFunctionType()) { | 
| 1168 | + | std::string funcType = simParams_->getSwitchingFunctionType(); | 
| 1169 | + | toUpper(funcType); | 
| 1170 | + | if (funcType == "CUBIC") { | 
| 1171 | + | ft = CUBIC; | 
| 1172 | + | } else { | 
| 1173 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { | 
| 1174 | + | ft = FIFTH_ORDER_POLY; | 
| 1175 | + | } else { | 
| 1176 | + | // throw error | 
| 1177 | + | sprintf( painCave.errMsg, | 
| 1178 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); | 
| 1179 | + | painCave.isFatal = 1; | 
| 1180 | + | simError(); | 
| 1181 | + | } | 
| 1182 | + | } | 
| 1183 | + | } | 
| 1184 | + |  | 
| 1185 | + | // send switching function notification to switcheroo | 
| 1186 | + | setFunctionType(&ft); | 
| 1187 | + |  | 
| 1188 | + | } | 
| 1189 | + |  | 
| 1190 | + | void SimInfo::setupAccumulateBoxDipole() { | 
| 1191 | + |  | 
| 1192 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 1193 | + | if ( simParams_->haveAccumulateBoxDipole() ) | 
| 1194 | + | if ( simParams_->getAccumulateBoxDipole() ) { | 
| 1195 | + | setAccumulateBoxDipole(); | 
| 1196 | + | calcBoxDipole_ = true; | 
| 1197 | + | } | 
| 1198 | + |  | 
| 1199 | + | } | 
| 1200 | + |  | 
| 1201 |  | void SimInfo::addProperty(GenericData* genData) { | 
| 1202 |  | properties_.addProperty(genData); | 
| 1203 |  | } | 
| 1254 |  | Molecule* mol; | 
| 1255 |  |  | 
| 1256 |  | Vector3d comVel(0.0); | 
| 1257 | < | double totalMass = 0.0; | 
| 1257 | > | RealType totalMass = 0.0; | 
| 1258 |  |  | 
| 1259 |  |  | 
| 1260 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1261 | < | double mass = mol->getMass(); | 
| 1261 | > | RealType mass = mol->getMass(); | 
| 1262 |  | totalMass += mass; | 
| 1263 |  | comVel += mass * mol->getComVel(); | 
| 1264 |  | } | 
| 1265 |  |  | 
| 1266 |  | #ifdef IS_MPI | 
| 1267 | < | double tmpMass = totalMass; | 
| 1267 | > | RealType tmpMass = totalMass; | 
| 1268 |  | Vector3d tmpComVel(comVel); | 
| 1269 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1270 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1269 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1270 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1271 |  | #endif | 
| 1272 |  |  | 
| 1273 |  | comVel /= totalMass; | 
| 1280 |  | Molecule* mol; | 
| 1281 |  |  | 
| 1282 |  | Vector3d com(0.0); | 
| 1283 | < | double totalMass = 0.0; | 
| 1283 | > | RealType totalMass = 0.0; | 
| 1284 |  |  | 
| 1285 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1286 | < | double mass = mol->getMass(); | 
| 1286 | > | RealType mass = mol->getMass(); | 
| 1287 |  | totalMass += mass; | 
| 1288 |  | com += mass * mol->getCom(); | 
| 1289 |  | } | 
| 1290 |  |  | 
| 1291 |  | #ifdef IS_MPI | 
| 1292 | < | double tmpMass = totalMass; | 
| 1292 | > | RealType tmpMass = totalMass; | 
| 1293 |  | Vector3d tmpCom(com); | 
| 1294 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1295 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1294 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1295 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1296 |  | #endif | 
| 1297 |  |  | 
| 1298 |  | com /= totalMass; | 
| 1316 |  | Molecule* mol; | 
| 1317 |  |  | 
| 1318 |  |  | 
| 1319 | < | double totalMass = 0.0; | 
| 1319 | > | RealType totalMass = 0.0; | 
| 1320 |  |  | 
| 1321 |  |  | 
| 1322 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1323 | < | double mass = mol->getMass(); | 
| 1323 | > | RealType mass = mol->getMass(); | 
| 1324 |  | totalMass += mass; | 
| 1325 |  | com += mass * mol->getCom(); | 
| 1326 |  | comVel += mass * mol->getComVel(); | 
| 1327 |  | } | 
| 1328 |  |  | 
| 1329 |  | #ifdef IS_MPI | 
| 1330 | < | double tmpMass = totalMass; | 
| 1330 | > | RealType tmpMass = totalMass; | 
| 1331 |  | Vector3d tmpCom(com); | 
| 1332 |  | Vector3d tmpComVel(comVel); | 
| 1333 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1334 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1335 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1333 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1334 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1335 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1336 |  | #endif | 
| 1337 |  |  | 
| 1338 |  | com /= totalMass; | 
| 1351 |  | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | 
| 1352 |  |  | 
| 1353 |  |  | 
| 1354 | < | double xx = 0.0; | 
| 1355 | < | double yy = 0.0; | 
| 1356 | < | double zz = 0.0; | 
| 1357 | < | double xy = 0.0; | 
| 1358 | < | double xz = 0.0; | 
| 1359 | < | double yz = 0.0; | 
| 1354 | > | RealType xx = 0.0; | 
| 1355 | > | RealType yy = 0.0; | 
| 1356 | > | RealType zz = 0.0; | 
| 1357 | > | RealType xy = 0.0; | 
| 1358 | > | RealType xz = 0.0; | 
| 1359 | > | RealType yz = 0.0; | 
| 1360 |  | Vector3d com(0.0); | 
| 1361 |  | Vector3d comVel(0.0); | 
| 1362 |  |  | 
| 1368 |  | Vector3d thisq(0.0); | 
| 1369 |  | Vector3d thisv(0.0); | 
| 1370 |  |  | 
| 1371 | < | double thisMass = 0.0; | 
| 1371 | > | RealType thisMass = 0.0; | 
| 1372 |  |  | 
| 1373 |  |  | 
| 1374 |  |  | 
| 1406 |  | #ifdef IS_MPI | 
| 1407 |  | Mat3x3d tmpI(inertiaTensor); | 
| 1408 |  | Vector3d tmpAngMom; | 
| 1409 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1410 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1409 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1410 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1411 |  | #endif | 
| 1412 |  |  | 
| 1413 |  | return; | 
| 1428 |  | Vector3d thisr(0.0); | 
| 1429 |  | Vector3d thisp(0.0); | 
| 1430 |  |  | 
| 1431 | < | double thisMass; | 
| 1431 | > | RealType thisMass; | 
| 1432 |  |  | 
| 1433 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1434 |  | thisMass = mol->getMass(); | 
| 1441 |  |  | 
| 1442 |  | #ifdef IS_MPI | 
| 1443 |  | Vector3d tmpAngMom; | 
| 1444 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1444 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1445 |  | #endif | 
| 1446 |  |  | 
| 1447 |  | return angularMomentum; | 
| 1448 |  | } | 
| 1449 |  |  | 
| 1450 | < |  | 
| 1450 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | 
| 1451 | > | return IOIndexToIntegrableObject.at(index); | 
| 1452 | > | } | 
| 1453 | > |  | 
| 1454 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { | 
| 1455 | > | IOIndexToIntegrableObject= v; | 
| 1456 | > | } | 
| 1457 | > |  | 
| 1458 | > | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes | 
| 1459 | > | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 | 
| 1460 | > | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to | 
| 1461 | > | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. | 
| 1462 | > | */ | 
| 1463 | > | void SimInfo::getGyrationalVolume(RealType &volume){ | 
| 1464 | > | Mat3x3d intTensor; | 
| 1465 | > | RealType det; | 
| 1466 | > | Vector3d dummyAngMom; | 
| 1467 | > | RealType sysconstants; | 
| 1468 | > | RealType geomCnst; | 
| 1469 | > |  | 
| 1470 | > | geomCnst = 3.0/2.0; | 
| 1471 | > | /* Get the inertial tensor and angular momentum for free*/ | 
| 1472 | > | getInertiaTensor(intTensor,dummyAngMom); | 
| 1473 | > |  | 
| 1474 | > | det = intTensor.determinant(); | 
| 1475 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | 
| 1476 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); | 
| 1477 | > | return; | 
| 1478 | > | } | 
| 1479 | > |  | 
| 1480 | > | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ | 
| 1481 | > | Mat3x3d intTensor; | 
| 1482 | > | Vector3d dummyAngMom; | 
| 1483 | > | RealType sysconstants; | 
| 1484 | > | RealType geomCnst; | 
| 1485 | > |  | 
| 1486 | > | geomCnst = 3.0/2.0; | 
| 1487 | > | /* Get the inertial tensor and angular momentum for free*/ | 
| 1488 | > | getInertiaTensor(intTensor,dummyAngMom); | 
| 1489 | > |  | 
| 1490 | > | detI = intTensor.determinant(); | 
| 1491 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | 
| 1492 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); | 
| 1493 | > | return; | 
| 1494 | > | } | 
| 1495 | > | /* | 
| 1496 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { | 
| 1497 | > | assert( v.size() == nAtoms_ + nRigidBodies_); | 
| 1498 | > | sdByGlobalIndex_ = v; | 
| 1499 | > | } | 
| 1500 | > |  | 
| 1501 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { | 
| 1502 | > | //assert(index < nAtoms_ + nRigidBodies_); | 
| 1503 | > | return sdByGlobalIndex_.at(index); | 
| 1504 | > | } | 
| 1505 | > | */ | 
| 1506 |  | }//end namespace oopse | 
| 1507 |  |  |