| 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" | 
| 65 | 
  | 
#include "selection/SelectionManager.hpp" | 
| 66 | 
+ | 
#include "io/ForceFieldOptions.hpp" | 
| 67 | 
+ | 
#include "UseTheForce/ForceField.hpp" | 
| 68 | 
  | 
 | 
| 69 | 
  | 
#ifdef IS_MPI | 
| 70 | 
  | 
#include "UseTheForce/mpiComponentPlan.h" | 
| 72 | 
  | 
#endif  | 
| 73 | 
  | 
 | 
| 74 | 
  | 
namespace oopse { | 
| 75 | 
+ | 
  std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { | 
| 76 | 
+ | 
    std::map<int, std::set<int> >::iterator i = container.find(index); | 
| 77 | 
+ | 
    std::set<int> result; | 
| 78 | 
+ | 
    if (i != container.end()) { | 
| 79 | 
+ | 
        result = i->second; | 
| 80 | 
+ | 
    } | 
| 81 | 
  | 
 | 
| 82 | 
< | 
  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,  | 
| 83 | 
< | 
                   ForceField* ff, Globals* simParams) :  | 
| 84 | 
< | 
    stamps_(stamps), forceField_(ff), simParams_(simParams),  | 
| 82 | 
> | 
    return result; | 
| 83 | 
> | 
  } | 
| 84 | 
> | 
   | 
| 85 | 
> | 
  SimInfo::SimInfo(ForceField* ff, Globals* simParams) :  | 
| 86 | 
> | 
    forceField_(ff), simParams_(simParams),  | 
| 87 | 
  | 
    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 88 | 
  | 
    nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),  | 
| 89 | 
  | 
    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 91 | 
  | 
    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0), | 
| 92 | 
  | 
    sman_(NULL), fortranInitialized_(false) { | 
| 93 | 
  | 
 | 
| 84 | 
– | 
             | 
| 85 | 
– | 
      std::vector<std::pair<MoleculeStamp*, int> >::iterator i; | 
| 94 | 
  | 
      MoleculeStamp* molStamp; | 
| 95 | 
  | 
      int nMolWithSameStamp; | 
| 96 | 
  | 
      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 98 | 
  | 
      CutoffGroupStamp* cgStamp;     | 
| 99 | 
  | 
      RigidBodyStamp* rbStamp; | 
| 100 | 
  | 
      int nRigidAtoms = 0; | 
| 101 | 
< | 
     | 
| 102 | 
< | 
      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { | 
| 103 | 
< | 
        molStamp = i->first; | 
| 104 | 
< | 
        nMolWithSameStamp = i->second; | 
| 101 | 
> | 
      std::vector<Component*> components = simParams->getComponents(); | 
| 102 | 
> | 
       | 
| 103 | 
> | 
      for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 104 | 
> | 
        molStamp = (*i)->getMoleculeStamp(); | 
| 105 | 
> | 
        nMolWithSameStamp = (*i)->getNMol(); | 
| 106 | 
  | 
         | 
| 107 | 
  | 
        addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 108 | 
  | 
 | 
| 109 | 
  | 
        //calculate atoms in molecules | 
| 110 | 
  | 
        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;    | 
| 111 | 
  | 
 | 
| 103 | 
– | 
 | 
| 112 | 
  | 
        //calculate atoms in cutoff groups | 
| 113 | 
  | 
        int nAtomsInGroups = 0; | 
| 114 | 
  | 
        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 115 | 
  | 
         | 
| 116 | 
  | 
        for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 117 | 
< | 
          cgStamp = molStamp->getCutoffGroup(j); | 
| 117 | 
> | 
          cgStamp = molStamp->getCutoffGroupStamp(j); | 
| 118 | 
  | 
          nAtomsInGroups += cgStamp->getNMembers(); | 
| 119 | 
  | 
        } | 
| 120 | 
  | 
 | 
| 127 | 
  | 
        int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 128 | 
  | 
         | 
| 129 | 
  | 
        for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 130 | 
< | 
          rbStamp = molStamp->getRigidBody(j); | 
| 130 | 
> | 
          rbStamp = molStamp->getRigidBodyStamp(j); | 
| 131 | 
  | 
          nAtomsInRigidBodies += rbStamp->getNMembers(); | 
| 132 | 
  | 
        } | 
| 133 | 
  | 
 | 
| 166 | 
  | 
    } | 
| 167 | 
  | 
    molecules_.clear(); | 
| 168 | 
  | 
        | 
| 161 | 
– | 
    delete stamps_; | 
| 169 | 
  | 
    delete sman_; | 
| 170 | 
  | 
    delete simParams_; | 
| 171 | 
  | 
    delete forceField_; | 
| 272 | 
  | 
          } | 
| 273 | 
  | 
        } | 
| 274 | 
  | 
             | 
| 275 | 
< | 
      }//end for (integrableObject) | 
| 276 | 
< | 
    }// end for (mol) | 
| 275 | 
> | 
      } | 
| 276 | 
> | 
    } | 
| 277 | 
  | 
     | 
| 278 | 
  | 
    // n_constraints is local, so subtract them on each processor | 
| 279 | 
  | 
    ndf_local -= nConstraints_; | 
| 352 | 
  | 
    int b; | 
| 353 | 
  | 
    int c; | 
| 354 | 
  | 
    int d; | 
| 355 | 
+ | 
 | 
| 356 | 
+ | 
    std::map<int, std::set<int> > atomGroups; | 
| 357 | 
+ | 
 | 
| 358 | 
+ | 
    Molecule::RigidBodyIterator rbIter; | 
| 359 | 
+ | 
    RigidBody* rb; | 
| 360 | 
+ | 
    Molecule::IntegrableObjectIterator ii; | 
| 361 | 
+ | 
    StuntDouble* integrableObject; | 
| 362 | 
  | 
     | 
| 363 | 
+ | 
    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 364 | 
+ | 
           integrableObject = mol->nextIntegrableObject(ii)) { | 
| 365 | 
+ | 
 | 
| 366 | 
+ | 
      if (integrableObject->isRigidBody()) { | 
| 367 | 
+ | 
          rb = static_cast<RigidBody*>(integrableObject); | 
| 368 | 
+ | 
          std::vector<Atom*> atoms = rb->getAtoms(); | 
| 369 | 
+ | 
          std::set<int> rigidAtoms; | 
| 370 | 
+ | 
          for (int i = 0; i < atoms.size(); ++i) { | 
| 371 | 
+ | 
            rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 372 | 
+ | 
          } | 
| 373 | 
+ | 
          for (int i = 0; i < atoms.size(); ++i) { | 
| 374 | 
+ | 
            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 375 | 
+ | 
          }       | 
| 376 | 
+ | 
      } else { | 
| 377 | 
+ | 
        std::set<int> oneAtomSet; | 
| 378 | 
+ | 
        oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 379 | 
+ | 
        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));         | 
| 380 | 
+ | 
      } | 
| 381 | 
+ | 
    }   | 
| 382 | 
+ | 
 | 
| 383 | 
+ | 
     | 
| 384 | 
+ | 
     | 
| 385 | 
  | 
    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 386 | 
  | 
      a = bond->getAtomA()->getGlobalIndex(); | 
| 387 | 
  | 
      b = bond->getAtomB()->getGlobalIndex();         | 
| 392 | 
  | 
      a = bend->getAtomA()->getGlobalIndex(); | 
| 393 | 
  | 
      b = bend->getAtomB()->getGlobalIndex();         | 
| 394 | 
  | 
      c = bend->getAtomC()->getGlobalIndex(); | 
| 395 | 
+ | 
      std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 396 | 
+ | 
      std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 397 | 
+ | 
      std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 398 | 
  | 
 | 
| 399 | 
< | 
      exclude_.addPair(a, b); | 
| 400 | 
< | 
      exclude_.addPair(a, c); | 
| 401 | 
< | 
      exclude_.addPair(b, c);         | 
| 399 | 
> | 
      exclude_.addPairs(rigidSetA, rigidSetB); | 
| 400 | 
> | 
      exclude_.addPairs(rigidSetA, rigidSetC); | 
| 401 | 
> | 
      exclude_.addPairs(rigidSetB, rigidSetC); | 
| 402 | 
> | 
       | 
| 403 | 
> | 
      //exclude_.addPair(a, b); | 
| 404 | 
> | 
      //exclude_.addPair(a, c); | 
| 405 | 
> | 
      //exclude_.addPair(b, c);         | 
| 406 | 
  | 
    } | 
| 407 | 
  | 
 | 
| 408 | 
  | 
    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 410 | 
  | 
      b = torsion->getAtomB()->getGlobalIndex();         | 
| 411 | 
  | 
      c = torsion->getAtomC()->getGlobalIndex();         | 
| 412 | 
  | 
      d = torsion->getAtomD()->getGlobalIndex();         | 
| 413 | 
+ | 
      std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 414 | 
+ | 
      std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 415 | 
+ | 
      std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 416 | 
+ | 
      std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 417 | 
  | 
 | 
| 418 | 
+ | 
      exclude_.addPairs(rigidSetA, rigidSetB); | 
| 419 | 
+ | 
      exclude_.addPairs(rigidSetA, rigidSetC); | 
| 420 | 
+ | 
      exclude_.addPairs(rigidSetA, rigidSetD); | 
| 421 | 
+ | 
      exclude_.addPairs(rigidSetB, rigidSetC); | 
| 422 | 
+ | 
      exclude_.addPairs(rigidSetB, rigidSetD); | 
| 423 | 
+ | 
      exclude_.addPairs(rigidSetC, rigidSetD); | 
| 424 | 
+ | 
 | 
| 425 | 
+ | 
      /* | 
| 426 | 
+ | 
      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 427 | 
+ | 
      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 428 | 
+ | 
      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 429 | 
+ | 
      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 430 | 
+ | 
      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 431 | 
+ | 
      exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 432 | 
+ | 
          | 
| 433 | 
+ | 
       | 
| 434 | 
  | 
      exclude_.addPair(a, b); | 
| 435 | 
  | 
      exclude_.addPair(a, c); | 
| 436 | 
  | 
      exclude_.addPair(a, d); | 
| 437 | 
  | 
      exclude_.addPair(b, c); | 
| 438 | 
  | 
      exclude_.addPair(b, d); | 
| 439 | 
  | 
      exclude_.addPair(c, d);         | 
| 440 | 
+ | 
      */ | 
| 441 | 
  | 
    } | 
| 442 | 
  | 
 | 
| 379 | 
– | 
    Molecule::RigidBodyIterator rbIter; | 
| 380 | 
– | 
    RigidBody* rb; | 
| 443 | 
  | 
    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 444 | 
  | 
      std::vector<Atom*> atoms = rb->getAtoms(); | 
| 445 | 
  | 
      for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 464 | 
  | 
    int b; | 
| 465 | 
  | 
    int c; | 
| 466 | 
  | 
    int d; | 
| 467 | 
+ | 
 | 
| 468 | 
+ | 
    std::map<int, std::set<int> > atomGroups; | 
| 469 | 
+ | 
 | 
| 470 | 
+ | 
    Molecule::RigidBodyIterator rbIter; | 
| 471 | 
+ | 
    RigidBody* rb; | 
| 472 | 
+ | 
    Molecule::IntegrableObjectIterator ii; | 
| 473 | 
+ | 
    StuntDouble* integrableObject; | 
| 474 | 
  | 
     | 
| 475 | 
+ | 
    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 476 | 
+ | 
           integrableObject = mol->nextIntegrableObject(ii)) { | 
| 477 | 
+ | 
 | 
| 478 | 
+ | 
      if (integrableObject->isRigidBody()) { | 
| 479 | 
+ | 
          rb = static_cast<RigidBody*>(integrableObject); | 
| 480 | 
+ | 
          std::vector<Atom*> atoms = rb->getAtoms(); | 
| 481 | 
+ | 
          std::set<int> rigidAtoms; | 
| 482 | 
+ | 
          for (int i = 0; i < atoms.size(); ++i) { | 
| 483 | 
+ | 
            rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 484 | 
+ | 
          } | 
| 485 | 
+ | 
          for (int i = 0; i < atoms.size(); ++i) { | 
| 486 | 
+ | 
            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 487 | 
+ | 
          }       | 
| 488 | 
+ | 
      } else { | 
| 489 | 
+ | 
        std::set<int> oneAtomSet; | 
| 490 | 
+ | 
        oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 491 | 
+ | 
        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));         | 
| 492 | 
+ | 
      } | 
| 493 | 
+ | 
    }   | 
| 494 | 
+ | 
 | 
| 495 | 
+ | 
     | 
| 496 | 
  | 
    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 497 | 
  | 
      a = bond->getAtomA()->getGlobalIndex(); | 
| 498 | 
  | 
      b = bond->getAtomB()->getGlobalIndex();         | 
| 504 | 
  | 
      b = bend->getAtomB()->getGlobalIndex();         | 
| 505 | 
  | 
      c = bend->getAtomC()->getGlobalIndex(); | 
| 506 | 
  | 
 | 
| 507 | 
< | 
      exclude_.removePair(a, b); | 
| 508 | 
< | 
      exclude_.removePair(a, c); | 
| 509 | 
< | 
      exclude_.removePair(b, c);         | 
| 507 | 
> | 
      std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 508 | 
> | 
      std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 509 | 
> | 
      std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 510 | 
> | 
 | 
| 511 | 
> | 
      exclude_.removePairs(rigidSetA, rigidSetB); | 
| 512 | 
> | 
      exclude_.removePairs(rigidSetA, rigidSetC); | 
| 513 | 
> | 
      exclude_.removePairs(rigidSetB, rigidSetC); | 
| 514 | 
> | 
       | 
| 515 | 
> | 
      //exclude_.removePair(a, b); | 
| 516 | 
> | 
      //exclude_.removePair(a, c); | 
| 517 | 
> | 
      //exclude_.removePair(b, c);         | 
| 518 | 
  | 
    } | 
| 519 | 
  | 
 | 
| 520 | 
  | 
    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 523 | 
  | 
      c = torsion->getAtomC()->getGlobalIndex();         | 
| 524 | 
  | 
      d = torsion->getAtomD()->getGlobalIndex();         | 
| 525 | 
  | 
 | 
| 526 | 
+ | 
      std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 527 | 
+ | 
      std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 528 | 
+ | 
      std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 529 | 
+ | 
      std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 530 | 
+ | 
 | 
| 531 | 
+ | 
      exclude_.removePairs(rigidSetA, rigidSetB); | 
| 532 | 
+ | 
      exclude_.removePairs(rigidSetA, rigidSetC); | 
| 533 | 
+ | 
      exclude_.removePairs(rigidSetA, rigidSetD); | 
| 534 | 
+ | 
      exclude_.removePairs(rigidSetB, rigidSetC); | 
| 535 | 
+ | 
      exclude_.removePairs(rigidSetB, rigidSetD); | 
| 536 | 
+ | 
      exclude_.removePairs(rigidSetC, rigidSetD); | 
| 537 | 
+ | 
 | 
| 538 | 
+ | 
      /* | 
| 539 | 
+ | 
      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 540 | 
+ | 
      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 541 | 
+ | 
      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 542 | 
+ | 
      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 543 | 
+ | 
      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 544 | 
+ | 
      exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 545 | 
+ | 
 | 
| 546 | 
+ | 
       | 
| 547 | 
  | 
      exclude_.removePair(a, b); | 
| 548 | 
  | 
      exclude_.removePair(a, c); | 
| 549 | 
  | 
      exclude_.removePair(a, d); | 
| 550 | 
  | 
      exclude_.removePair(b, c); | 
| 551 | 
  | 
      exclude_.removePair(b, d); | 
| 552 | 
  | 
      exclude_.removePair(c, d);         | 
| 553 | 
+ | 
      */ | 
| 554 | 
  | 
    } | 
| 555 | 
  | 
 | 
| 436 | 
– | 
    Molecule::RigidBodyIterator rbIter; | 
| 437 | 
– | 
    RigidBody* rb; | 
