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

Comparing trunk/src/brains/SimCreator.cpp (file contents):
Revision 1810 by gezelter, Thu Nov 8 14:23:43 2012 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

# Line 1 | Line 1
1   /*
2 < * copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
5   * non-exclusive, royalty free, license to use, modify and
# Line 35 | Line 35
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 44 | Line 44
44   * @file SimCreator.cpp
45   * @author tlin
46   * @date 11/03/2004
47 – * @time 13:51am
47   * @version 1.0
48   */
49 +
50 + #ifdef IS_MPI
51 + #include "mpi.h"
52 + #include "math/ParallelRandNumGen.hpp"
53 + #endif
54 +
55   #include <exception>
56   #include <iostream>
57   #include <sstream>
# Line 84 | Line 89
89   #include "types/FixedChargeAdapter.hpp"
90   #include "types/FluctuatingChargeAdapter.hpp"
91  
87 – #ifdef IS_MPI
88 – #include "mpi.h"
89 – #include "math/ParallelRandNumGen.hpp"
90 – #endif
92  
93   namespace OpenMD {
94    
# Line 105 | Line 106 | namespace OpenMD {
106          MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
107   #endif                
108          SimplePreprocessor preprocessor;
109 <        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock,
109 <                                ppStream);
109 >        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
110                  
111   #ifdef IS_MPI            
112          //brocasting the stream size
113          streamSize = ppStream.str().size() +1;
114          MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
115 <        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
116 <                              streamSize, MPI::CHAR, masterNode);
117 <                
115 >        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI::CHAR, masterNode);
116 >                          
117        } else {
119 –
118          MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
119  
120          //get stream size
# Line 130 | Line 128 | namespace OpenMD {
128                  
129          ppStream.str(buf);
130          delete [] buf;
133 –
131        }
132   #endif            
133        // Create a scanner that reads from the input stream
# Line 256 | Line 253 | namespace OpenMD {
253      int metaDataBlockStart = -1;
254      int metaDataBlockEnd = -1;
255      int i, j;
256 <    streamoff mdOffset(0);
256 >    streamoff mdOffset;
257      int mdFileVersion;
258  
259      // Create a string for embedding the version information in the MetaData
# Line 515 | Line 512 | namespace OpenMD {
512    
513   #ifdef IS_MPI
514    void SimCreator::divideMolecules(SimInfo *info) {
518 –    RealType numerator;
519 –    RealType denominator;
520 –    RealType precast;
521 –    RealType x;
522 –    RealType y;
515      RealType a;
524 –    int old_atoms;
525 –    int add_atoms;
526 –    int new_atoms;
527 –    int nTarget;
528 –    int done;
529 –    int i;
530 –    int loops;
531 –    int which_proc;
516      int nProcessors;
517      std::vector<int> atomsPerProc;
518      int nGlobalMols = info->getNGlobalMolecules();
519 <    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
519 >    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an
520 >                                                    // error
521 >                                                    // condition:
522      
523      nProcessors = MPI::COMM_WORLD.Get_size();
524      
# Line 543 | Line 529 | namespace OpenMD {
529                "\tthe number of molecules.  This will not result in a \n"
530                "\tusable division of atoms for force decomposition.\n"
531                "\tEither try a smaller number of processors, or run the\n"
532 <              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
532 >              "\tsingle-processor version of OpenMD.\n", nProcessors,
533 >              nGlobalMols);
534        
535        painCave.isFatal = 1;
536        simError();
537      }
538      
552 –    int seedValue;
539      Globals * simParams = info->getSimParams();
540 <    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
540 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel
541 >                             //random number generator
542      if (simParams->haveSeed()) {
543 <      seedValue = simParams->getSeed();
543 >      int seedValue = simParams->getSeed();
544        myRandom = new SeqRandNumGen(seedValue);
545      }else {
546        myRandom = new SeqRandNumGen();
# Line 566 | Line 553 | namespace OpenMD {
553      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
554      
555      if (worldRank == 0) {
556 <      numerator = info->getNGlobalAtoms();
557 <      denominator = nProcessors;
558 <      precast = numerator / denominator;
559 <      nTarget = (int)(precast + 0.5);
556 >      RealType numerator = info->getNGlobalAtoms();
557 >      RealType denominator = nProcessors;
558 >      RealType precast = numerator / denominator;
559 >      int nTarget = (int)(precast + 0.5);
560        
561 <      for(i = 0; i < nGlobalMols; i++) {
561 >      for(int i = 0; i < nGlobalMols; i++) {
562  
563 <        done = 0;
564 <        loops = 0;
563 >        int done = 0;
564 >        int loops = 0;
565          
566          while (!done) {
567            loops++;
568            
569            // Pick a processor at random
570            
571 <          which_proc = (int) (myRandom->rand() * nProcessors);
571 >          int which_proc = (int) (myRandom->rand() * nProcessors);
572            
573            //get the molecule stamp first
574            int stampId = info->getMoleculeStampId(i);
575            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
576            
577            // How many atoms does this processor have so far?
578 <          old_atoms = atomsPerProc[which_proc];
579 <          add_atoms = moleculeStamp->getNAtoms();
580 <          new_atoms = old_atoms + add_atoms;
578 >          int old_atoms = atomsPerProc[which_proc];
579 >          int add_atoms = moleculeStamp->getNAtoms();
580 >          int new_atoms = old_atoms + add_atoms;
581            
582            // If we've been through this loop too many times, we need
583            // to just give up and assign the molecule to this processor
# Line 634 | Line 621 | namespace OpenMD {
621            //           Pacc(x) = exp(- a * x)
622            // where a = penalty / (average atoms per molecule)
623            
624 <          x = (RealType)(new_atoms - nTarget);
625 <          y = myRandom->rand();
624 >          RealType x = (RealType)(new_atoms - nTarget);
625 >          RealType y = myRandom->rand();
626            
627            if (y < exp(- a * x)) {
628              molToProcMap[i] = which_proc;
# Line 704 | Line 691 | namespace OpenMD {
691      set<AtomType*>::iterator i;
692      bool hasDirectionalAtoms = false;
693      bool hasFixedCharge = false;
694 <    bool hasMultipoles = false;    
694 >    bool hasDipoles = false;    
695 >    bool hasQuadrupoles = false;    
696      bool hasPolarizable = false;    
697      bool hasFluctuatingCharge = false;    
698      bool hasMetallic = false;
# Line 726 | Line 714 | namespace OpenMD {
714        if (da.isDirectional()){
715          hasDirectionalAtoms = true;
716        }
717 <      if (ma.isMultipole()){
718 <        hasMultipoles = true;
717 >      if (ma.isDipole()){
718 >        hasDipoles = true;
719        }
720 +      if (ma.isQuadrupole()){
721 +        hasQuadrupoles = true;
722 +      }
723        if (ea.isEAM() || sca.isSuttonChen()){
724          hasMetallic = true;
725        }
# Line 752 | Line 743 | namespace OpenMD {
743          storageLayout |= DataStorage::dslTorque;
744        }
745      }
746 <    if (hasMultipoles) {
747 <      storageLayout |= DataStorage::dslElectroFrame;
746 >    if (hasDipoles) {
747 >      storageLayout |= DataStorage::dslDipole;
748      }
749 +    if (hasQuadrupoles) {
750 +      storageLayout |= DataStorage::dslQuadrupole;
751 +    }
752      if (hasFixedCharge || hasFluctuatingCharge) {
753        storageLayout |= DataStorage::dslSkippedCharge;
754      }
# Line 789 | Line 783 | namespace OpenMD {
783        }
784      }
785  
786 <    if (simParams->getOutputElectricField()) {
786 >    if (simParams->getOutputElectricField() | simParams->haveElectricField()) {
787        storageLayout |= DataStorage::dslElectricField;
788      }
789 +
790      if (simParams->getOutputFluctuatingCharges()) {
791        storageLayout |= DataStorage::dslFlucQPosition;
792        storageLayout |= DataStorage::dslFlucQVelocity;
793        storageLayout |= DataStorage::dslFlucQForce;
794      }
795  
796 +    info->setStorageLayout(storageLayout);
797 +
798      return storageLayout;
799    }
800  
# Line 806 | Line 803 | namespace OpenMD {
803      Molecule::AtomIterator ai;
804      Molecule::RigidBodyIterator ri;
805      Molecule::CutoffGroupIterator ci;
806 +    Molecule::BondIterator boi;
807 +    Molecule::BendIterator bei;
808 +    Molecule::TorsionIterator ti;
809 +    Molecule::InversionIterator ii;
810      Molecule::IntegrableObjectIterator  ioi;
811 <    Molecule * mol;
812 <    Atom * atom;
813 <    RigidBody * rb;
814 <    CutoffGroup * cg;
811 >    Molecule* mol;
812 >    Atom* atom;
813 >    RigidBody* rb;
814 >    CutoffGroup* cg;
815 >    Bond* bond;
816 >    Bend* bend;
817 >    Torsion* torsion;
818 >    Inversion* inversion;
819      int beginAtomIndex;
820      int beginRigidBodyIndex;
821      int beginCutoffGroupIndex;
822 +    int beginBondIndex;
823 +    int beginBendIndex;
824 +    int beginTorsionIndex;
825 +    int beginInversionIndex;
826      int nGlobalAtoms = info->getNGlobalAtoms();
827      int nGlobalRigidBodies = info->getNGlobalRigidBodies();
828      
829      beginAtomIndex = 0;
830 <    //rigidbody's index begins right after atom's
830 >    // The rigid body indices begin immediately after the atom indices:
831      beginRigidBodyIndex = info->getNGlobalAtoms();
832      beginCutoffGroupIndex = 0;
833 <
833 >    beginBondIndex = 0;
834 >    beginBendIndex = 0;
835 >    beginTorsionIndex = 0;
836 >    beginInversionIndex = 0;
837 >  
838      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
839        
840   #ifdef IS_MPI      
# Line 830 | Line 843 | namespace OpenMD {
843          // stuff to do if I own this molecule
844          mol = info->getMoleculeByGlobalIndex(i);
845  
846 <        //local index(index in DataStorge) of atom is important
847 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
846 >        // The local index(index in DataStorge) of the atom is important:
847 >        for(atom = mol->beginAtom(ai); atom != NULL;
848 >            atom = mol->nextAtom(ai)) {
849            atom->setGlobalIndex(beginAtomIndex++);
850          }
851          
# Line 840 | Line 854 | namespace OpenMD {
854            rb->setGlobalIndex(beginRigidBodyIndex++);
855          }
856          
857 <        //local index of cutoff group is trivial, it only depends on
858 <        //the order of travesing
857 >        // The local index of other objects only depends on the order
858 >        // of traversal:
859          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
860              cg = mol->nextCutoffGroup(ci)) {
861            cg->setGlobalIndex(beginCutoffGroupIndex++);
862          }        
863 +        for(bond = mol->beginBond(boi); bond != NULL;
864 +            bond = mol->nextBond(boi)) {
865 +          bond->setGlobalIndex(beginBondIndex++);
866 +        }        
867 +        for(bend = mol->beginBend(bei); bend != NULL;
868 +            bend = mol->nextBend(bei)) {
869 +          bend->setGlobalIndex(beginBendIndex++);
870 +        }        
871 +        for(torsion = mol->beginTorsion(ti); torsion != NULL;
872 +            torsion = mol->nextTorsion(ti)) {
873 +          torsion->setGlobalIndex(beginTorsionIndex++);
874 +        }        
875 +        for(inversion = mol->beginInversion(ii); inversion != NULL;
876 +            inversion = mol->nextInversion(ii)) {
877 +          inversion->setGlobalIndex(beginInversionIndex++);
878 +        }        
879          
880   #ifdef IS_MPI        
881        }  else {
# Line 858 | Line 888 | namespace OpenMD {
888          beginAtomIndex += stamp->getNAtoms();
889          beginRigidBodyIndex += stamp->getNRigidBodies();
890          beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
891 +        beginBondIndex += stamp->getNBonds();
892 +        beginBendIndex += stamp->getNBends();
893 +        beginTorsionIndex += stamp->getNTorsions();
894 +        beginInversionIndex += stamp->getNInversions();
895        }
896   #endif          
897  
# Line 865 | Line 899 | namespace OpenMD {
899  
900      //fill globalGroupMembership
901      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
902 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
903 <      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
904 <        
902 >    for(mol = info->beginMolecule(mi); mol != NULL;
903 >        mol = info->nextMolecule(mi)) {        
904 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
905 >           cg = mol->nextCutoffGroup(ci)) {        
906          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
907            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
908          }
# Line 960 | Line 995 | namespace OpenMD {
995    }
996    
997    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
998 <
998 >    
999      DumpReader reader(info, mdFileName);
1000      int nframes = reader.getNFrames();
1001 <
1001 >    
1002      if (nframes > 0) {
1003        reader.readFrame(nframes - 1);
1004      } else {

Diff Legend

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