| 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/notifyCutoffs_interface.h" | 
| 62 | > | #include "UseTheForce/DarkSide/electrostatic_interface.h" | 
| 63 | > | #include "UseTheForce/DarkSide/switcheroo_interface.h" | 
| 64 |  | #include "utils/MemoryUtils.hpp" | 
| 65 |  | #include "utils/simError.h" | 
| 66 |  | #include "selection/SelectionManager.hpp" | 
| 67 | + | #include "io/ForceFieldOptions.hpp" | 
| 68 | + | #include "UseTheForce/ForceField.hpp" | 
| 69 |  |  | 
| 70 |  | #ifdef IS_MPI | 
| 71 |  | #include "UseTheForce/mpiComponentPlan.h" | 
| 73 |  | #endif | 
| 74 |  |  | 
| 75 |  | namespace oopse { | 
| 76 | + | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { | 
| 77 | + | std::map<int, std::set<int> >::iterator i = container.find(index); | 
| 78 | + | std::set<int> result; | 
| 79 | + | if (i != container.end()) { | 
| 80 | + | result = i->second; | 
| 81 | + | } | 
| 82 |  |  | 
| 83 | < | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, | 
| 84 | < | ForceField* ff, Globals* simParams) : | 
| 85 | < | stamps_(stamps), forceField_(ff), simParams_(simParams), | 
| 86 | < | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 83 | > | return result; | 
| 84 | > | } | 
| 85 | > |  | 
| 86 | > | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | 
| 87 | > | forceField_(ff), simParams_(simParams), | 
| 88 | > | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 89 |  | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | 
| 90 |  | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 91 |  | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0), | 
| 92 |  | nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0), | 
| 93 | < | sman_(NULL), fortranInitialized_(false) { | 
| 93 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { | 
| 94 |  |  | 
| 79 | – |  | 
| 80 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; | 
| 95 |  | MoleculeStamp* molStamp; | 
| 96 |  | int nMolWithSameStamp; | 
| 97 |  | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 98 | < | int nGroups = 0;          //total cutoff groups defined in meta-data file | 
| 98 | > | int nGroups = 0;      //total cutoff groups defined in meta-data file | 
| 99 |  | CutoffGroupStamp* cgStamp; | 
| 100 |  | RigidBodyStamp* rbStamp; | 
| 101 |  | int nRigidAtoms = 0; | 
| 102 | < |  | 
| 103 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { | 
| 104 | < | molStamp = i->first; | 
| 105 | < | nMolWithSameStamp = i->second; | 
| 102 | > | std::vector<Component*> components = simParams->getComponents(); | 
| 103 | > |  | 
| 104 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 105 | > | molStamp = (*i)->getMoleculeStamp(); | 
| 106 | > | nMolWithSameStamp = (*i)->getNMol(); | 
| 107 |  |  | 
| 108 |  | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 109 |  |  | 
| 110 |  | //calculate atoms in molecules | 
| 111 |  | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 112 |  |  | 
| 98 | – |  | 
| 113 |  | //calculate atoms in cutoff groups | 
| 114 |  | int nAtomsInGroups = 0; | 
| 115 |  | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 116 |  |  | 
| 117 |  | for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 118 | < | cgStamp = molStamp->getCutoffGroup(j); | 
| 118 | > | cgStamp = molStamp->getCutoffGroupStamp(j); | 
| 119 |  | nAtomsInGroups += cgStamp->getNMembers(); | 
| 120 |  | } | 
| 121 |  |  | 
| 122 |  | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | 
| 123 | + |  | 
| 124 |  | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | 
| 125 |  |  | 
| 126 |  | //calculate atoms in rigid bodies | 
| 128 |  | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 129 |  |  | 
| 130 |  | for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 131 | < | rbStamp = molStamp->getRigidBody(j); | 
| 131 | > | rbStamp = molStamp->getRigidBodyStamp(j); | 
| 132 |  | nAtomsInRigidBodies += rbStamp->getNMembers(); | 
| 133 |  | } | 
| 134 |  |  | 
| 137 |  |  | 
| 138 |  | } | 
| 139 |  |  | 
| 140 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff group | 
| 141 | < | //therefore the total number of cutoff groups in the system is equal to | 
| 142 | < | //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data | 
| 143 | < | //file plus the number of cutoff groups defined in meta-data file | 
| 140 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff | 
| 141 | > | //group therefore the total number of cutoff groups in the system is | 
| 142 | > | //equal to the total number of atoms minus number of atoms belong to | 
| 143 | > | //cutoff group defined in meta-data file plus the number of cutoff | 
| 144 | > | //groups defined in meta-data file | 
| 145 |  | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | 
| 146 |  |  | 
| 147 | < | //every free atom (atom does not belong to rigid bodies) is an integrable object | 
| 148 | < | //therefore the total number of  integrable objects in the system is equal to | 
| 149 | < | //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data | 
| 150 | < | //file plus the number of  rigid bodies defined in meta-data file | 
| 151 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; | 
| 152 | < |  | 
| 147 | > | //every free atom (atom does not belong to rigid bodies) is an | 
| 148 | > | //integrable object therefore the total number of integrable objects | 
| 149 | > | //in the system is equal to the total number of atoms minus number of | 
| 150 | > | //atoms belong to rigid body defined in meta-data file plus the number | 
| 151 | > | //of rigid bodies defined in meta-data file | 
| 152 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms | 
| 153 | > | + nGlobalRigidBodies_; | 
| 154 | > |  | 
| 155 |  | nGlobalMols_ = molStampIds_.size(); | 
| 156 |  |  | 
| 157 |  | #ifdef IS_MPI | 
| 167 |  | } | 
| 168 |  | molecules_.clear(); | 
| 169 |  |  | 
| 152 | – | delete stamps_; | 
| 170 |  | delete sman_; | 
| 171 |  | delete simParams_; | 
| 172 |  | delete forceField_; | 
| 273 |  | } | 
| 274 |  | } | 
| 275 |  |  | 
| 276 | < | }//end for (integrableObject) | 
| 277 | < | }// end for (mol) | 
| 276 | > | } | 
| 277 | > | } | 
| 278 |  |  | 
| 279 |  | // n_constraints is local, so subtract them on each processor | 
| 280 |  | ndf_local -= nConstraints_; | 
| 291 |  |  | 
| 292 |  | } | 
| 293 |  |  | 
| 294 | + | int SimInfo::getFdf() { | 
| 295 | + | #ifdef IS_MPI | 
| 296 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 297 | + | #else | 
| 298 | + | fdf_ = fdf_local; | 
| 299 | + | #endif | 
| 300 | + | return fdf_; | 
| 301 | + | } | 
| 302 | + |  | 
| 303 |  | void SimInfo::calcNdfRaw() { | 
| 304 |  | int ndfRaw_local; | 
| 305 |  |  | 
| 362 |  | int b; | 
| 363 |  | int c; | 
| 364 |  | int d; | 
| 365 | + |  | 
| 366 | + | std::map<int, std::set<int> > atomGroups; | 
| 367 | + |  | 
| 368 | + | Molecule::RigidBodyIterator rbIter; | 
| 369 | + | RigidBody* rb; | 
| 370 | + | Molecule::IntegrableObjectIterator ii; | 
| 371 | + | StuntDouble* integrableObject; | 
| 372 |  |  | 
| 373 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 374 | + | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 375 | + |  | 
| 376 | + | if (integrableObject->isRigidBody()) { | 
| 377 | + | rb = static_cast<RigidBody*>(integrableObject); | 
| 378 | + | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 379 | + | std::set<int> rigidAtoms; | 
| 380 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 381 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 382 | + | } | 
| 383 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 384 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 385 | + | } | 
| 386 | + | } else { | 
| 387 | + | std::set<int> oneAtomSet; | 
| 388 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 389 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 390 | + | } | 
| 391 | + | } | 
| 392 | + |  | 
| 393 | + |  | 
| 394 | + |  | 
| 395 |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 396 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 397 |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 402 |  | a = bend->getAtomA()->getGlobalIndex(); | 
| 403 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 404 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 405 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 406 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 407 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 408 |  |  | 
| 409 | < | exclude_.addPair(a, b); | 
| 410 | < | exclude_.addPair(a, c); | 
| 411 | < | exclude_.addPair(b, c); | 
| 409 | > | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 410 | > | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 411 | > | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 412 | > |  | 
| 413 | > | //exclude_.addPair(a, b); | 
| 414 | > | //exclude_.addPair(a, c); | 
| 415 | > | //exclude_.addPair(b, c); | 
| 416 |  | } | 
| 417 |  |  | 
| 418 |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 420 |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 421 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 422 |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 423 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 424 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 425 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 426 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 427 |  |  | 
| 428 | + | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 429 | + | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 430 | + | exclude_.addPairs(rigidSetA, rigidSetD); | 
| 431 | + | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 432 | + | exclude_.addPairs(rigidSetB, rigidSetD); | 
| 433 | + | exclude_.addPairs(rigidSetC, rigidSetD); | 
| 434 | + |  | 
| 435 | + | /* | 
| 436 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 437 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 438 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 439 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 440 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 441 | + | exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 442 | + |  | 
| 443 | + |  | 
| 444 |  | exclude_.addPair(a, b); | 
| 445 |  | exclude_.addPair(a, c); | 
| 446 |  | exclude_.addPair(a, d); | 
| 447 |  | exclude_.addPair(b, c); | 
| 448 |  | exclude_.addPair(b, d); | 
| 449 |  | exclude_.addPair(c, d); | 
| 450 | + | */ | 
| 451 |  | } | 
| 452 |  |  | 
| 370 | – | Molecule::RigidBodyIterator rbIter; | 
| 371 | – | RigidBody* rb; | 
| 453 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 454 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 455 |  | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 474 |  | int b; | 
| 475 |  | int c; | 
| 476 |  | int d; | 
| 477 | + |  | 
| 478 | + | std::map<int, std::set<int> > atomGroups; | 
| 479 | + |  | 
| 480 | + | Molecule::RigidBodyIterator rbIter; | 
| 481 | + | RigidBody* rb; | 
| 482 | + | Molecule::IntegrableObjectIterator ii; | 
| 483 | + | StuntDouble* integrableObject; | 
| 484 |  |  | 
| 485 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 486 | + | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 487 | + |  | 
| 488 | + | if (integrableObject->isRigidBody()) { | 
| 489 | + | rb = static_cast<RigidBody*>(integrableObject); | 
| 490 | + | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 491 | + | std::set<int> rigidAtoms; | 
| 492 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 493 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 494 | + | } | 
| 495 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 496 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 497 | + | } | 
| 498 | + | } else { | 
| 499 | + | std::set<int> oneAtomSet; | 
| 500 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 501 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 502 | + | } | 
| 503 | + | } | 
| 504 | + |  | 
| 505 | + |  | 
| 506 |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 507 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 508 |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 514 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 515 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 516 |  |  | 
| 517 | < | exclude_.removePair(a, b); | 
| 518 | < | exclude_.removePair(a, c); | 
| 519 | < | exclude_.removePair(b, c); | 
| 517 | > | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 518 | > | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 519 | > | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 520 | > |  | 
| 521 | > | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 522 | > | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 523 | > | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 524 | > |  | 
| 525 | > | //exclude_.removePair(a, b); | 
| 526 | > | //exclude_.removePair(a, c); | 
| 527 | > | //exclude_.removePair(b, c); | 
| 528 |  | } | 
| 529 |  |  | 
| 530 |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 533 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 534 |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 535 |  |  | 
| 536 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 537 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 538 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 539 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 540 | + |  | 
| 541 | + | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 542 | + | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 543 | + | exclude_.removePairs(rigidSetA, rigidSetD); | 
| 544 | + | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 545 | + | exclude_.removePairs(rigidSetB, rigidSetD); | 
| 546 | + | exclude_.removePairs(rigidSetC, rigidSetD); | 
| 547 | + |  | 
| 548 | + | /* | 
| 549 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 550 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 551 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 552 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 553 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 554 | + | exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 555 | + |  | 
| 556 | + |  | 
| 557 |  | exclude_.removePair(a, b); | 
| 558 |  | exclude_.removePair(a, c); | 
| 559 |  | exclude_.removePair(a, d); | 
| 560 |  | exclude_.removePair(b, c); | 
| 561 |  | exclude_.removePair(b, d); | 
| 562 |  | exclude_.removePair(c, d); | 
| 563 | + | */ | 
| 564 |  | } | 
| 565 |  |  | 
| 427 | – | Molecule::RigidBodyIterator rbIter; | 
| 428 | – | RigidBody* rb; | 
| 566 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 567 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 568 |  | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 600 |  | //setup fortran force field | 
| 601 |  | /** @deprecate */ | 
| 602 |  | int isError = 0; | 
| 603 | < | initFortranFF( &fInfo_.SIM_uses_RF, &fInfo_.SIM_uses_UW, | 
| 604 | < | &fInfo_.SIM_uses_DW, &isError ); | 
| 603 | > |  | 
| 604 | > | setupCutoff(); | 
| 605 | > |  | 
| 606 | > | setupElectrostaticSummationMethod( isError ); | 
| 607 | > | setupSwitchingFunction(); | 
| 608 | > | setupAccumulateBoxDipole(); | 
| 609 | > |  | 
| 610 |  | if(isError){ | 
| 611 |  | sprintf( painCave.errMsg, | 
| 612 |  | "ForceField error: There was an error initializing the forceField in fortran.\n" ); | 
| 613 |  | painCave.isFatal = 1; | 
| 614 |  | simError(); | 
| 615 |  | } | 
| 474 | – |  | 
| 475 | – |  | 
| 476 | – | setupCutoff(); | 
| 616 |  |  | 
| 617 |  | calcNdf(); | 
| 618 |  | calcNdfRaw(); | 
| 647 |  | int useLennardJones = 0; | 
| 648 |  | int useElectrostatic = 0; | 
| 649 |  | int useEAM = 0; | 
| 650 | + | int useSC = 0; | 
| 651 |  | int useCharge = 0; | 
| 652 |  | int useDirectional = 0; | 
| 653 |  | int useDipole = 0; | 
| 659 |  | int useDirectionalAtom = 0; | 
| 660 |  | int useElectrostatics = 0; | 
| 661 |  | //usePBC and useRF are from simParams | 
| 662 | < | int usePBC = simParams_->getPBC(); | 
| 663 | < | int useRF = simParams_->getUseRF(); | 
| 664 | < | int useUW = simParams_->getUseUndampedWolf(); | 
| 665 | < | int useDW = simParams_->getUseDampedWolf(); | 
| 662 | > | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | 
| 663 | > | int useRF; | 
| 664 | > | int useSF; | 
| 665 | > | int useSP; | 
| 666 | > | int useBoxDipole; | 
| 667 | > | std::string myMethod; | 
| 668 |  |  | 
| 669 | + | // set the useRF logical | 
| 670 | + | useRF = 0; | 
| 671 | + | useSF = 0; | 
| 672 | + | useSP = 0; | 
| 673 | + |  | 
| 674 | + |  | 
| 675 | + | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 676 | + | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 677 | + | toUpper(myMethod); | 
| 678 | + | if (myMethod == "REACTION_FIELD"){ | 
| 679 | + | useRF = 1; | 
| 680 | + | } else if (myMethod == "SHIFTED_FORCE"){ | 
| 681 | + | useSF = 1; | 
| 682 | + | } else if (myMethod == "SHIFTED_POTENTIAL"){ | 
| 683 | + | useSP = 1; | 
| 684 | + | } | 
| 685 | + | } | 
| 686 | + |  | 
| 687 | + | if (simParams_->haveAccumulateBoxDipole()) | 
| 688 | + | if (simParams_->getAccumulateBoxDipole()) | 
| 689 | + | useBoxDipole = 1; | 
| 690 | + |  | 
| 691 |  | //loop over all of the atom types | 
| 692 |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 693 |  | useLennardJones |= (*i)->isLennardJones(); | 
| 694 |  | useElectrostatic |= (*i)->isElectrostatic(); | 
| 695 |  | useEAM |= (*i)->isEAM(); | 
| 696 | + | useSC |= (*i)->isSC(); | 
| 697 |  | useCharge |= (*i)->isCharge(); | 
| 698 |  | useDirectional |= (*i)->isDirectional(); | 
| 699 |  | useDipole |= (*i)->isDipole(); | 
| 744 |  | temp = useEAM; | 
| 745 |  | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 746 |  |  | 
| 747 | + | temp = useSC; | 
| 748 | + | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 749 | + |  | 
| 750 |  | temp = useShape; | 
| 751 |  | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 752 |  |  | 
| 756 |  | temp = useRF; | 
| 757 |  | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 758 |  |  | 
| 759 | < | temp = useUW; | 
| 760 | < | MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 759 | > | temp = useSF; | 
| 760 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 761 |  |  | 
| 762 | < | temp = useDW; | 
| 763 | < | MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 764 | < |  | 
| 762 | > | temp = useSP; | 
| 763 | > | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 764 | > |  | 
| 765 | > | temp = useBoxDipole; | 
| 766 | > | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 767 | > |  | 
| 768 |  | #endif | 
| 769 |  |  | 
| 770 |  | fInfo_.SIM_uses_PBC = usePBC; | 
| 777 |  | fInfo_.SIM_uses_StickyPower = useStickyPower; | 
| 778 |  | fInfo_.SIM_uses_GayBerne = useGayBerne; | 
| 779 |  | fInfo_.SIM_uses_EAM = useEAM; | 
| 780 | + | fInfo_.SIM_uses_SC = useSC; | 
| 781 |  | fInfo_.SIM_uses_Shapes = useShape; | 
| 782 |  | fInfo_.SIM_uses_FLARB = useFLARB; | 
| 783 |  | fInfo_.SIM_uses_RF = useRF; | 
| 784 | < | fInfo_.SIM_uses_UW = useUW; | 
| 785 | < | fInfo_.SIM_uses_DW = useDW; | 
| 786 | < |  | 
| 615 | < | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { | 
| 616 | < |  | 
| 617 | < | if (simParams_->haveDielectric()) { | 
| 618 | < | fInfo_.dielect = simParams_->getDielectric(); | 
| 619 | < | } else { | 
| 620 | < | sprintf(painCave.errMsg, | 
| 621 | < | "SimSetup Error: No Dielectric constant was set.\n" | 
| 622 | < | "\tYou are trying to use Reaction Field without" | 
| 623 | < | "\tsetting a dielectric constant!\n"); | 
| 624 | < | painCave.isFatal = 1; | 
| 625 | < | simError(); | 
| 626 | < | } | 
| 627 | < |  | 
| 628 | < | } else { | 
| 629 | < | fInfo_.dielect = 0.0; | 
| 630 | < | } | 
| 631 | < |  | 
| 784 | > | fInfo_.SIM_uses_SF = useSF; | 
| 785 | > | fInfo_.SIM_uses_SP = useSP; | 
| 786 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; | 
| 787 |  | } | 
| 788 |  |  | 
| 789 |  | void SimInfo::setupFortranSim() { | 
| 800 |  | } | 
| 801 |  |  | 
| 802 |  | //calculate mass ratio of cutoff group | 
| 803 | < | std::vector<double> mfact; | 
| 803 | > | std::vector<RealType> mfact; | 
| 804 |  | SimInfo::MoleculeIterator mi; | 
| 805 |  | Molecule* mol; | 
| 806 |  | Molecule::CutoffGroupIterator ci; | 
| 807 |  | CutoffGroup* cg; | 
| 808 |  | Molecule::AtomIterator ai; | 
| 809 |  | Atom* atom; | 
| 810 | < | double totalMass; | 
| 810 | > | RealType totalMass; | 
| 811 |  |  | 
| 812 |  | //to avoid memory reallocation, reserve enough space for mfact | 
| 813 |  | mfact.reserve(getNCutoffGroups()); | 
| 817 |  |  | 
| 818 |  | totalMass = cg->getMass(); | 
| 819 |  | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | 
| 820 | < | mfact.push_back(atom->getMass()/totalMass); | 
| 820 | > | // Check for massless groups - set mfact to 1 if true | 
| 821 | > | if (totalMass != 0) | 
| 822 | > | mfact.push_back(atom->getMass()/totalMass); | 
| 823 | > | else | 
| 824 | > | mfact.push_back( 1.0 ); | 
| 825 |  | } | 
| 826 |  |  | 
| 827 |  | } | 
| 930 |  |  | 
| 931 |  | #endif | 
| 932 |  |  | 
| 933 | < | double SimInfo::calcMaxCutoffRadius() { | 
| 933 | > | void SimInfo::setupCutoff() { | 
| 934 | > |  | 
| 935 | > | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | 
| 936 |  |  | 
| 937 | + | // Check the cutoff policy | 
| 938 | + | int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | 
| 939 |  |  | 
| 940 | < | std::set<AtomType*> atomTypes; | 
| 941 | < | std::set<AtomType*>::iterator i; | 
| 942 | < | std::vector<double> cutoffRadius; | 
| 943 | < |  | 
| 944 | < | //get the unique atom types | 
| 782 | < | atomTypes = getUniqueAtomTypes(); | 
| 783 | < |  | 
| 784 | < | //query the max cutoff radius among these atom types | 
| 785 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 786 | < | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); | 
| 940 | > | std::string myPolicy; | 
| 941 | > | if (forceFieldOptions_.haveCutoffPolicy()){ | 
| 942 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); | 
| 943 | > | }else if (simParams_->haveCutoffPolicy()) { | 
| 944 | > | myPolicy = simParams_->getCutoffPolicy(); | 
| 945 |  | } | 
| 946 |  |  | 
| 947 | < | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); | 
| 948 | < | #ifdef IS_MPI | 
| 791 | < | //pick the max cutoff radius among the processors | 
| 792 | < | #endif | 
| 793 | < |  | 
| 794 | < | return maxCutoffRadius; | 
| 795 | < | } | 
| 796 | < |  | 
| 797 | < | void SimInfo::getCutoff(double& rcut, double& rsw) { | 
| 798 | < |  | 
| 799 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 800 | < |  | 
| 801 | < | if (!simParams_->haveRcut()){ | 
| 802 | < | sprintf(painCave.errMsg, | 
| 803 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 804 | < | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 805 | < | "\tfor the cutoffRadius.\n"); | 
| 806 | < | painCave.isFatal = 0; | 
| 807 | < | simError(); | 
| 808 | < | rcut = 15.0; | 
| 809 | < | } else{ | 
| 810 | < | rcut = simParams_->getRcut(); | 
| 811 | < | } | 
| 812 | < |  | 
| 813 | < | if (!simParams_->haveRsw()){ | 
| 814 | < | sprintf(painCave.errMsg, | 
| 815 | < | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 816 | < | "\tOOPSE will use a default value of\n" | 
| 817 | < | "\t0.95 * cutoffRadius for the switchingRadius\n"); | 
| 818 | < | painCave.isFatal = 0; | 
| 819 | < | simError(); | 
| 820 | < | rsw = 0.95 * rcut; | 
| 821 | < | } else{ | 
| 822 | < | rsw = simParams_->getRsw(); | 
| 823 | < | } | 
| 824 | < |  | 
| 825 | < | } else { | 
| 826 | < | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | 
| 827 | < | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | 
| 828 | < |  | 
| 829 | < | if (simParams_->haveRcut()) { | 
| 830 | < | rcut = simParams_->getRcut(); | 
| 831 | < | } else { | 
| 832 | < | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | 
| 833 | < | rcut = calcMaxCutoffRadius(); | 
| 834 | < | } | 
| 835 | < |  | 
| 836 | < | if (simParams_->haveRsw()) { | 
| 837 | < | rsw  = simParams_->getRsw(); | 
| 838 | < | } else { | 
| 839 | < | rsw = rcut; | 
| 840 | < | } | 
| 841 | < |  | 
| 842 | < | } | 
| 843 | < | } | 
| 844 | < |  | 
| 845 | < | void SimInfo::setupCutoff() { | 
| 846 | < | getCutoff(rcut_, rsw_); | 
| 847 | < | double rnblist = rcut_ + 1; // skin of neighbor list | 
| 848 | < |  | 
| 849 | < | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | 
| 850 | < |  | 
| 851 | < | int cp =  TRADITIONAL_CUTOFF_POLICY; | 
| 852 | < | if (simParams_->haveCutoffPolicy()) { | 
| 853 | < | std::string myPolicy = simParams_->getCutoffPolicy(); | 
| 947 | > | if (!myPolicy.empty()){ | 
| 948 | > | toUpper(myPolicy); | 
| 949 |  | if (myPolicy == "MIX") { | 
| 950 |  | cp = MIX_CUTOFF_POLICY; | 
| 951 |  | } else { | 
| 963 |  | } | 
| 964 |  | } | 
| 965 |  | } | 
| 966 | + | } | 
| 967 | + | notifyFortranCutoffPolicy(&cp); | 
| 968 | + |  | 
| 969 | + | // Check the Skin Thickness for neighborlists | 
| 970 | + | RealType skin; | 
| 971 | + | if (simParams_->haveSkinThickness()) { | 
| 972 | + | skin = simParams_->getSkinThickness(); | 
| 973 | + | notifyFortranSkinThickness(&skin); | 
| 974 | + | } | 
| 975 | + |  | 
| 976 | + | // Check if the cutoff was set explicitly: | 
| 977 | + | if (simParams_->haveCutoffRadius()) { | 
| 978 | + | rcut_ = simParams_->getCutoffRadius(); | 
| 979 | + | if (simParams_->haveSwitchingRadius()) { | 
| 980 | + | rsw_  = simParams_->getSwitchingRadius(); | 
| 981 | + | } else { | 
| 982 | + | if (fInfo_.SIM_uses_Charges | | 
| 983 | + | fInfo_.SIM_uses_Dipoles | | 
| 984 | + | fInfo_.SIM_uses_RF) { | 
| 985 | + |  | 
| 986 | + | rsw_ = 0.85 * rcut_; | 
| 987 | + | sprintf(painCave.errMsg, | 
| 988 | + | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 989 | + | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" | 
| 990 | + | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 991 | + | painCave.isFatal = 0; | 
| 992 | + | simError(); | 
| 993 | + | } else { | 
| 994 | + | rsw_ = rcut_; | 
| 995 | + | sprintf(painCave.errMsg, | 
| 996 | + | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 997 | + | "\tOOPSE will use the same value as the cutoffRadius.\n" | 
| 998 | + | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 999 | + | painCave.isFatal = 0; | 
| 1000 | + | simError(); | 
| 1001 | + | } | 
| 1002 | + | } | 
| 1003 | + |  | 
| 1004 | + | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1005 | + |  | 
| 1006 | + | } else { | 
| 1007 | + |  | 
| 1008 | + | // For electrostatic atoms, we'll assume a large safe value: | 
| 1009 | + | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 1010 | + | sprintf(painCave.errMsg, | 
| 1011 | + | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 1012 | + | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 1013 | + | "\tfor the cutoffRadius.\n"); | 
| 1014 | + | painCave.isFatal = 0; | 
| 1015 | + | simError(); | 
| 1016 | + | rcut_ = 15.0; | 
| 1017 | + |  | 
| 1018 | + | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1019 | + | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1020 | + | toUpper(myMethod); | 
| 1021 | + | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 1022 | + | if (simParams_->haveSwitchingRadius()){ | 
| 1023 | + | sprintf(painCave.errMsg, | 
| 1024 | + | "SimInfo Warning: A value was set for the switchingRadius\n" | 
| 1025 | + | "\teven though the electrostaticSummationMethod was\n" | 
| 1026 | + | "\tset to %s\n", myMethod.c_str()); | 
| 1027 | + | painCave.isFatal = 1; | 
| 1028 | + | simError(); | 
| 1029 | + | } | 
| 1030 | + | } | 
| 1031 | + | } | 
| 1032 | + |  | 
| 1033 | + | if (simParams_->haveSwitchingRadius()){ | 
| 1034 | + | rsw_ = simParams_->getSwitchingRadius(); | 
| 1035 | + | } else { | 
| 1036 | + | sprintf(painCave.errMsg, | 
| 1037 | + | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 1038 | + | "\tOOPSE will use a default value of\n" | 
| 1039 | + | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 1040 | + | painCave.isFatal = 0; | 
| 1041 | + | simError(); | 
| 1042 | + | rsw_ = 0.85 * rcut_; | 
| 1043 | + | } | 
| 1044 | + | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1045 | + | } else { | 
| 1046 | + | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | 
| 1047 | + | // We'll punt and let fortran figure out the cutoffs later. | 
| 1048 | + |  | 
| 1049 | + | notifyFortranYouAreOnYourOwn(); | 
| 1050 | + |  | 
| 1051 | + | } | 
| 1052 |  | } | 
| 1053 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | 
| 1053 | > | } | 
| 1054 | > |  | 
| 1055 | > | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | 
| 1056 | > |  | 
| 1057 | > | int errorOut; | 
| 1058 | > | int esm =  NONE; | 
| 1059 | > | int sm = UNDAMPED; | 
| 1060 | > | RealType alphaVal; | 
| 1061 | > | RealType dielectric; | 
| 1062 | > |  | 
| 1063 | > | errorOut = isError; | 
| 1064 | > |  | 
| 1065 | > | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1066 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1067 | > | toUpper(myMethod); | 
| 1068 | > | if (myMethod == "NONE") { | 
| 1069 | > | esm = NONE; | 
| 1070 | > | } else { | 
| 1071 | > | if (myMethod == "SWITCHING_FUNCTION") { | 
| 1072 | > | esm = SWITCHING_FUNCTION; | 
| 1073 | > | } else { | 
| 1074 | > | if (myMethod == "SHIFTED_POTENTIAL") { | 
| 1075 | > | esm = SHIFTED_POTENTIAL; | 
| 1076 | > | } else { | 
| 1077 | > | if (myMethod == "SHIFTED_FORCE") { | 
| 1078 | > | esm = SHIFTED_FORCE; | 
| 1079 | > | } else { | 
| 1080 | > | if (myMethod == "REACTION_FIELD") { | 
| 1081 | > | esm = REACTION_FIELD; | 
| 1082 | > | dielectric = simParams_->getDielectric(); | 
| 1083 | > | if (!simParams_->haveDielectric()) { | 
| 1084 | > | // throw warning | 
| 1085 | > | sprintf( painCave.errMsg, | 
| 1086 | > | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" | 
| 1087 | > | "\tA default value of %f will be used for the dielectric.\n", dielectric); | 
| 1088 | > | painCave.isFatal = 0; | 
| 1089 | > | simError(); | 
| 1090 | > | } | 
| 1091 | > | } else { | 
| 1092 | > | // throw error | 
| 1093 | > | sprintf( painCave.errMsg, | 
| 1094 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" | 
| 1095 | > | "\t(Input file specified %s .)\n" | 
| 1096 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" | 
| 1097 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" | 
| 1098 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); | 
| 1099 | > | painCave.isFatal = 1; | 
| 1100 | > | simError(); | 
| 1101 | > | } | 
| 1102 | > | } | 
| 1103 | > | } | 
| 1104 | > | } | 
| 1105 | > | } | 
| 1106 | > | } | 
| 1107 | > |  | 
| 1108 | > | if (simParams_->haveElectrostaticScreeningMethod()) { | 
| 1109 | > | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); | 
| 1110 | > | toUpper(myScreen); | 
| 1111 | > | if (myScreen == "UNDAMPED") { | 
| 1112 | > | sm = UNDAMPED; | 
| 1113 | > | } else { | 
| 1114 | > | if (myScreen == "DAMPED") { | 
| 1115 | > | sm = DAMPED; | 
| 1116 | > | if (!simParams_->haveDampingAlpha()) { | 
| 1117 | > | // first set a cutoff dependent alpha value | 
| 1118 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang | 
| 1119 | > | alphaVal = 0.5125 - rcut_* 0.025; | 
| 1120 | > | // for values rcut > 20.5, alpha is zero | 
| 1121 | > | if (alphaVal < 0) alphaVal = 0; | 
| 1122 | > |  | 
| 1123 | > | // throw warning | 
| 1124 | > | sprintf( painCave.errMsg, | 
| 1125 | > | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1126 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); | 
| 1127 | > | painCave.isFatal = 0; | 
| 1128 | > | simError(); | 
| 1129 | > | } | 
| 1130 | > | } else { | 
| 1131 | > | // throw error | 
| 1132 | > | sprintf( painCave.errMsg, | 
| 1133 | > | "SimInfo error: Unknown electrostaticScreeningMethod.\n" | 
| 1134 | > | "\t(Input file specified %s .)\n" | 
| 1135 | > | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" | 
| 1136 | > | "or \"damped\".\n", myScreen.c_str() ); | 
| 1137 | > | painCave.isFatal = 1; | 
| 1138 | > | simError(); | 
| 1139 | > | } | 
| 1140 | > | } | 
| 1141 | > | } | 
| 1142 | > |  | 
| 1143 | > | // let's pass some summation method variables to fortran | 
| 1144 | > | setElectrostaticSummationMethod( &esm ); | 
| 1145 | > | setFortranElectrostaticMethod( &esm ); | 
| 1146 | > | setScreeningMethod( &sm ); | 
| 1147 | > | setDampingAlpha( &alphaVal ); | 
| 1148 | > | setReactionFieldDielectric( &dielectric ); | 
| 1149 | > | initFortranFF( &errorOut ); | 
| 1150 |  | } | 
| 1151 |  |  | 
| 1152 | + | void SimInfo::setupSwitchingFunction() { | 
| 1153 | + | int ft = CUBIC; | 
| 1154 | + |  | 
| 1155 | + | if (simParams_->haveSwitchingFunctionType()) { | 
| 1156 | + | std::string funcType = simParams_->getSwitchingFunctionType(); | 
| 1157 | + | toUpper(funcType); | 
| 1158 | + | if (funcType == "CUBIC") { | 
| 1159 | + | ft = CUBIC; | 
| 1160 | + | } else { | 
| 1161 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { | 
| 1162 | + | ft = FIFTH_ORDER_POLY; | 
| 1163 | + | } else { | 
| 1164 | + | // throw error | 
| 1165 | + | sprintf( painCave.errMsg, | 
| 1166 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); | 
| 1167 | + | painCave.isFatal = 1; | 
| 1168 | + | simError(); | 
| 1169 | + | } | 
| 1170 | + | } | 
| 1171 | + | } | 
| 1172 | + |  | 
| 1173 | + | // send switching function notification to switcheroo | 
| 1174 | + | setFunctionType(&ft); | 
| 1175 | + |  | 
| 1176 | + | } | 
| 1177 | + |  | 
| 1178 | + | void SimInfo::setupAccumulateBoxDipole() { | 
| 1179 | + |  | 
| 1180 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 1181 | + | if ( simParams_->haveAccumulateBoxDipole() ) | 
| 1182 | + | if ( simParams_->getAccumulateBoxDipole() ) { | 
| 1183 | + | setAccumulateBoxDipole(); | 
| 1184 | + | calcBoxDipole_ = true; | 
| 1185 | + | } | 
| 1186 | + |  | 
| 1187 | + | } | 
| 1188 | + |  | 
| 1189 |  | void SimInfo::addProperty(GenericData* genData) { | 
| 1190 |  | properties_.addProperty(genData); | 
| 1191 |  | } | 
| 1242 |  | Molecule* mol; | 
| 1243 |  |  | 
| 1244 |  | Vector3d comVel(0.0); | 
| 1245 | < | double totalMass = 0.0; | 
| 1245 | > | RealType totalMass = 0.0; | 
| 1246 |  |  | 
| 1247 |  |  | 
| 1248 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1249 | < | double mass = mol->getMass(); | 
| 1249 | > | RealType mass = mol->getMass(); | 
| 1250 |  | totalMass += mass; | 
| 1251 |  | comVel += mass * mol->getComVel(); | 
| 1252 |  | } | 
| 1253 |  |  | 
| 1254 |  | #ifdef IS_MPI | 
| 1255 | < | double tmpMass = totalMass; | 
| 1255 | > | RealType tmpMass = totalMass; | 
| 1256 |  | Vector3d tmpComVel(comVel); | 
| 1257 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1258 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1257 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1258 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1259 |  | #endif | 
| 1260 |  |  | 
| 1261 |  | comVel /= totalMass; | 
| 1268 |  | Molecule* mol; | 
| 1269 |  |  | 
| 1270 |  | Vector3d com(0.0); | 
| 1271 | < | double totalMass = 0.0; | 
| 1271 | > | RealType totalMass = 0.0; | 
| 1272 |  |  | 
| 1273 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1274 | < | double mass = mol->getMass(); | 
| 1274 | > | RealType mass = mol->getMass(); | 
| 1275 |  | totalMass += mass; | 
| 1276 |  | com += mass * mol->getCom(); | 
| 1277 |  | } | 
| 1278 |  |  | 
| 1279 |  | #ifdef IS_MPI | 
| 1280 | < | double tmpMass = totalMass; | 
| 1280 | > | RealType tmpMass = totalMass; | 
| 1281 |  | Vector3d tmpCom(com); | 
| 1282 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1283 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1282 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1283 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1284 |  | #endif | 
| 1285 |  |  | 
| 1286 |  | com /= totalMass; | 
| 1304 |  | Molecule* mol; | 
| 1305 |  |  | 
| 1306 |  |  | 
| 1307 | < | double totalMass = 0.0; | 
| 1307 | > | RealType totalMass = 0.0; | 
| 1308 |  |  | 
| 1309 |  |  | 
| 1310 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1311 | < | double mass = mol->getMass(); | 
| 1311 | > | RealType mass = mol->getMass(); | 
| 1312 |  | totalMass += mass; | 
| 1313 |  | com += mass * mol->getCom(); | 
| 1314 |  | comVel += mass * mol->getComVel(); | 
| 1315 |  | } | 
| 1316 |  |  | 
| 1317 |  | #ifdef IS_MPI | 
| 1318 | < | double tmpMass = totalMass; | 
| 1318 | > | RealType tmpMass = totalMass; | 
| 1319 |  | Vector3d tmpCom(com); | 
| 1320 |  | Vector3d tmpComVel(comVel); | 
| 1321 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1322 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1323 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1321 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1322 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1323 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1324 |  | #endif | 
| 1325 |  |  | 
| 1326 |  | com /= totalMass; | 
| 1339 |  | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | 
| 1340 |  |  | 
| 1341 |  |  | 
| 1342 | < | double xx = 0.0; | 
| 1343 | < | double yy = 0.0; | 
| 1344 | < | double zz = 0.0; | 
| 1345 | < | double xy = 0.0; | 
| 1346 | < | double xz = 0.0; | 
| 1347 | < | double yz = 0.0; | 
| 1342 | > | RealType xx = 0.0; | 
| 1343 | > | RealType yy = 0.0; | 
| 1344 | > | RealType zz = 0.0; | 
| 1345 | > | RealType xy = 0.0; | 
| 1346 | > | RealType xz = 0.0; | 
| 1347 | > | RealType yz = 0.0; | 
| 1348 |  | Vector3d com(0.0); | 
| 1349 |  | Vector3d comVel(0.0); | 
| 1350 |  |  | 
| 1356 |  | Vector3d thisq(0.0); | 
| 1357 |  | Vector3d thisv(0.0); | 
| 1358 |  |  | 
| 1359 | < | double thisMass = 0.0; | 
| 1359 | > | RealType thisMass = 0.0; | 
| 1360 |  |  | 
| 1361 |  |  | 
| 1362 |  |  | 
| 1394 |  | #ifdef IS_MPI | 
| 1395 |  | Mat3x3d tmpI(inertiaTensor); | 
| 1396 |  | Vector3d tmpAngMom; | 
| 1397 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1398 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1397 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1398 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1399 |  | #endif | 
| 1400 |  |  | 
| 1401 |  | return; | 
| 1416 |  | Vector3d thisr(0.0); | 
| 1417 |  | Vector3d thisp(0.0); | 
| 1418 |  |  | 
| 1419 | < | double thisMass; | 
| 1419 | > | RealType thisMass; | 
| 1420 |  |  | 
| 1421 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1422 |  | thisMass = mol->getMass(); | 
| 1429 |  |  | 
| 1430 |  | #ifdef IS_MPI | 
| 1431 |  | Vector3d tmpAngMom; | 
| 1432 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1432 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1433 |  | #endif | 
| 1434 |  |  | 
| 1435 |  | return angularMomentum; | 
| 1436 |  | } | 
| 1437 |  |  | 
| 1438 | < |  | 
| 1438 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | 
| 1439 | > | return IOIndexToIntegrableObject.at(index); | 
| 1440 | > | } | 
| 1441 | > |  | 
| 1442 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { | 
| 1443 | > | IOIndexToIntegrableObject= v; | 
| 1444 | > | } | 
| 1445 | > |  | 
| 1446 | > | /* | 
| 1447 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { | 
| 1448 | > | assert( v.size() == nAtoms_ + nRigidBodies_); | 
| 1449 | > | sdByGlobalIndex_ = v; | 
| 1450 | > | } | 
| 1451 | > |  | 
| 1452 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { | 
| 1453 | > | //assert(index < nAtoms_ + nRigidBodies_); | 
| 1454 | > | return sdByGlobalIndex_.at(index); | 
| 1455 | > | } | 
| 1456 | > | */ | 
| 1457 |  | }//end namespace oopse | 
| 1458 |  |  |