| 48 |  |  | 
| 49 |  | #include <algorithm> | 
| 50 |  | #include <set> | 
| 51 | + | #include <map> | 
| 52 |  |  | 
| 53 |  | #include "brains/SimInfo.hpp" | 
| 54 |  | #include "math/Vector3.hpp" | 
| 59 |  | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" | 
| 60 |  | #include "UseTheForce/doForces_interface.h" | 
| 61 |  | #include "UseTheForce/DarkSide/electrostatic_interface.h" | 
| 61 | – | #include "UseTheForce/notifyCutoffs_interface.h" | 
| 62 |  | #include "UseTheForce/DarkSide/switcheroo_interface.h" | 
| 63 |  | #include "utils/MemoryUtils.hpp" | 
| 64 |  | #include "utils/simError.h" | 
| 70 |  | #endif | 
| 71 |  |  | 
| 72 |  | namespace oopse { | 
| 73 | + | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { | 
| 74 | + | std::map<int, std::set<int> >::iterator i = container.find(index); | 
| 75 | + | std::set<int> result; | 
| 76 | + | if (i != container.end()) { | 
| 77 | + | result = i->second; | 
| 78 | + | } | 
| 79 |  |  | 
| 80 | < | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, | 
| 81 | < | ForceField* ff, Globals* simParams) : | 
| 82 | < | stamps_(stamps), forceField_(ff), simParams_(simParams), | 
| 80 | > | return result; | 
| 81 | > | } | 
| 82 | > |  | 
| 83 | > | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | 
| 84 | > | forceField_(ff), simParams_(simParams), | 
| 85 |  | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 86 |  | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | 
| 87 |  | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 89 |  | nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0), | 
| 90 |  | sman_(NULL), fortranInitialized_(false) { | 
| 91 |  |  | 
| 84 | – |  | 
| 85 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; | 
| 92 |  | MoleculeStamp* molStamp; | 
| 93 |  | int nMolWithSameStamp; | 
| 94 |  | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 96 |  | CutoffGroupStamp* cgStamp; | 
| 97 |  | RigidBodyStamp* rbStamp; | 
| 98 |  | int nRigidAtoms = 0; | 
| 99 | < |  | 
| 100 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { | 
| 101 | < | molStamp = i->first; | 
| 102 | < | nMolWithSameStamp = i->second; | 
| 99 | > | std::vector<Component*> components = simParams->getComponents(); | 
| 100 | > |  | 
| 101 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 102 | > | molStamp = (*i)->getMoleculeStamp(); | 
| 103 | > | nMolWithSameStamp = (*i)->getNMol(); | 
| 104 |  |  | 
| 105 |  | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 106 |  |  | 
| 107 |  | //calculate atoms in molecules | 
| 108 |  | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 109 |  |  | 
| 103 | – |  | 
| 110 |  | //calculate atoms in cutoff groups | 
| 111 |  | int nAtomsInGroups = 0; | 
| 112 |  | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 113 |  |  | 
| 114 |  | for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 115 | < | cgStamp = molStamp->getCutoffGroup(j); | 
| 115 | > | cgStamp = molStamp->getCutoffGroupStamp(j); | 
| 116 |  | nAtomsInGroups += cgStamp->getNMembers(); | 
| 117 |  | } | 
| 118 |  |  | 
| 125 |  | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 126 |  |  | 
| 127 |  | for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 128 | < | rbStamp = molStamp->getRigidBody(j); | 
| 128 | > | rbStamp = molStamp->getRigidBodyStamp(j); | 
| 129 |  | nAtomsInRigidBodies += rbStamp->getNMembers(); | 
| 130 |  | } | 
| 131 |  |  | 
| 164 |  | } | 
| 165 |  | molecules_.clear(); | 
| 166 |  |  | 
| 161 | – | delete stamps_; | 
| 167 |  | delete sman_; | 
| 168 |  | delete simParams_; | 
| 169 |  | delete forceField_; | 
| 270 |  | } | 
| 271 |  | } | 
| 272 |  |  | 
| 273 | < | }//end for (integrableObject) | 
| 274 | < | }// end for (mol) | 
| 273 | > | } | 
| 274 | > | } | 
| 275 |  |  | 
| 276 |  | // n_constraints is local, so subtract them on each processor | 
| 277 |  | ndf_local -= nConstraints_; | 
| 350 |  | int b; | 
| 351 |  | int c; | 
| 352 |  | int d; | 
| 353 | + |  | 
| 354 | + | std::map<int, std::set<int> > atomGroups; | 
| 355 | + |  | 
| 356 | + | Molecule::RigidBodyIterator rbIter; | 
| 357 | + | RigidBody* rb; | 
| 358 | + | Molecule::IntegrableObjectIterator ii; | 
| 359 | + | StuntDouble* integrableObject; | 
| 360 |  |  | 
| 361 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 362 | + | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 363 | + |  | 
| 364 | + | if (integrableObject->isRigidBody()) { | 
| 365 | + | rb = static_cast<RigidBody*>(integrableObject); | 
| 366 | + | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 367 | + | std::set<int> rigidAtoms; | 
| 368 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 369 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 370 | + | } | 
| 371 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 372 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 373 | + | } | 
| 374 | + | } else { | 
| 375 | + | std::set<int> oneAtomSet; | 
| 376 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 377 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 378 | + | } | 
| 379 | + | } | 
| 380 | + |  | 
| 381 | + |  | 
| 382 | + |  | 
| 383 |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 384 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 385 |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 390 |  | a = bend->getAtomA()->getGlobalIndex(); | 
| 391 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 392 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 393 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 394 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 395 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 396 |  |  | 
| 397 | < | exclude_.addPair(a, b); | 
| 398 | < | exclude_.addPair(a, c); | 
| 399 | < | exclude_.addPair(b, c); | 
| 397 | > | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 398 | > | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 399 | > | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 400 | > |  | 
| 401 | > | //exclude_.addPair(a, b); | 
| 402 | > | //exclude_.addPair(a, c); | 
| 403 | > | //exclude_.addPair(b, c); | 
| 404 |  | } | 
| 405 |  |  | 
| 406 |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 408 |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 409 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 410 |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 411 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 412 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 413 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 414 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 415 |  |  | 
| 416 | + | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 417 | + | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 418 | + | exclude_.addPairs(rigidSetA, rigidSetD); | 
| 419 | + | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 420 | + | exclude_.addPairs(rigidSetB, rigidSetD); | 
| 421 | + | exclude_.addPairs(rigidSetC, rigidSetD); | 
| 422 | + |  | 
| 423 | + | /* | 
| 424 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 425 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 426 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 427 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 428 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 429 | + | exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 430 | + |  | 
| 431 | + |  | 
| 432 |  | exclude_.addPair(a, b); | 
| 433 |  | exclude_.addPair(a, c); | 
| 434 |  | exclude_.addPair(a, d); | 
| 435 |  | exclude_.addPair(b, c); | 
| 436 |  | exclude_.addPair(b, d); | 
| 437 |  | exclude_.addPair(c, d); | 
| 438 | + | */ | 
| 439 |  | } | 
| 440 |  |  | 
| 379 | – | Molecule::RigidBodyIterator rbIter; | 
| 380 | – | RigidBody* rb; | 
| 441 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 442 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 443 |  | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 462 |  | int b; | 
| 463 |  | int c; | 
| 464 |  | int d; | 
| 465 | + |  | 
| 466 | + | std::map<int, std::set<int> > atomGroups; | 
| 467 | + |  | 
| 468 | + | Molecule::RigidBodyIterator rbIter; | 
| 469 | + | RigidBody* rb; | 
| 470 | + | Molecule::IntegrableObjectIterator ii; | 
| 471 | + | StuntDouble* integrableObject; | 
| 472 |  |  | 
| 473 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 474 | + | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 475 | + |  | 
| 476 | + | if (integrableObject->isRigidBody()) { | 
| 477 | + | rb = static_cast<RigidBody*>(integrableObject); | 
| 478 | + | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 479 | + | std::set<int> rigidAtoms; | 
| 480 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 481 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 482 | + | } | 
| 483 | + | for (int i = 0; i < atoms.size(); ++i) { | 
| 484 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 485 | + | } | 
| 486 | + | } else { | 
| 487 | + | std::set<int> oneAtomSet; | 
| 488 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 489 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 490 | + | } | 
| 491 | + | } | 
| 492 | + |  | 
| 493 | + |  | 
| 494 |  | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 495 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 496 |  | b = bond->getAtomB()->getGlobalIndex(); | 
| 502 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 503 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 504 |  |  | 
| 505 | < | exclude_.removePair(a, b); | 
| 506 | < | exclude_.removePair(a, c); | 
| 507 | < | exclude_.removePair(b, c); | 
| 505 | > | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 506 | > | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 507 | > | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 508 | > |  | 
| 509 | > | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 510 | > | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 511 | > | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 512 | > |  | 
| 513 | > | //exclude_.removePair(a, b); | 
| 514 | > | //exclude_.removePair(a, c); | 
| 515 | > | //exclude_.removePair(b, c); | 
| 516 |  | } | 
| 517 |  |  | 
| 518 |  | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 520 |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 521 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 522 |  | d = torsion->getAtomD()->getGlobalIndex(); | 
| 523 | + |  | 
| 524 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 525 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 526 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 527 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 528 | + |  | 
| 529 | + | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 530 | + | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 531 | + | exclude_.removePairs(rigidSetA, rigidSetD); | 
| 532 | + | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 533 | + | exclude_.removePairs(rigidSetB, rigidSetD); | 
| 534 | + | exclude_.removePairs(rigidSetC, rigidSetD); | 
| 535 | + |  | 
| 536 | + | /* | 
| 537 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 538 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 539 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 540 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 541 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 542 | + | exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 543 |  |  | 
| 544 | + |  | 
| 545 |  | exclude_.removePair(a, b); | 
| 546 |  | exclude_.removePair(a, c); | 
| 547 |  | exclude_.removePair(a, d); | 
| 548 |  | exclude_.removePair(b, c); | 
| 549 |  | exclude_.removePair(b, d); | 
| 550 |  | exclude_.removePair(c, d); | 
| 551 | + | */ | 
| 552 |  | } | 
| 553 |  |  | 
| 436 | – | Molecule::RigidBodyIterator rbIter; | 
| 437 | – | RigidBody* rb; | 
| 554 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 555 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 556 |  | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 915 |  |  | 
| 916 |  |  | 
| 917 |  | } | 
| 802 | – |  | 
| 803 | – | #endif | 
| 804 | – |  | 
| 805 | – | double SimInfo::calcMaxCutoffRadius() { | 
| 918 |  |  | 
| 807 | – |  | 
| 808 | – | std::set<AtomType*> atomTypes; | 
| 809 | – | std::set<AtomType*>::iterator i; | 
| 810 | – | std::vector<double> cutoffRadius; | 
| 811 | – |  | 
| 812 | – | //get the unique atom types | 
| 813 | – | atomTypes = getUniqueAtomTypes(); | 
| 814 | – |  | 
| 815 | – | //query the max cutoff radius among these atom types | 
| 816 | – | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 817 | – | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); | 
| 818 | – | } | 
| 819 | – |  | 
| 820 | – | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); | 
| 821 | – | #ifdef IS_MPI | 
| 822 | – | //pick the max cutoff radius among the processors | 
| 919 |  | #endif | 
| 920 |  |  | 
| 921 | < | return maxCutoffRadius; | 
| 826 | < | } | 
| 827 | < |  | 
| 828 | < | void SimInfo::getCutoff(double& rcut, double& rsw) { | 
| 921 | > | void SimInfo::setupCutoff() { | 
| 922 |  |  | 
| 923 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 831 | < |  | 
| 832 | < | if (!simParams_->haveCutoffRadius()){ | 
| 833 | < | sprintf(painCave.errMsg, | 
| 834 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 835 | < | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 836 | < | "\tfor the cutoffRadius.\n"); | 
| 837 | < | painCave.isFatal = 0; | 
| 838 | < | simError(); | 
| 839 | < | rcut = 15.0; | 
| 840 | < | } else{ | 
| 841 | < | rcut = simParams_->getCutoffRadius(); | 
| 842 | < | } | 
| 843 | < |  | 
| 844 | < | if (!simParams_->haveSwitchingRadius()){ | 
| 845 | < | sprintf(painCave.errMsg, | 
| 846 | < | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 847 | < | "\tOOPSE will use a default value of\n" | 
| 848 | < | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 849 | < | painCave.isFatal = 0; | 
| 850 | < | simError(); | 
| 851 | < | rsw = 0.85 * rcut; | 
| 852 | < | } else{ | 
| 853 | < | rsw = simParams_->getSwitchingRadius(); | 
| 854 | < | } | 
| 855 | < |  | 
| 856 | < | } else { | 
| 857 | < | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | 
| 858 | < | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | 
| 859 | < |  | 
| 860 | < | if (simParams_->haveCutoffRadius()) { | 
| 861 | < | rcut = simParams_->getCutoffRadius(); | 
| 862 | < | } else { | 
| 863 | < | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | 
| 864 | < | rcut = calcMaxCutoffRadius(); | 
| 865 | < | } | 
| 866 | < |  | 
| 867 | < | if (simParams_->haveSwitchingRadius()) { | 
| 868 | < | rsw  = simParams_->getSwitchingRadius(); | 
| 869 | < | } else { | 
| 870 | < | rsw = rcut; | 
| 871 | < | } | 
| 872 | < |  | 
| 873 | < | } | 
| 874 | < | } | 
| 875 | < |  | 
| 876 | < | void SimInfo::setupCutoff() { | 
| 877 | < | getCutoff(rcut_, rsw_); | 
| 878 | < | double rnblist = rcut_ + 1; // skin of neighbor list | 
| 879 | < |  | 
| 880 | < | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | 
| 881 | < |  | 
| 923 | > | // Check the cutoff policy | 
| 924 |  | int cp =  TRADITIONAL_CUTOFF_POLICY; | 
| 925 |  | if (simParams_->haveCutoffPolicy()) { | 
| 926 |  | std::string myPolicy = simParams_->getCutoffPolicy(); | 
| 942 |  | } | 
| 943 |  | } | 
| 944 |  | } | 
| 945 | < | } | 
| 945 | > | } | 
| 946 | > | notifyFortranCutoffPolicy(&cp); | 
| 947 |  |  | 
| 948 | < |  | 
| 948 | > | // Check the Skin Thickness for neighborlists | 
| 949 | > | double skin; | 
| 950 |  | if (simParams_->haveSkinThickness()) { | 
| 951 | < | double skinThickness = simParams_->getSkinThickness(); | 
| 952 | < | } | 
| 951 | > | skin = simParams_->getSkinThickness(); | 
| 952 | > | notifyFortranSkinThickness(&skin); | 
| 953 | > | } | 
| 954 | > |  | 
| 955 | > | // Check if the cutoff was set explicitly: | 
| 956 | > | if (simParams_->haveCutoffRadius()) { | 
| 957 | > | rcut_ = simParams_->getCutoffRadius(); | 
| 958 | > | if (simParams_->haveSwitchingRadius()) { | 
| 959 | > | rsw_  = simParams_->getSwitchingRadius(); | 
| 960 | > | } else { | 
| 961 | > | rsw_ = rcut_; | 
| 962 | > | } | 
| 963 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 964 | > |  | 
| 965 | > | } else { | 
| 966 | > |  | 
| 967 | > | // For electrostatic atoms, we'll assume a large safe value: | 
| 968 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 969 | > | sprintf(painCave.errMsg, | 
| 970 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 971 | > | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 972 | > | "\tfor the cutoffRadius.\n"); | 
| 973 | > | painCave.isFatal = 0; | 
| 974 | > | simError(); | 
| 975 | > | rcut_ = 15.0; | 
| 976 | > |  | 
| 977 | > | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 978 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 979 | > | toUpper(myMethod); | 
| 980 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 981 | > | if (simParams_->haveSwitchingRadius()){ | 
| 982 | > | sprintf(painCave.errMsg, | 
| 983 | > | "SimInfo Warning: A value was set for the switchingRadius\n" | 
| 984 | > | "\teven though the electrostaticSummationMethod was\n" | 
| 985 | > | "\tset to %s\n", myMethod.c_str()); | 
| 986 | > | painCave.isFatal = 1; | 
| 987 | > | simError(); | 
| 988 | > | } | 
| 989 | > | } | 
| 990 | > | } | 
| 991 | > |  | 
| 992 | > | if (simParams_->haveSwitchingRadius()){ | 
| 993 | > | rsw_ = simParams_->getSwitchingRadius(); | 
| 994 | > | } else { | 
| 995 | > | sprintf(painCave.errMsg, | 
| 996 | > | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 997 | > | "\tOOPSE will use a default value of\n" | 
| 998 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 999 | > | painCave.isFatal = 0; | 
| 1000 | > | simError(); | 
| 1001 | > | rsw_ = 0.85 * rcut_; | 
| 1002 | > | } | 
| 1003 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1004 | > | } else { | 
| 1005 | > | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | 
| 1006 | > | // We'll punt and let fortran figure out the cutoffs later. | 
| 1007 | > |  | 
| 1008 | > | notifyFortranYouAreOnYourOwn(); | 
| 1009 |  |  | 
| 1010 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | 
| 1011 | < | // also send cutoff notification to electrostatics | 
| 912 | < | setElectrostaticCutoffRadius(&rcut_, &rsw_); | 
| 1010 | > | } | 
| 1011 | > | } | 
| 1012 |  | } | 
| 1013 |  |  | 
| 1014 |  | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | 
| 1043 |  | } else { | 
| 1044 |  | // throw error | 
| 1045 |  | sprintf( painCave.errMsg, | 
| 1046 | < | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() ); | 
| 1046 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" | 
| 1047 | > | "\t(Input file specified %s .)\n" | 
| 1048 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" | 
| 1049 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" | 
| 1050 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); | 
| 1051 |  | painCave.isFatal = 1; | 
| 1052 |  | simError(); | 
| 1053 |  | } | 
| 1068 |  | if (!simParams_->haveDampingAlpha()) { | 
| 1069 |  | //throw error | 
| 1070 |  | sprintf( painCave.errMsg, | 
| 1071 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); | 
| 1071 | > | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1072 | > | "\tA default value of %f (1/ang) will be used.\n", alphaVal); | 
| 1073 |  | painCave.isFatal = 0; | 
| 1074 |  | simError(); | 
| 1075 |  | } | 
| 1076 |  | } else { | 
| 1077 |  | // throw error | 
| 1078 |  | sprintf( painCave.errMsg, | 
| 1079 | < | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); | 
| 1079 | > | "SimInfo error: Unknown electrostaticScreeningMethod.\n" | 
| 1080 | > | "\t(Input file specified %s .)\n" | 
| 1081 | > | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" | 
| 1082 | > | "or \"damped\".\n", myScreen.c_str() ); | 
| 1083 |  | painCave.isFatal = 1; | 
| 1084 |  | simError(); | 
| 1085 |  | } | 
| 1087 |  | } | 
| 1088 |  |  | 
| 1089 |  | // let's pass some summation method variables to fortran | 
| 1090 | < | setElectrostaticSummationMethod( &esm ); | 
| 1090 | > | setElectrostaticSumMethod( &esm ); | 
| 1091 | > | setFortranElectrostaticMethod( &esm ); | 
| 1092 |  | setScreeningMethod( &sm ); | 
| 1093 |  | setDampingAlpha( &alphaVal ); | 
| 1094 |  | setReactionFieldDielectric( &dielectric ); | 
| 1095 | < | initFortranFF( &esm, &errorOut ); | 
| 1095 | > | initFortranFF( &errorOut ); | 
| 1096 |  | } | 
| 1097 |  |  | 
| 1098 |  | void SimInfo::setupSwitchingFunction() { |