| 556 | 
  | 
    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 557 | 
  | 
      std::vector<Atom*> atoms = rb->getAtoms(); | 
| 558 | 
  | 
      for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 919 | 
  | 
  } | 
| 920 | 
  | 
 | 
| 921 | 
  | 
#endif | 
| 804 | 
– | 
 | 
| 805 | 
– | 
  double SimInfo::calcMaxCutoffRadius() { | 
| 922 | 
  | 
 | 
| 923 | 
+ | 
  void SimInfo::setupCutoff() {            | 
| 924 | 
+ | 
     | 
| 925 | 
+ | 
    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | 
| 926 | 
  | 
 | 
| 927 | 
< | 
    std::set<AtomType*> atomTypes; | 
| 928 | 
< | 
    std::set<AtomType*>::iterator i; | 
| 810 | 
< | 
    std::vector<double> cutoffRadius; | 
| 927 | 
> | 
    // Check the cutoff policy | 
| 928 | 
> | 
    int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | 
| 929 | 
  | 
 | 
| 930 | 
< | 
    //get the unique atom types | 
| 931 | 
< | 
    atomTypes = getUniqueAtomTypes(); | 
| 932 | 
< | 
 | 
| 933 | 
< | 
    //query the max cutoff radius among these atom types | 
| 934 | 
< | 
    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 817 | 
< | 
      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); | 
| 930 | 
> | 
    std::string myPolicy; | 
| 931 | 
> | 
    if (forceFieldOptions_.haveCutoffPolicy()){ | 
| 932 | 
> | 
      myPolicy = forceFieldOptions_.getCutoffPolicy(); | 
| 933 | 
> | 
    }else if (simParams_->haveCutoffPolicy()) { | 
| 934 | 
> | 
      myPolicy = simParams_->getCutoffPolicy(); | 
| 935 | 
  | 
    } | 
| 936 | 
  | 
 | 
| 937 | 
< | 
    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); | 
| 821 | 
< | 
#ifdef IS_MPI | 
| 822 | 
< | 
    //pick the max cutoff radius among the processors | 
| 823 | 
< | 
#endif | 
| 824 | 
< | 
 | 
| 825 | 
< | 
    return maxCutoffRadius; | 
| 826 | 
< | 
  } | 
| 827 | 
< | 
 | 
| 828 | 
< | 
  void SimInfo::getCutoff(double& rcut, double& rsw) { | 
| 829 | 
< | 
     | 
| 830 | 
< | 
    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 | 
< | 
     | 
| 882 | 
< | 
    int cp =  TRADITIONAL_CUTOFF_POLICY; | 
| 883 | 
< | 
    if (simParams_->haveCutoffPolicy()) { | 
| 884 | 
< | 
      std::string myPolicy = simParams_->getCutoffPolicy(); | 
| 937 | 
> | 
    if (!myPolicy.empty()){ | 
| 938 | 
  | 
      toUpper(myPolicy); | 
| 939 | 
  | 
      if (myPolicy == "MIX") { | 
| 940 | 
  | 
        cp = MIX_CUTOFF_POLICY; | 
| 953 | 
  | 
          }      | 
| 954 | 
  | 
        }            | 
| 955 | 
  | 
      } | 
| 956 | 
< | 
    } | 
| 956 | 
> | 
    }            | 
| 957 | 
> | 
    notifyFortranCutoffPolicy(&cp); | 
| 958 | 
  | 
 | 
| 959 | 
< | 
 | 
| 959 | 
> | 
    // Check the Skin Thickness for neighborlists | 
| 960 | 
> | 
    double skin; | 
| 961 | 
  | 
    if (simParams_->haveSkinThickness()) { | 
| 962 | 
< | 
      double skinThickness = simParams_->getSkinThickness(); | 
| 963 | 
< | 
    } | 
| 962 | 
> | 
      skin = simParams_->getSkinThickness(); | 
| 963 | 
> | 
      notifyFortranSkinThickness(&skin); | 
| 964 | 
> | 
    }             | 
| 965 | 
> | 
         | 
| 966 | 
> | 
    // Check if the cutoff was set explicitly: | 
