| 6 |  | * redistribute this software in source and binary code form, provided | 
| 7 |  | * that the following conditions are met: | 
| 8 |  | * | 
| 9 | < | * 1. Acknowledgement of the program authors must be made in any | 
| 10 | < | *    publication of scientific results based in part on use of the | 
| 11 | < | *    program.  An acceptable form of acknowledgement is citation of | 
| 12 | < | *    the article in which the program was described (Matthew | 
| 13 | < | *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher | 
| 14 | < | *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented | 
| 15 | < | *    Parallel Simulation Engine for Molecular Dynamics," | 
| 16 | < | *    J. Comput. Chem. 26, pp. 252-271 (2005)) | 
| 17 | < | * | 
| 18 | < | * 2. Redistributions of source code must retain the above copyright | 
| 9 | > | * 1. Redistributions of source code must retain the above copyright | 
| 10 |  | *    notice, this list of conditions and the following disclaimer. | 
| 11 |  | * | 
| 12 | < | * 3. Redistributions in binary form must reproduce the above copyright | 
| 12 | > | * 2. Redistributions in binary form must reproduce the above copyright | 
| 13 |  | *    notice, this list of conditions and the following disclaimer in the | 
| 14 |  | *    documentation and/or other materials provided with the | 
| 15 |  | *    distribution. | 
| 28 |  | * arising out of the use of or inability to use software, even if the | 
| 29 |  | * University of Notre Dame has been advised of the possibility of | 
| 30 |  | * such damages. | 
| 31 | + | * | 
| 32 | + | * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your | 
| 33 | + | * research, please cite the appropriate papers when you publish your | 
| 34 | + | * work.  Good starting points are: | 
| 35 | + | * | 
| 36 | + | * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | 
| 37 | + | * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | 
| 38 | + | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). | 
| 39 | + | * [4]  Vardeman & Gezelter, in progress (2009). | 
| 40 |  | */ | 
| 41 |  |  | 
| 42 |  | /** | 
| 74 |  | #include "UseTheForce/DarkSide/simParallel_interface.h" | 
| 75 |  | #endif | 
| 76 |  |  | 
| 77 | < | namespace oopse { | 
| 77 | > | namespace OpenMD { | 
| 78 |  | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { | 
| 79 |  | std::map<int, std::set<int> >::iterator i = container.find(index); | 
| 80 |  | std::set<int> result; | 
| 90 |  | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 91 |  | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | 
| 92 |  | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 93 | < | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0), | 
| 94 | < | nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0), | 
| 95 | < | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), | 
| 96 | < | useAtomicVirial_(true) { | 
| 93 | > | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0), | 
| 94 | > | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | 
| 95 | > | nConstraints_(0), sman_(NULL), fortranInitialized_(false), | 
| 96 | > | calcBoxDipole_(false), useAtomicVirial_(true) { | 
| 97 |  |  | 
| 98 | + |  | 
| 99 |  | MoleculeStamp* molStamp; | 
| 100 |  | int nMolWithSameStamp; | 
| 101 |  | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 103 |  | CutoffGroupStamp* cgStamp; | 
| 104 |  | RigidBodyStamp* rbStamp; | 
| 105 |  | int nRigidAtoms = 0; | 
| 106 | + |  | 
| 107 |  | std::vector<Component*> components = simParams->getComponents(); | 
| 108 |  |  | 
| 109 |  | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 158 |  | + nGlobalRigidBodies_; | 
| 159 |  |  | 
| 160 |  | nGlobalMols_ = molStampIds_.size(); | 
| 159 | – |  | 
| 160 | – | #ifdef IS_MPI | 
| 161 |  | molToProcMap_.resize(nGlobalMols_); | 
| 162 | – | #endif | 
| 163 | – |  | 
| 162 |  | } | 
| 163 |  |  | 
| 164 |  | SimInfo::~SimInfo() { | 
| 196 |  | nBonds_ += mol->getNBonds(); | 
| 197 |  | nBends_ += mol->getNBends(); | 
| 198 |  | nTorsions_ += mol->getNTorsions(); | 
| 199 | + | nInversions_ += mol->getNInversions(); | 
| 200 |  | nRigidBodies_ += mol->getNRigidBodies(); | 
| 201 |  | nIntegrableObjects_ += mol->getNIntegrableObjects(); | 
| 202 |  | nCutoffGroups_ += mol->getNCutoffGroups(); | 
| 203 |  | nConstraints_ += mol->getNConstraintPairs(); | 
| 204 |  |  | 
| 205 | < | addExcludePairs(mol); | 
| 206 | < |  | 
| 205 | > | addInteractionPairs(mol); | 
| 206 | > |  | 
| 207 |  | return true; | 
| 208 |  | } else { | 
| 209 |  | return false; | 
| 222 |  | nBonds_ -= mol->getNBonds(); | 
| 223 |  | nBends_ -= mol->getNBends(); | 
| 224 |  | nTorsions_ -= mol->getNTorsions(); | 
| 225 | + | nInversions_ -= mol->getNInversions(); | 
| 226 |  | nRigidBodies_ -= mol->getNRigidBodies(); | 
| 227 |  | nIntegrableObjects_ -= mol->getNIntegrableObjects(); | 
| 228 |  | nCutoffGroups_ -= mol->getNCutoffGroups(); | 
| 229 |  | nConstraints_ -= mol->getNConstraintPairs(); | 
| 230 |  |  | 
| 231 | < | removeExcludePairs(mol); | 
| 231 | > | removeInteractionPairs(mol); | 
| 232 |  | molecules_.erase(mol->getGlobalIndex()); | 
| 233 |  |  | 
| 234 |  | delete mol; | 
| 354 |  |  | 
| 355 |  | } | 
| 356 |  |  | 
| 357 | < | void SimInfo::addExcludePairs(Molecule* mol) { | 
| 357 | > | void SimInfo::addInteractionPairs(Molecule* mol) { | 
| 358 | > | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | 
| 359 |  | std::vector<Bond*>::iterator bondIter; | 
| 360 |  | std::vector<Bend*>::iterator bendIter; | 
| 361 |  | std::vector<Torsion*>::iterator torsionIter; | 
| 362 | + | std::vector<Inversion*>::iterator inversionIter; | 
| 363 |  | Bond* bond; | 
| 364 |  | Bend* bend; | 
| 365 |  | Torsion* torsion; | 
| 366 | + | Inversion* inversion; | 
| 367 |  | int a; | 
| 368 |  | int b; | 
| 369 |  | int c; | 
| 370 |  | int d; | 
| 371 |  |  | 
| 372 | < | std::map<int, std::set<int> > atomGroups; | 
| 372 | > | // atomGroups can be used to add special interaction maps between | 
| 373 | > | // groups of atoms that are in two separate rigid bodies. | 
| 374 | > | // However, most site-site interactions between two rigid bodies | 
| 375 | > | // are probably not special, just the ones between the physically | 
| 376 | > | // bonded atoms.  Interactions *within* a single rigid body should | 
| 377 | > | // always be excluded.  These are done at the bottom of this | 
| 378 | > | // function. | 
| 379 |  |  | 
| 380 | + | std::map<int, std::set<int> > atomGroups; | 
| 381 |  | Molecule::RigidBodyIterator rbIter; | 
| 382 |  | RigidBody* rb; | 
| 383 |  | Molecule::IntegrableObjectIterator ii; | 
| 384 |  | StuntDouble* integrableObject; | 
| 385 |  |  | 
| 386 | < | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 387 | < | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 388 | < |  | 
| 386 | > | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 387 | > | integrableObject != NULL; | 
| 388 | > | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 389 | > |  | 
| 390 |  | if (integrableObject->isRigidBody()) { | 
| 391 | < | rb = static_cast<RigidBody*>(integrableObject); | 
| 392 | < | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 393 | < | std::set<int> rigidAtoms; | 
| 394 | < | for (int i = 0; i < atoms.size(); ++i) { | 
| 395 | < | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 396 | < | } | 
| 397 | < | for (int i = 0; i < atoms.size(); ++i) { | 
| 398 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 399 | < | } | 
| 391 | > | rb = static_cast<RigidBody*>(integrableObject); | 
| 392 | > | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 393 | > | std::set<int> rigidAtoms; | 
| 394 | > | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 395 | > | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 396 | > | } | 
| 397 | > | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 398 | > | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 399 | > | } | 
| 400 |  | } else { | 
| 401 |  | std::set<int> oneAtomSet; | 
| 402 |  | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 403 |  | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); | 
| 404 |  | } | 
| 405 |  | } | 
| 406 | + |  | 
| 407 | + | for (bond= mol->beginBond(bondIter); bond != NULL; | 
| 408 | + | bond = mol->nextBond(bondIter)) { | 
| 409 |  |  | 
| 396 | – |  | 
| 397 | – |  | 
| 398 | – | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 410 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 411 | < | b = bond->getAtomB()->getGlobalIndex(); | 
| 412 | < | exclude_.addPair(a, b); | 
| 411 | > | b = bond->getAtomB()->getGlobalIndex(); | 
| 412 | > |  | 
| 413 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 414 | > | oneTwoInteractions_.addPair(a, b); | 
| 415 | > | } else { | 
| 416 | > | excludedInteractions_.addPair(a, b); | 
| 417 | > | } | 
| 418 |  | } | 
| 419 |  |  | 
| 420 | < | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | 
| 420 | > | for (bend= mol->beginBend(bendIter); bend != NULL; | 
| 421 | > | bend = mol->nextBend(bendIter)) { | 
| 422 | > |  | 
| 423 |  | a = bend->getAtomA()->getGlobalIndex(); | 
| 424 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 425 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 408 | – | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 409 | – | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 410 | – | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 411 | – |  | 
| 412 | – | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 413 | – | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 414 | – | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 426 |  |  | 
| 427 | < | //exclude_.addPair(a, b); | 
| 428 | < | //exclude_.addPair(a, c); | 
| 429 | < | //exclude_.addPair(b, c); | 
| 427 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 428 | > | oneTwoInteractions_.addPair(a, b); | 
| 429 | > | oneTwoInteractions_.addPair(b, c); | 
| 430 | > | } else { | 
| 431 | > | excludedInteractions_.addPair(a, b); | 
| 432 | > | excludedInteractions_.addPair(b, c); | 
| 433 | > | } | 
| 434 | > |  | 
| 435 | > | if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { | 
| 436 | > | oneThreeInteractions_.addPair(a, c); | 
| 437 | > | } else { | 
| 438 | > | excludedInteractions_.addPair(a, c); | 
| 439 | > | } | 
| 440 |  | } | 
| 441 |  |  | 
| 442 | < | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 442 | > | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; | 
| 443 | > | torsion = mol->nextTorsion(torsionIter)) { | 
| 444 | > |  | 
| 445 |  | a = torsion->getAtomA()->getGlobalIndex(); | 
| 446 |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 447 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 448 | < | d = torsion->getAtomD()->getGlobalIndex(); | 
| 426 | < | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 427 | < | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 428 | < | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 429 | < | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 448 | > | d = torsion->getAtomD()->getGlobalIndex(); | 
| 449 |  |  | 
| 450 | < | exclude_.addPairs(rigidSetA, rigidSetB); | 
| 451 | < | exclude_.addPairs(rigidSetA, rigidSetC); | 
| 452 | < | exclude_.addPairs(rigidSetA, rigidSetD); | 
| 453 | < | exclude_.addPairs(rigidSetB, rigidSetC); | 
| 454 | < | exclude_.addPairs(rigidSetB, rigidSetD); | 
| 455 | < | exclude_.addPairs(rigidSetC, rigidSetD); | 
| 450 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 451 | > | oneTwoInteractions_.addPair(a, b); | 
| 452 | > | oneTwoInteractions_.addPair(b, c); | 
| 453 | > | oneTwoInteractions_.addPair(c, d); | 
| 454 | > | } else { | 
| 455 | > | excludedInteractions_.addPair(a, b); | 
| 456 | > | excludedInteractions_.addPair(b, c); | 
| 457 | > | excludedInteractions_.addPair(c, d); | 
| 458 | > | } | 
| 459 |  |  | 
| 460 | < | /* | 
| 461 | < | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 462 | < | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 463 | < | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 464 | < | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 465 | < | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 466 | < | exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 467 | < |  | 
| 468 | < |  | 
| 469 | < | exclude_.addPair(a, b); | 
| 470 | < | exclude_.addPair(a, c); | 
| 471 | < | exclude_.addPair(a, d); | 
| 472 | < | exclude_.addPair(b, c); | 
| 451 | < | exclude_.addPair(b, d); | 
| 452 | < | exclude_.addPair(c, d); | 
| 453 | < | */ | 
| 460 | > | if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { | 
| 461 | > | oneThreeInteractions_.addPair(a, c); | 
| 462 | > | oneThreeInteractions_.addPair(b, d); | 
| 463 | > | } else { | 
| 464 | > | excludedInteractions_.addPair(a, c); | 
| 465 | > | excludedInteractions_.addPair(b, d); | 
| 466 | > | } | 
| 467 | > |  | 
| 468 | > | if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) { | 
| 469 | > | oneFourInteractions_.addPair(a, d); | 
| 470 | > | } else { | 
| 471 | > | excludedInteractions_.addPair(a, d); | 
| 472 | > | } | 
| 473 |  | } | 
| 474 |  |  | 
| 475 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 475 | > | for (inversion= mol->beginInversion(inversionIter); inversion != NULL; | 
| 476 | > | inversion = mol->nextInversion(inversionIter)) { | 
| 477 | > |  | 
| 478 | > | a = inversion->getAtomA()->getGlobalIndex(); | 
| 479 | > | b = inversion->getAtomB()->getGlobalIndex(); | 
| 480 | > | c = inversion->getAtomC()->getGlobalIndex(); | 
| 481 | > | d = inversion->getAtomD()->getGlobalIndex(); | 
| 482 | > |  | 
| 483 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 484 | > | oneTwoInteractions_.addPair(a, b); | 
| 485 | > | oneTwoInteractions_.addPair(a, c); | 
| 486 | > | oneTwoInteractions_.addPair(a, d); | 
| 487 | > | } else { | 
| 488 | > | excludedInteractions_.addPair(a, b); | 
| 489 | > | excludedInteractions_.addPair(a, c); | 
| 490 | > | excludedInteractions_.addPair(a, d); | 
| 491 | > | } | 
| 492 | > |  | 
| 493 | > | if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { | 
| 494 | > | oneThreeInteractions_.addPair(b, c); | 
| 495 | > | oneThreeInteractions_.addPair(b, d); | 
| 496 | > | oneThreeInteractions_.addPair(c, d); | 
| 497 | > | } else { | 
| 498 | > | excludedInteractions_.addPair(b, c); | 
| 499 | > | excludedInteractions_.addPair(b, d); | 
| 500 | > | excludedInteractions_.addPair(c, d); | 
| 501 | > | } | 
| 502 | > | } | 
| 503 | > |  | 
| 504 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 505 | > | rb = mol->nextRigidBody(rbIter)) { | 
| 506 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 507 | < | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 508 | < | for (int j = i + 1; j < atoms.size(); ++j) { | 
| 507 | > | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | 
| 508 | > | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | 
| 509 |  | a = atoms[i]->getGlobalIndex(); | 
| 510 |  | b = atoms[j]->getGlobalIndex(); | 
| 511 | < | exclude_.addPair(a, b); | 
| 511 | > | excludedInteractions_.addPair(a, b); | 
| 512 |  | } | 
| 513 |  | } | 
| 514 |  | } | 
| 515 |  |  | 
| 516 |  | } | 
| 517 |  |  | 
| 518 | < | void SimInfo::removeExcludePairs(Molecule* mol) { | 
| 518 | > | void SimInfo::removeInteractionPairs(Molecule* mol) { | 
| 519 | > | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | 
| 520 |  | std::vector<Bond*>::iterator bondIter; | 
| 521 |  | std::vector<Bend*>::iterator bendIter; | 
| 522 |  | std::vector<Torsion*>::iterator torsionIter; | 
| 523 | + | std::vector<Inversion*>::iterator inversionIter; | 
| 524 |  | Bond* bond; | 
| 525 |  | Bend* bend; | 
| 526 |  | Torsion* torsion; | 
| 527 | + | Inversion* inversion; | 
| 528 |  | int a; | 
| 529 |  | int b; | 
| 530 |  | int c; | 
| 531 |  | int d; | 
| 532 |  |  | 
| 533 |  | std::map<int, std::set<int> > atomGroups; | 
| 482 | – |  | 
| 534 |  | Molecule::RigidBodyIterator rbIter; | 
| 535 |  | RigidBody* rb; | 
| 536 |  | Molecule::IntegrableObjectIterator ii; | 
| 537 |  | StuntDouble* integrableObject; | 
| 538 |  |  | 
| 539 | < | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 540 | < | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 541 | < |  | 
| 539 | > | for (integrableObject = mol->beginIntegrableObject(ii); | 
| 540 | > | integrableObject != NULL; | 
| 541 | > | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 542 | > |  | 
| 543 |  | if (integrableObject->isRigidBody()) { | 
| 544 | < | rb = static_cast<RigidBody*>(integrableObject); | 
| 545 | < | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 546 | < | std::set<int> rigidAtoms; | 
| 547 | < | for (int i = 0; i < atoms.size(); ++i) { | 
| 548 | < | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 549 | < | } | 
| 550 | < | for (int i = 0; i < atoms.size(); ++i) { | 
| 551 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 552 | < | } | 
| 544 | > | rb = static_cast<RigidBody*>(integrableObject); | 
| 545 | > | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 546 | > | std::set<int> rigidAtoms; | 
| 547 | > | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 548 | > | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | 
| 549 | > | } | 
| 550 | > | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | 
| 551 | > | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); | 
| 552 | > | } | 
| 553 |  | } else { | 
| 554 |  | std::set<int> oneAtomSet; | 
| 555 |  | oneAtomSet.insert(integrableObject->getGlobalIndex()); | 
| 557 |  | } | 
| 558 |  | } | 
| 559 |  |  | 
| 560 | < |  | 
| 561 | < | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 560 | > | for (bond= mol->beginBond(bondIter); bond != NULL; | 
| 561 | > | bond = mol->nextBond(bondIter)) { | 
| 562 | > |  | 
| 563 |  | a = bond->getAtomA()->getGlobalIndex(); | 
| 564 | < | b = bond->getAtomB()->getGlobalIndex(); | 
| 565 | < | exclude_.removePair(a, b); | 
| 564 | > | b = bond->getAtomB()->getGlobalIndex(); | 
| 565 | > |  | 
| 566 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 567 | > | oneTwoInteractions_.removePair(a, b); | 
| 568 | > | } else { | 
| 569 | > | excludedInteractions_.removePair(a, b); | 
| 570 | > | } | 
| 571 |  | } | 
| 572 |  |  | 
| 573 | < | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | 
| 573 | > | for (bend= mol->beginBend(bendIter); bend != NULL; | 
| 574 | > | bend = mol->nextBend(bendIter)) { | 
| 575 | > |  | 
| 576 |  | a = bend->getAtomA()->getGlobalIndex(); | 
| 577 |  | b = bend->getAtomB()->getGlobalIndex(); | 
| 578 |  | c = bend->getAtomC()->getGlobalIndex(); | 
| 519 | – |  | 
| 520 | – | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 521 | – | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 522 | – | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 523 | – |  | 
| 524 | – | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 525 | – | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 526 | – | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 579 |  |  | 
| 580 | < | //exclude_.removePair(a, b); | 
| 581 | < | //exclude_.removePair(a, c); | 
| 582 | < | //exclude_.removePair(b, c); | 
| 580 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 581 | > | oneTwoInteractions_.removePair(a, b); | 
| 582 | > | oneTwoInteractions_.removePair(b, c); | 
| 583 | > | } else { | 
| 584 | > | excludedInteractions_.removePair(a, b); | 
| 585 | > | excludedInteractions_.removePair(b, c); | 
| 586 | > | } | 
| 587 | > |  | 
| 588 | > | if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { | 
| 589 | > | oneThreeInteractions_.removePair(a, c); | 
| 590 | > | } else { | 
| 591 | > | excludedInteractions_.removePair(a, c); | 
| 592 | > | } | 
| 593 |  | } | 
| 594 |  |  | 
| 595 | < | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 595 | > | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; | 
| 596 | > | torsion = mol->nextTorsion(torsionIter)) { | 
| 597 | > |  | 
| 598 |  | a = torsion->getAtomA()->getGlobalIndex(); | 
| 599 |  | b = torsion->getAtomB()->getGlobalIndex(); | 
| 600 |  | c = torsion->getAtomC()->getGlobalIndex(); | 
| 601 | < | d = torsion->getAtomD()->getGlobalIndex(); | 
| 601 | > | d = torsion->getAtomD()->getGlobalIndex(); | 
| 602 | > |  | 
| 603 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 604 | > | oneTwoInteractions_.removePair(a, b); | 
| 605 | > | oneTwoInteractions_.removePair(b, c); | 
| 606 | > | oneTwoInteractions_.removePair(c, d); | 
| 607 | > | } else { | 
| 608 | > | excludedInteractions_.removePair(a, b); | 
| 609 | > | excludedInteractions_.removePair(b, c); | 
| 610 | > | excludedInteractions_.removePair(c, d); | 
| 611 | > | } | 
| 612 |  |  | 
| 613 | < | std::set<int> rigidSetA = getRigidSet(a, atomGroups); | 
| 614 | < | std::set<int> rigidSetB = getRigidSet(b, atomGroups); | 
| 615 | < | std::set<int> rigidSetC = getRigidSet(c, atomGroups); | 
| 616 | < | std::set<int> rigidSetD = getRigidSet(d, atomGroups); | 
| 613 | > | if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { | 
| 614 | > | oneThreeInteractions_.removePair(a, c); | 
| 615 | > | oneThreeInteractions_.removePair(b, d); | 
| 616 | > | } else { | 
| 617 | > | excludedInteractions_.removePair(a, c); | 
| 618 | > | excludedInteractions_.removePair(b, d); | 
| 619 | > | } | 
| 620 |  |  | 
| 621 | < | exclude_.removePairs(rigidSetA, rigidSetB); | 
| 622 | < | exclude_.removePairs(rigidSetA, rigidSetC); | 
| 623 | < | exclude_.removePairs(rigidSetA, rigidSetD); | 
| 624 | < | exclude_.removePairs(rigidSetB, rigidSetC); | 
| 625 | < | exclude_.removePairs(rigidSetB, rigidSetD); | 
| 626 | < | exclude_.removePairs(rigidSetC, rigidSetD); | 
| 621 | > | if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) { | 
| 622 | > | oneFourInteractions_.removePair(a, d); | 
| 623 | > | } else { | 
| 624 | > | excludedInteractions_.removePair(a, d); | 
| 625 | > | } | 
| 626 | > | } | 
| 627 |  |  | 
| 628 | < | /* | 
| 629 | < | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); | 
| 553 | < | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 554 | < | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 555 | < | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); | 
| 556 | < | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 557 | < | exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); | 
| 628 | > | for (inversion= mol->beginInversion(inversionIter); inversion != NULL; | 
| 629 | > | inversion = mol->nextInversion(inversionIter)) { | 
| 630 |  |  | 
| 631 | < |  | 
| 632 | < | exclude_.removePair(a, b); | 
| 633 | < | exclude_.removePair(a, c); | 
| 634 | < | exclude_.removePair(a, d); | 
| 563 | < | exclude_.removePair(b, c); | 
| 564 | < | exclude_.removePair(b, d); | 
| 565 | < | exclude_.removePair(c, d); | 
| 566 | < | */ | 
| 567 | < | } | 
| 631 | > | a = inversion->getAtomA()->getGlobalIndex(); | 
| 632 | > | b = inversion->getAtomB()->getGlobalIndex(); | 
| 633 | > | c = inversion->getAtomC()->getGlobalIndex(); | 
| 634 | > | d = inversion->getAtomD()->getGlobalIndex(); | 
| 635 |  |  | 
| 636 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 636 | > | if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { | 
| 637 | > | oneTwoInteractions_.removePair(a, b); | 
| 638 | > | oneTwoInteractions_.removePair(a, c); | 
| 639 | > | oneTwoInteractions_.removePair(a, d); | 
| 640 | > | } else { | 
| 641 | > | excludedInteractions_.removePair(a, b); | 
| 642 | > | excludedInteractions_.removePair(a, c); | 
| 643 | > | excludedInteractions_.removePair(a, d); | 
| 644 | > | } | 
| 645 | > |  | 
| 646 | > | if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { | 
| 647 | > | oneThreeInteractions_.removePair(b, c); | 
| 648 | > | oneThreeInteractions_.removePair(b, d); | 
| 649 | > | oneThreeInteractions_.removePair(c, d); | 
| 650 | > | } else { | 
| 651 | > | excludedInteractions_.removePair(b, c); | 
| 652 | > | excludedInteractions_.removePair(b, d); | 
| 653 | > | excludedInteractions_.removePair(c, d); | 
| 654 | > | } | 
| 655 | > | } | 
| 656 | > |  | 
| 657 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 658 | > | rb = mol->nextRigidBody(rbIter)) { | 
| 659 |  | std::vector<Atom*> atoms = rb->getAtoms(); | 
| 660 | < | for (int i = 0; i < atoms.size() -1 ; ++i) { | 
| 661 | < | for (int j = i + 1; j < atoms.size(); ++j) { | 
| 660 | > | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | 
| 661 | > | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | 
| 662 |  | a = atoms[i]->getGlobalIndex(); | 
| 663 |  | b = atoms[j]->getGlobalIndex(); | 
| 664 | < | exclude_.removePair(a, b); | 
| 664 | > | excludedInteractions_.removePair(a, b); | 
| 665 |  | } | 
| 666 |  | } | 
| 667 |  | } | 
| 668 | < |  | 
| 668 | > |  | 
| 669 |  | } | 
| 670 | < |  | 
| 671 | < |  | 
| 670 | > |  | 
| 671 | > |  | 
| 672 |  | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { | 
| 673 |  | int curStampId; | 
| 674 | < |  | 
| 674 | > |  | 
| 675 |  | //index from 0 | 
| 676 |  | curStampId = moleculeStamps_.size(); | 
| 677 |  |  | 
| 763 |  | useRF = 0; | 
| 764 |  | useSF = 0; | 
| 765 |  | useSP = 0; | 
| 766 | + | useBoxDipole = 0; | 
| 767 |  |  | 
| 768 |  |  | 
| 769 |  | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 888 |  |  | 
| 889 |  | void SimInfo::setupFortranSim() { | 
| 890 |  | int isError; | 
| 891 | < | int nExclude; | 
| 891 | > | int nExclude, nOneTwo, nOneThree, nOneFour; | 
| 892 |  | std::vector<int> fortranGlobalGroupMembership; | 
| 893 |  |  | 
| 804 | – | nExclude = exclude_.getSize(); | 
| 894 |  | isError = 0; | 
| 895 |  |  | 
| 896 |  | //globalGroupMembership_ is filled by SimCreator | 
| 922 |  | else | 
| 923 |  | mfact.push_back( 1.0 ); | 
| 924 |  | } | 
| 836 | – |  | 
| 925 |  | } | 
| 926 |  | } | 
| 927 |  |  | 
| 945 |  | } | 
| 946 |  |  | 
| 947 |  | //setup fortran simulation | 
| 860 | – | int nGlobalExcludes = 0; | 
| 861 | – | int* globalExcludes = NULL; | 
| 862 | – | int* excludeList = exclude_.getExcludeList(); | 
| 863 | – | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList , | 
| 864 | – | &nGlobalExcludes, globalExcludes, &molMembershipArray[0], | 
| 865 | – | &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); | 
| 948 |  |  | 
| 949 | < | if( isError ){ | 
| 949 | > | nExclude = excludedInteractions_.getSize(); | 
| 950 | > | nOneTwo = oneTwoInteractions_.getSize(); | 
| 951 | > | nOneThree = oneThreeInteractions_.getSize(); | 
| 952 | > | nOneFour = oneFourInteractions_.getSize(); | 
| 953 |  |  | 
| 954 | + | int* excludeList = excludedInteractions_.getPairList(); | 
| 955 | + | int* oneTwoList = oneTwoInteractions_.getPairList(); | 
| 956 | + | int* oneThreeList = oneThreeInteractions_.getPairList(); | 
| 957 | + | int* oneFourList = oneFourInteractions_.getPairList(); | 
| 958 | + |  | 
| 959 | + | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], | 
| 960 | + | &nExclude, excludeList, | 
| 961 | + | &nOneTwo, oneTwoList, | 
| 962 | + | &nOneThree, oneThreeList, | 
| 963 | + | &nOneFour, oneFourList, | 
| 964 | + | &molMembershipArray[0], &mfact[0], &nCutoffGroups_, | 
| 965 | + | &fortranGlobalGroupMembership[0], &isError); | 
| 966 | + |  | 
| 967 | + | if( isError ){ | 
| 968 | + |  | 
| 969 |  | sprintf( painCave.errMsg, | 
| 970 |  | "There was an error setting the simulation information in fortran.\n" ); | 
| 971 |  | painCave.isFatal = 1; | 
| 972 | < | painCave.severity = OOPSE_ERROR; | 
| 972 | > | painCave.severity = OPENMD_ERROR; | 
| 973 |  | simError(); | 
| 974 |  | } | 
| 975 | < |  | 
| 976 | < | #ifdef IS_MPI | 
| 975 | > |  | 
| 976 | > |  | 
| 977 |  | sprintf( checkPointMsg, | 
| 978 |  | "succesfully sent the simulation information to fortran.\n"); | 
| 979 | < | MPIcheckPoint(); | 
| 980 | < | #endif // is_mpi | 
| 981 | < |  | 
| 979 | > |  | 
| 980 | > | errorCheckPoint(); | 
| 981 | > |  | 
| 982 |  | // Setup number of neighbors in neighbor list if present | 
| 983 |  | if (simParams_->haveNeighborListNeighbors()) { | 
| 984 |  | int nlistNeighbors = simParams_->getNeighborListNeighbors(); | 
| 989 |  | } | 
| 990 |  |  | 
| 991 |  |  | 
| 892 | – | #ifdef IS_MPI | 
| 992 |  | void SimInfo::setupFortranParallel() { | 
| 993 | < |  | 
| 993 | > | #ifdef IS_MPI | 
| 994 |  | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex | 
| 995 |  | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); | 
| 996 |  | std::vector<int> localToGlobalCutoffGroupIndex; | 
| 1040 |  | } | 
| 1041 |  |  | 
| 1042 |  | sprintf(checkPointMsg, " mpiRefresh successful.\n"); | 
| 1043 | < | MPIcheckPoint(); | 
| 1043 | > | errorCheckPoint(); | 
| 1044 |  |  | 
| 1045 | < |  | 
| 1045 | > | #endif | 
| 1046 |  | } | 
| 1047 |  |  | 
| 949 | – | #endif | 
| 950 | – |  | 
| 1048 |  | void SimInfo::setupCutoff() { | 
| 1049 |  |  | 
| 1050 |  | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | 
| 1053 |  | int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | 
| 1054 |  |  | 
| 1055 |  | // Set LJ shifting bools to false | 
| 1056 | < | ljsp_ = false; | 
| 1057 | < | ljsf_ = false; | 
| 1056 | > | ljsp_ = 0; | 
| 1057 | > | ljsf_ = 0; | 
| 1058 |  |  | 
| 1059 |  | std::string myPolicy; | 
| 1060 |  | if (forceFieldOptions_.haveCutoffPolicy()){ | 
| 1105 |  | rsw_ = 0.85 * rcut_; | 
| 1106 |  | sprintf(painCave.errMsg, | 
| 1107 |  | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 1108 | < | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" | 
| 1108 | > | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" | 
| 1109 |  | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1110 |  | painCave.isFatal = 0; | 
| 1111 |  | simError(); | 
| 1113 |  | rsw_ = rcut_; | 
| 1114 |  | sprintf(painCave.errMsg, | 
| 1115 |  | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 1116 | < | "\tOOPSE will use the same value as the cutoffRadius.\n" | 
| 1116 | > | "\tOpenMD will use the same value as the cutoffRadius.\n" | 
| 1117 |  | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1118 |  | painCave.isFatal = 0; | 
| 1119 |  | simError(); | 
| 1125 |  | toUpper(myMethod); | 
| 1126 |  |  | 
| 1127 |  | if (myMethod == "SHIFTED_POTENTIAL") { | 
| 1128 | < | ljsp_ = true; | 
| 1128 | > | ljsp_ = 1; | 
| 1129 |  | } else if (myMethod == "SHIFTED_FORCE") { | 
| 1130 | < | ljsf_ = true; | 
| 1130 | > | ljsf_ = 1; | 
| 1131 |  | } | 
| 1132 |  | } | 
| 1133 | + |  | 
| 1134 |  | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); | 
| 1135 |  |  | 
| 1136 |  | } else { | 
| 1139 |  | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 1140 |  | sprintf(painCave.errMsg, | 
| 1141 |  | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 1142 | < | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 1142 | > | "\tOpenMD will use a default value of 15.0 angstroms" | 
| 1143 |  | "\tfor the cutoffRadius.\n"); | 
| 1144 |  | painCave.isFatal = 0; | 
| 1145 |  | simError(); | 
| 1152 |  | // For the time being, we're tethering the LJ shifted behavior to the | 
| 1153 |  | // electrostaticSummationMethod keyword options | 
| 1154 |  | if (myMethod == "SHIFTED_POTENTIAL") { | 
| 1155 | < | ljsp_ = true; | 
| 1155 | > | ljsp_ = 1; | 
| 1156 |  | } else if (myMethod == "SHIFTED_FORCE") { | 
| 1157 | < | ljsf_ = true; | 
| 1157 | > | ljsf_ = 1; | 
| 1158 |  | } | 
| 1159 |  | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 1160 |  | if (simParams_->haveSwitchingRadius()){ | 
| 1173 |  | } else { | 
| 1174 |  | sprintf(painCave.errMsg, | 
| 1175 |  | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 1176 | < | "\tOOPSE will use a default value of\n" | 
| 1176 | > | "\tOpenMD will use a default value of\n" | 
| 1177 |  | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 1178 |  | painCave.isFatal = 0; | 
| 1179 |  | simError(); | 
| 1634 |  | return sdByGlobalIndex_.at(index); | 
| 1635 |  | } | 
| 1636 |  | */ | 
| 1637 | < | }//end namespace oopse | 
| 1637 | > | }//end namespace OpenMD | 
| 1638 |  |  |