ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/SimCreator.cpp
(Generate patch)

Comparing branches/development/src/brains/SimCreator.cpp (file contents):
Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC vs.
Revision 1798 by gezelter, Thu Sep 13 14:10:11 2012 UTC

# Line 56 | Line 56
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"
# Line 76 | Line 76
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"
# Line 248 | Line 256 | namespace OpenMD {
256      int metaDataBlockStart = -1;
257      int metaDataBlockEnd = -1;
258      int i;
259 <    int mdOffset;
259 >    streamoff mdOffset;
260      int mdFileVersion;
261  
262 +
263   #ifdef IS_MPI            
264      const int masterNode = 0;
265      if (worldRank == masterNode) {
266   #endif
267  
268 <      std::ifstream mdFile_(mdFileName.c_str());
268 >      std::ifstream mdFile_;
269 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
270        
271        if (mdFile_.fail()) {
272          sprintf(painCave.errMsg,
# Line 363 | Line 373 | namespace OpenMD {
373                                     metaDataBlockStart + 1);
374      
375      //create the force field
376 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
376 >    ForceField * ff = new ForceField(simParams->getForceField());
377  
378      if (ff == NULL) {
379        sprintf(painCave.errMsg,
# Line 414 | Line 424 | namespace OpenMD {
424      //create the molecules
425      createMolecules(info);
426      
427 +    //find the storage layout
428 +
429 +    int storageLayout = computeStorageLayout(info);
430 +
431      //allocate memory for DataStorage(circular reference, need to
432      //break it)
433 <    info->setSnapshotManager(new SimSnapshotManager(info));
433 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
434      
435      //set the global index of atoms, rigidbodies and cutoffgroups
436      //(only need to be set once, the global index will never change
# Line 493 | Line 507 | namespace OpenMD {
507      int nGlobalMols = info->getNGlobalMolecules();
508      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
509      
510 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
510 >    nProcessors = MPI::COMM_WORLD.Get_size();
511      
512      if (nProcessors > nGlobalMols) {
513        sprintf(painCave.errMsg,
# Line 608 | Line 622 | namespace OpenMD {
622        delete myRandom;
623        
624        // Spray out this nonsense to all other processors:
625 <      
612 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
625 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
626      } else {
627        
628        // Listen to your marching orders from processor 0:
629 <      
617 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
629 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
630      }
631      
632      info->setMolToProcMap(molToProcMap);
# Line 653 | Line 665 | namespace OpenMD {
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 hasDipoles = false;    
677 +    bool hasQuadrupoles = false;    
678 +    bool hasPolarizable = false;    
679 +    bool hasFluctuatingCharge = false;    
680 +    bool hasMetallic = false;
681 +    int storageLayout = 0;
682 +    storageLayout |= DataStorage::dslPosition;
683 +    storageLayout |= DataStorage::dslVelocity;
684 +    storageLayout |= DataStorage::dslForce;
685 +
686 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
687 +
688 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
689 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
690 +      EAMAdapter ea = EAMAdapter( (*i) );
691 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
692 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
693 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
694 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
695 +
696 +      if (da.isDirectional()){
697 +        hasDirectionalAtoms = true;
698 +      }
699 +      if (ma.isDipole()){
700 +        hasDipoles = true;
701 +      }
702 +      if (ma.isQuadrupole()){
703 +        hasQuadrupoles = true;
704 +      }
705 +      if (ea.isEAM() || sca.isSuttonChen()){
706 +        hasMetallic = true;
707 +      }
708 +      if ( fca.isFixedCharge() ){
709 +        hasFixedCharge = true;
710 +      }
711 +      if ( fqa.isFluctuatingCharge() ){
712 +        hasFluctuatingCharge = true;
713 +      }
714 +      if ( pa.isPolarizable() ){
715 +        hasPolarizable = true;
716 +      }
717 +    }
718 +    
719 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
720 +      storageLayout |= DataStorage::dslAmat;
721 +      if(storageLayout & DataStorage::dslVelocity) {
722 +        storageLayout |= DataStorage::dslAngularMomentum;
723 +      }
724 +      if (storageLayout & DataStorage::dslForce) {
725 +        storageLayout |= DataStorage::dslTorque;
726 +      }
727 +    }
728 +    if (hasDipoles) {
729 +      storageLayout |= DataStorage::dslDipole;
730 +    }
731 +    if (hasQuadrupoles) {
732 +      storageLayout |= DataStorage::dslQuadrupole;
733 +    }
734 +    if (hasFixedCharge || hasFluctuatingCharge) {
735 +      storageLayout |= DataStorage::dslSkippedCharge;
736 +    }
737 +    if (hasMetallic) {
738 +      storageLayout |= DataStorage::dslDensity;
739 +      storageLayout |= DataStorage::dslFunctional;
740 +      storageLayout |= DataStorage::dslFunctionalDerivative;
741 +    }
742 +    if (hasPolarizable) {
743 +      storageLayout |= DataStorage::dslElectricField;
744 +    }
745 +    if (hasFluctuatingCharge){
746 +      storageLayout |= DataStorage::dslFlucQPosition;
747 +      if(storageLayout & DataStorage::dslVelocity) {
748 +        storageLayout |= DataStorage::dslFlucQVelocity;
749 +      }
750 +      if (storageLayout & DataStorage::dslForce) {
751 +        storageLayout |= DataStorage::dslFlucQForce;
752 +      }
753 +    }
754 +    
755 +    // if the user has asked for them, make sure we've got the memory for the
756 +    // objects defined.
757 +
758 +    if (simParams->getOutputParticlePotential()) {
759 +      storageLayout |= DataStorage::dslParticlePot;
760 +    }
761 +
762 +    if (simParams->havePrintHeatFlux()) {
763 +      if (simParams->getPrintHeatFlux()) {
764 +        storageLayout |= DataStorage::dslParticlePot;
765 +      }
766 +    }
767 +
768 +    if (simParams->getOutputElectricField()) {
769 +      storageLayout |= DataStorage::dslElectricField;
770 +    }
771 +
772 +    if (simParams->getOutputFluctuatingCharges()) {
773 +      storageLayout |= DataStorage::dslFlucQPosition;
774 +      storageLayout |= DataStorage::dslFlucQVelocity;
775 +      storageLayout |= DataStorage::dslFlucQForce;
776 +    }
777 +
778 +    return storageLayout;
779 +  }
780 +
781    void SimCreator::setGlobalIndex(SimInfo *info) {
782      SimInfo::MoleculeIterator mi;
783      Molecule::AtomIterator ai;
# Line 669 | Line 794 | namespace OpenMD {
794      int nGlobalAtoms = info->getNGlobalAtoms();
795      
796      beginAtomIndex = 0;
797 <    beginRigidBodyIndex = 0;
797 >    //rigidbody's index begins right after atom's
798 >    beginRigidBodyIndex = info->getNGlobalAtoms();
799      beginCutoffGroupIndex = 0;
800  
801      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
# Line 732 | Line 858 | namespace OpenMD {
858      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
859      // docs said we could.
860      std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
861 <    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
862 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
861 >    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
862 >                              &tmpGroupMembership[0], nGlobalAtoms,
863 >                              MPI::INT, MPI::SUM);
864      info->setGlobalGroupMembership(tmpGroupMembership);
865   #else
866      info->setGlobalGroupMembership(globalGroupMembership);
# Line 750 | Line 877 | namespace OpenMD {
877      
878   #ifdef IS_MPI
879      std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
880 +    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
881 +                              nGlobalAtoms,
882 +                              MPI::INT, MPI::SUM);
883      
754 –    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
755 –                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
756 –    
884      info->setGlobalMolMembership(tmpMolMembership);
885   #else
886      info->setGlobalMolMembership(globalMolMembership);
# Line 769 | Line 896 | namespace OpenMD {
896      
897   #ifdef IS_MPI
898      std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
899 <    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
900 <                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
899 >    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
900 >                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
901   #else
902      std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
903   #endif    
# Line 787 | Line 914 | namespace OpenMD {
914      for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
915        int myGlobalIndex = mol->getGlobalIndex();
916        int globalIO = startingIOIndexForMol[myGlobalIndex];
917 <      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
918 <           integrableObject = mol->nextIntegrableObject(ioi)) {
919 <        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
920 <        IOIndexToIntegrableObject[globalIO] = integrableObject;
917 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
918 >           sd = mol->nextIntegrableObject(ioi)) {
919 >        sd->setGlobalIntegrableObjectIndex(globalIO);
920 >        IOIndexToIntegrableObject[globalIO] = sd;
921          globalIO++;
922        }
923      }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines