| 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).                         | 
| 39 | 
> | 
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | 
> | 
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 | 
  | 
 */ | 
| 42 | 
  | 
  | 
| 43 | 
  | 
/** | 
| 56 | 
  | 
#include "brains/SimCreator.hpp" | 
| 57 | 
  | 
#include "brains/SimSnapshotManager.hpp" | 
| 58 | 
  | 
#include "io/DumpReader.hpp" | 
| 59 | 
< | 
#include "UseTheForce/ForceFieldFactory.hpp" | 
| 59 | 
> | 
#include "brains/ForceField.hpp" | 
| 60 | 
  | 
#include "utils/simError.h" | 
| 61 | 
  | 
#include "utils/StringUtils.hpp" | 
| 62 | 
  | 
#include "math/SeqRandNumGen.hpp" | 
| 76 | 
  | 
#include "antlr/NoViableAltForCharException.hpp" | 
| 77 | 
  | 
#include "antlr/NoViableAltException.hpp" | 
| 78 | 
  | 
 | 
| 79 | 
+ | 
#include "types/DirectionalAdapter.hpp" | 
| 80 | 
+ | 
#include "types/MultipoleAdapter.hpp" | 
| 81 | 
+ | 
#include "types/EAMAdapter.hpp" | 
| 82 | 
+ | 
#include "types/SuttonChenAdapter.hpp" | 
| 83 | 
+ | 
#include "types/PolarizableAdapter.hpp" | 
| 84 | 
+ | 
#include "types/FixedChargeAdapter.hpp" | 
| 85 | 
+ | 
#include "types/FluctuatingChargeAdapter.hpp" | 
| 86 | 
+ | 
 | 
| 87 | 
  | 
#ifdef IS_MPI | 
| 88 | 
  | 
#include "mpi.h" | 
| 89 | 
  | 
#include "math/ParallelRandNumGen.hpp" | 
| 256 | 
  | 
    int metaDataBlockStart = -1; | 
| 257 | 
  | 
    int metaDataBlockEnd = -1; | 
| 258 | 
  | 
    int i; | 
| 259 | 
< | 
    int mdOffset; | 
| 259 | 
> | 
    streamoff mdOffset; | 
| 260 | 
  | 
    int mdFileVersion; | 
| 261 | 
  | 
 | 
| 262 | 
  | 
#ifdef IS_MPI             | 
| 371 | 
  | 
                                   metaDataBlockStart + 1); | 
| 372 | 
  | 
     | 
| 373 | 
  | 
    //create the force field | 
| 374 | 
< | 
    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField()); | 
| 374 | 
> | 
    ForceField * ff = new ForceField(simParams->getForceField()); | 
| 375 | 
  | 
 | 
| 376 | 
  | 
    if (ff == NULL) { | 
| 377 | 
  | 
      sprintf(painCave.errMsg,  | 
| 422 | 
  | 
    //create the molecules | 
| 423 | 
  | 
    createMolecules(info); | 
| 424 | 
  | 
     | 
| 425 | 
+ | 
    //find the storage layout | 
| 426 | 
+ | 
 | 
| 427 | 
+ | 
    int storageLayout = computeStorageLayout(info); | 
| 428 | 
+ | 
 | 
| 429 | 
  | 
    //allocate memory for DataStorage(circular reference, need to | 
| 430 | 
  | 
    //break it) | 
| 431 | 
< | 
    info->setSnapshotManager(new SimSnapshotManager(info)); | 
| 431 | 
> | 
    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout)); | 
| 432 | 
  | 
     | 
| 433 | 
  | 
    //set the global index of atoms, rigidbodies and cutoffgroups | 
| 434 | 
  | 
    //(only need to be set once, the global index will never change | 
| 665 | 
  | 
    } //end for(int i=0)    | 
| 666 | 
  | 
  } | 
| 667 | 
  | 
     | 
| 668 | 
+ | 
  int SimCreator::computeStorageLayout(SimInfo* info) { | 
| 669 | 
+ | 
 | 
| 670 | 
+ | 
    Globals* simParams = info->getSimParams(); | 
| 671 | 
+ | 
    int nRigidBodies = info->getNGlobalRigidBodies(); | 
| 672 | 
+ | 
    set<AtomType*> atomTypes = info->getSimulatedAtomTypes(); | 
| 673 | 
+ | 
    set<AtomType*>::iterator i; | 
| 674 | 
+ | 
    bool hasDirectionalAtoms = false; | 
| 675 | 
+ | 
    bool hasFixedCharge = false; | 
| 676 | 
+ | 
    bool hasMultipoles = false;     | 
| 677 | 
+ | 
    bool hasPolarizable = false;     | 
| 678 | 
+ | 
    bool hasFluctuatingCharge = false;     | 
| 679 | 
+ | 
    bool hasMetallic = false; | 
| 680 | 
+ | 
    int storageLayout = 0; | 
| 681 | 
+ | 
    storageLayout |= DataStorage::dslPosition; | 
| 682 | 
+ | 
    storageLayout |= DataStorage::dslVelocity; | 
| 683 | 
+ | 
    storageLayout |= DataStorage::dslForce; | 
| 684 | 
+ | 
 | 
| 685 | 
+ | 
    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 686 | 
+ | 
 | 
| 687 | 
+ | 
      DirectionalAdapter da = DirectionalAdapter( (*i) ); | 
| 688 | 
+ | 
      MultipoleAdapter ma = MultipoleAdapter( (*i) ); | 
| 689 | 
+ | 
      EAMAdapter ea = EAMAdapter( (*i) ); | 
| 690 | 
+ | 
      SuttonChenAdapter sca = SuttonChenAdapter( (*i) ); | 
| 691 | 
+ | 
      PolarizableAdapter pa = PolarizableAdapter( (*i) ); | 
| 692 | 
+ | 
      FixedChargeAdapter fca = FixedChargeAdapter( (*i) ); | 
| 693 | 
+ | 
      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) ); | 
| 694 | 
+ | 
 | 
| 695 | 
+ | 
      if (da.isDirectional()){ | 
| 696 | 
+ | 
        hasDirectionalAtoms = true; | 
| 697 | 
+ | 
      } | 
| 698 | 
+ | 
      if (ma.isMultipole()){ | 
| 699 | 
+ | 
        hasMultipoles = true; | 
| 700 | 
+ | 
      } | 
| 701 | 
+ | 
      if (ea.isEAM() || sca.isSuttonChen()){ | 
| 702 | 
+ | 
        hasMetallic = true; | 
| 703 | 
+ | 
      } | 
| 704 | 
+ | 
      if ( fca.isFixedCharge() ){ | 
| 705 | 
+ | 
        hasFixedCharge = true; | 
| 706 | 
+ | 
      } | 
| 707 | 
+ | 
      if ( fqa.isFluctuatingCharge() ){ | 
| 708 | 
+ | 
        hasFluctuatingCharge = true; | 
| 709 | 
+ | 
      } | 
| 710 | 
+ | 
      if ( pa.isPolarizable() ){ | 
| 711 | 
+ | 
        hasPolarizable = true; | 
| 712 | 
+ | 
      } | 
| 713 | 
+ | 
    } | 
| 714 | 
+ | 
     | 
| 715 | 
+ | 
    if (nRigidBodies > 0 || hasDirectionalAtoms) { | 
| 716 | 
+ | 
      storageLayout |= DataStorage::dslAmat; | 
| 717 | 
+ | 
      if(storageLayout & DataStorage::dslVelocity) { | 
| 718 | 
+ | 
        storageLayout |= DataStorage::dslAngularMomentum; | 
| 719 | 
+ | 
      } | 
| 720 | 
+ | 
      if (storageLayout & DataStorage::dslForce) { | 
| 721 | 
+ | 
        storageLayout |= DataStorage::dslTorque; | 
| 722 | 
+ | 
      } | 
| 723 | 
+ | 
    } | 
| 724 | 
+ | 
    if (hasFixedCharge || hasFluctuatingCharge) { | 
| 725 | 
+ | 
      storageLayout |= DataStorage::dslSkippedCharge; | 
| 726 | 
+ | 
    } | 
| 727 | 
+ | 
    if (hasMetallic) { | 
| 728 | 
+ | 
      storageLayout |= DataStorage::dslDensity; | 
| 729 | 
+ | 
      storageLayout |= DataStorage::dslFunctional; | 
| 730 | 
+ | 
      storageLayout |= DataStorage::dslFunctionalDerivative; | 
| 731 | 
+ | 
    } | 
| 732 | 
+ | 
    if (hasPolarizable) { | 
| 733 | 
+ | 
      storageLayout |= DataStorage::dslElectricField; | 
| 734 | 
+ | 
    } | 
| 735 | 
+ | 
    if (hasFluctuatingCharge){ | 
| 736 | 
+ | 
      storageLayout |= DataStorage::dslFlucQPosition; | 
| 737 | 
+ | 
      if(storageLayout & DataStorage::dslVelocity) { | 
| 738 | 
+ | 
        storageLayout |= DataStorage::dslFlucQVelocity; | 
| 739 | 
+ | 
      } | 
| 740 | 
+ | 
      if (storageLayout & DataStorage::dslForce) { | 
| 741 | 
+ | 
        storageLayout |= DataStorage::dslFlucQForce; | 
| 742 | 
+ | 
      } | 
| 743 | 
+ | 
    } | 
| 744 | 
+ | 
     | 
| 745 | 
+ | 
    // if the user has asked for them, make sure we've got the memory for the | 
| 746 | 
+ | 
    // objects defined. | 
| 747 | 
+ | 
 | 
| 748 | 
+ | 
    if (simParams->getOutputParticlePotential()) { | 
| 749 | 
+ | 
      storageLayout |= DataStorage::dslParticlePot; | 
| 750 | 
+ | 
    } | 
| 751 | 
+ | 
 | 
| 752 | 
+ | 
    if (simParams->havePrintHeatFlux()) { | 
| 753 | 
+ | 
      if (simParams->getPrintHeatFlux()) { | 
| 754 | 
+ | 
        storageLayout |= DataStorage::dslParticlePot; | 
| 755 | 
+ | 
      } | 
| 756 | 
+ | 
    } | 
| 757 | 
+ | 
 | 
| 758 | 
+ | 
    if (simParams->getOutputElectricField()) { | 
| 759 | 
+ | 
      storageLayout |= DataStorage::dslElectricField; | 
| 760 | 
+ | 
    } | 
| 761 | 
+ | 
 | 
| 762 | 
+ | 
    if (simParams->getOutputFluctuatingCharges()) { | 
| 763 | 
+ | 
      storageLayout |= DataStorage::dslFlucQPosition; | 
| 764 | 
+ | 
      storageLayout |= DataStorage::dslFlucQVelocity; | 
| 765 | 
+ | 
      storageLayout |= DataStorage::dslFlucQForce; | 
| 766 | 
+ | 
    } | 
| 767 | 
+ | 
 | 
| 768 | 
+ | 
    return storageLayout; | 
| 769 | 
+ | 
  } | 
| 770 | 
+ | 
 | 
| 771 | 
  | 
  void SimCreator::setGlobalIndex(SimInfo *info) { | 
| 772 | 
  | 
    SimInfo::MoleculeIterator mi; | 
| 773 | 
  | 
    Molecule::AtomIterator ai; | 
| 784 | 
  | 
    int nGlobalAtoms = info->getNGlobalAtoms(); | 
| 785 | 
  | 
     | 
| 786 | 
  | 
    beginAtomIndex = 0; | 
| 787 | 
< | 
    beginRigidBodyIndex = 0; | 
| 787 | 
> | 
    //rigidbody's index begins right after atom's | 
| 788 | 
> | 
    beginRigidBodyIndex = info->getNGlobalAtoms(); | 
| 789 | 
  | 
    beginCutoffGroupIndex = 0; | 
| 790 | 
  | 
 | 
| 791 | 
  | 
    for(int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 903 | 
  | 
    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 904 | 
  | 
      int myGlobalIndex = mol->getGlobalIndex(); | 
| 905 | 
  | 
      int globalIO = startingIOIndexForMol[myGlobalIndex]; | 
| 906 | 
< | 
      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL; | 
| 907 | 
< | 
           integrableObject = mol->nextIntegrableObject(ioi)) { | 
| 908 | 
< | 
        integrableObject->setGlobalIntegrableObjectIndex(globalIO); | 
| 909 | 
< | 
        IOIndexToIntegrableObject[globalIO] = integrableObject; | 
| 906 | 
> | 
      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; | 
| 907 | 
> | 
           sd = mol->nextIntegrableObject(ioi)) { | 
| 908 | 
> | 
        sd->setGlobalIntegrableObjectIndex(globalIO); | 
| 909 | 
> | 
        IOIndexToIntegrableObject[globalIO] = sd; | 
| 910 | 
  | 
        globalIO++; | 
| 911 | 
  | 
      } | 
| 912 | 
  | 
    } |