| 967 | 
> | 
    if (simParams_->haveCutoffRadius()) { | 
| 968 | 
> | 
      rcut_ = simParams_->getCutoffRadius(); | 
| 969 | 
> | 
      if (simParams_->haveSwitchingRadius()) { | 
| 970 | 
> | 
        rsw_  = simParams_->getSwitchingRadius(); | 
| 971 | 
> | 
      } else { | 
| 972 | 
> | 
        if (fInfo_.SIM_uses_Charges |  | 
| 973 | 
> | 
            fInfo_.SIM_uses_Dipoles |  | 
| 974 | 
> | 
            fInfo_.SIM_uses_RF) { | 
| 975 | 
> | 
           | 
| 976 | 
> | 
          rsw_ = 0.85 * rcut_; | 
| 977 | 
> | 
          sprintf(painCave.errMsg, | 
| 978 | 
> | 
                  "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 979 | 
> | 
                  "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" | 
| 980 | 
> | 
                  "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 981 | 
> | 
        painCave.isFatal = 0; | 
| 982 | 
> | 
        simError(); | 
| 983 | 
> | 
        } else { | 
| 984 | 
> | 
          rsw_ = rcut_; | 
| 985 | 
> | 
          sprintf(painCave.errMsg, | 
| 986 | 
> | 
                  "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 987 | 
> | 
                  "\tOOPSE will use the same value as the cutoffRadius.\n" | 
| 988 | 
> | 
                  "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 989 | 
> | 
          painCave.isFatal = 0; | 
| 990 | 
> | 
          simError(); | 
| 991 | 
> | 
        } | 
| 992 | 
> | 
      } | 
| 993 | 
> | 
       | 
| 994 | 
> | 
      notifyFortranCutoffs(&rcut_, &rsw_); | 
| 995 | 
> | 
       | 
| 996 | 
> | 
    } else { | 
| 997 | 
> | 
       | 
| 998 | 
> | 
      // For electrostatic atoms, we'll assume a large safe value: | 
| 999 | 
> | 
      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 1000 | 
> | 
        sprintf(painCave.errMsg, | 
| 1001 | 
> | 
                "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 1002 | 
> | 
                "\tOOPSE will use a default value of 15.0 angstroms" | 
| 1003 | 
> | 
                "\tfor the cutoffRadius.\n"); | 
| 1004 | 
> | 
        painCave.isFatal = 0; | 
| 1005 | 
> | 
        simError(); | 
| 1006 | 
> | 
        rcut_ = 15.0; | 
| 1007 | 
> | 
       | 
| 1008 | 
> | 
        if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1009 | 
> | 
          std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1010 | 
> | 
          toUpper(myMethod); | 
| 1011 | 
> | 
          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 1012 | 
> | 
            if (simParams_->haveSwitchingRadius()){ | 
| 1013 | 
> | 
              sprintf(painCave.errMsg, | 
| 1014 | 
> | 
                      "SimInfo Warning: A value was set for the switchingRadius\n" | 
| 1015 | 
> | 
                      "\teven though the electrostaticSummationMethod was\n" | 
| 1016 | 
> | 
                      "\tset to %s\n", myMethod.c_str()); | 
| 1017 | 
> | 
              painCave.isFatal = 1; | 
| 1018 | 
> | 
              simError();             | 
| 1019 | 
> | 
            }  | 
| 1020 | 
> | 
          } | 
| 1021 | 
> | 
        } | 
| 1022 | 
> | 
       | 
| 1023 | 
> | 
        if (simParams_->haveSwitchingRadius()){ | 
| 1024 | 
> | 
          rsw_ = simParams_->getSwitchingRadius(); | 
| 1025 | 
> | 
        } else {         | 
| 1026 | 
> | 
          sprintf(painCave.errMsg, | 
| 1027 | 
> | 
                  "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 1028 | 
> | 
                  "\tOOPSE will use a default value of\n" | 
| 1029 | 
> | 
                  "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 1030 | 
> | 
          painCave.isFatal = 0; | 
| 1031 | 
> | 
          simError(); | 
| 1032 | 
> | 
          rsw_ = 0.85 * rcut_; | 
| 1033 | 
> | 
        } | 
| 1034 | 
> | 
        notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1035 | 
> | 
      } else { | 
| 1036 | 
> | 
        // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | 
| 1037 | 
> | 
        // We'll punt and let fortran figure out the cutoffs later. | 
| 1038 | 
> | 
         | 
| 1039 | 
> | 
        notifyFortranYouAreOnYourOwn(); | 
| 1040 | 
  | 
 | 
| 1041 | 
< | 
    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | 
| 1042 | 
< | 
    // also send cutoff notification to electrostatics | 
| 912 | 
< | 
    setElectrostaticCutoffRadius(&rcut_, &rsw_); | 
| 1041 | 
> | 
      } | 
| 1042 | 
> | 
    } | 
| 1043 | 
  | 
  } | 
| 1044 | 
  | 
 | 
| 1045 | 
  | 
  void SimInfo::setupElectrostaticSummationMethod( int isError ) {     | 
| 1074 | 
  | 
              } else { | 
| 1075 | 
  | 
                // throw error         | 
| 1076 | 
  | 
                sprintf( painCave.errMsg, | 
| 1077 | 
< | 
                         "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() ); | 
| 1077 | 
> | 
                         "SimInfo error: Unknown electrostaticSummationMethod.\n" | 
| 1078 | 
> | 
                         "\t(Input file specified %s .)\n" | 
| 1079 | 
> | 
                         "\telectrostaticSummationMethod must be one of: \"none\",\n" | 
| 1080 | 
> | 
                         "\t\"shifted_potential\", \"shifted_force\", or \n" | 
| 1081 | 
> | 
                         "\t\"reaction_field\".\n", myMethod.c_str() ); | 
| 1082 | 
  | 
                painCave.isFatal = 1; | 
| 1083 | 
  | 
                simError(); | 
| 1084 | 
  | 
              }      | 
| 1099 | 
  | 
          if (!simParams_->haveDampingAlpha()) { | 
| 1100 | 
  | 
            //throw error | 
| 1101 | 
  | 
            sprintf( painCave.errMsg, | 
| 1102 | 
< | 
                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); | 
| 1102 | 
> | 
                     "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1103 | 
> | 
                     "\tA default value of %f (1/ang) will be used.\n", alphaVal); | 
| 1104 | 
  | 
            painCave.isFatal = 0; | 
| 1105 | 
  | 
            simError(); | 
| 1106 | 
  | 
          } | 
| 1107 | 
  | 
        } else { | 
| 1108 | 
  | 
          // throw error         | 
| 1109 | 
  | 
          sprintf( painCave.errMsg, | 
| 1110 | 
< | 
                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); | 
| 1110 | 
> | 
                   "SimInfo error: Unknown electrostaticScreeningMethod.\n" | 
| 1111 | 
> | 
                   "\t(Input file specified %s .)\n" | 
| 1112 | 
> | 
                   "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" | 
| 1113 | 
> | 
                   "or \"damped\".\n", myScreen.c_str() ); | 
| 1114 | 
  | 
          painCave.isFatal = 1; | 
| 1115 | 
  | 
          simError(); | 
| 1116 | 
  | 
        } | 
| 1119 | 
  | 
     | 
| 1120 | 
  | 
    // let's pass some summation method variables to fortran | 
| 1121 | 
  | 
    setElectrostaticSummationMethod( &esm ); | 
| 1122 | 
+ | 
    setFortranElectrostaticMethod( &esm ); | 
| 1123 | 
  | 
    setScreeningMethod( &sm ); | 
| 1124 | 
  | 
    setDampingAlpha( &alphaVal ); | 
| 1125 | 
  | 
    setReactionFieldDielectric( &dielectric ); | 
| 1126 | 
< | 
    initFortranFF( &esm, &errorOut ); | 
| 1126 | 
> | 
    initFortranFF( &errorOut ); | 
| 1127 | 
  | 
  } | 
| 1128 | 
  | 
 | 
| 1129 | 
  | 
  void SimInfo::setupSwitchingFunction() {     |