src/brains/SimCreator.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
/**
 * @file SimCreator.cpp
 * @author tlin
 * @date 11/03/2004
 * @time 13:51am
 * @version 1.0
 */

#include "brains/MoleculeCreator.hpp"
#include "brains/SimCreator.hpp"
#include "brains/SimSnapshotManager.hpp"
#include "io/DumpReader.hpp"
#include "io/parse_me.h"
#include "UseTheForce/ForceFieldFactory.hpp"
#include "utils/simError.h"
#include "utils/StringUtils.hpp"
#include "math/SeqRandNumGen.hpp"
#ifdef IS_MPI
#include "io/mpiBASS.h"
#include "math/ParallelRandNumGen.hpp"
#endif

namespace oopse {
  
  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, 
                             Globals* simParams){
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif // is_mpi
      
      set_interface_stamps(stamps, simParams);
      
#ifdef IS_MPI
      
      mpiEventInit();
      
#endif
      
      yacc_BASS(mdFileName.c_str());
      
#ifdef IS_MPI
      
      throwMPIEvent(NULL);
    } else {
      set_interface_stamps(stamps, simParams);
      mpiEventInit();
      MPIcheckPoint();
      mpiEventLoop();
    }
    
#endif
    
  }
  
  SimInfo*  SimCreator::createSim(const std::string & mdFileName, 
                                  bool loadInitCoords) {
    
    MakeStamps * stamps = new MakeStamps();
    
    Globals * simParams = new Globals();
    
    //parse meta-data file
    parseFile(mdFileName, stamps, simParams);
    
    //create the force field
    ForceField * ff = ForceFieldFactory::getInstance()
      ->createForceField(simParams->getForceField());
    
    if (ff == NULL) {
      sprintf(painCave.errMsg, 
              "ForceField Factory can not create %s force field\n",
              simParams->getForceField().c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    if (simParams->haveForceFieldFileName()) {
      ff->setForceFieldFileName(simParams->getForceFieldFileName());
    }
    
    std::string forcefieldFileName;
    forcefieldFileName = ff->getForceFieldFileName();
    
    if (simParams->haveForceFieldVariant()) {
      //If the force field has variant, the variant force field name will be
      //Base.variant.frc. For exampel EAM.u6.frc
      
      std::string variant = simParams->getForceFieldVariant();
      
      std::string::size_type pos = forcefieldFileName.rfind(".frc");
      variant = "." + variant;
      if (pos != std::string::npos) {
        forcefieldFileName.insert(pos, variant);
      } else {
        //If the default force field file name does not containt .frc suffix, just append the .variant
        forcefieldFileName.append(variant);
      }
    } 
    
    ff->parse(forcefieldFileName);
    
    //extract the molecule stamps
    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
    compList(stamps, simParams, moleculeStampPairs);
    
    //create SimInfo
    SimInfo * info = new SimInfo(stamps, moleculeStampPairs, ff, simParams);
     
    //gather parameters (SimCreator only retrieves part of the parameters)
    gatherParameters(info, mdFileName);
    
    //divide the molecules and determine the global index of molecules
#ifdef IS_MPI
    divideMolecules(info);
#endif 
    
    //create the molecules
    createMolecules(info);
    
    
    //allocate memory for DataStorage(circular reference, need to break it)
    info->setSnapshotManager(new SimSnapshotManager(info));
    
    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
    //global index will never change again). Local indices of atoms and rigidbodies are already set by 
    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager. 
    setGlobalIndex(info);
    
    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
    //atoms don't have the global index yet  (their global index are all initialized to -1).
    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
    //we can determine the beginning global indices of atoms before they get created.
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
      info->addExcludePairs(mol);
    }
    
    if (loadInitCoords)
      loadCoordinates(info);    
    
    return info;
  }
  
  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
    
    //figure out the output file names
    std::string prefix;
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif // is_mpi
      Globals * simParams = info->getSimParams();
      if (simParams->haveFinalConfig()) {
        prefix = getPrefix(simParams->getFinalConfig());
      } else {
        prefix = getPrefix(mdfile);
      }
      
      info->setFinalConfigFileName(prefix + ".eor");
      info->setDumpFileName(prefix + ".dump");
      info->setStatFileName(prefix + ".stat");
      info->setRestFileName(prefix + ".zang");
      
#ifdef IS_MPI
      
    }
    
#endif
    
  }
  
#ifdef IS_MPI
  void SimCreator::divideMolecules(SimInfo *info) {
    double numerator;
    double denominator;
    double precast;
    double x;
    double y;
    double a;
    int old_atoms;
    int add_atoms;
    int new_atoms;
    int nTarget;
    int done;
    int i;
    int j;
    int loops;
    int which_proc;
    int nProcessors;
    std::vector<int> atomsPerProc;
    int nGlobalMols = info->getNGlobalMolecules();
    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
    
    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
    
    if (nProcessors > nGlobalMols) {
      sprintf(painCave.errMsg,
              "nProcessors (%d) > nMol (%d)\n"
              "\tThe number of processors is larger than\n"
              "\tthe number of molecules.  This will not result in a \n"
              "\tusable division of atoms for force decomposition.\n"
              "\tEither try a smaller number of processors, or run the\n"
              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
      
      painCave.isFatal = 1;
      simError();
    }
    
    int seedValue;
    Globals * simParams = info->getSimParams();
    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
    if (simParams->haveSeed()) {
      seedValue = simParams->getSeed();
      myRandom = new SeqRandNumGen(seedValue);
    }else {
      myRandom = new SeqRandNumGen();
    }   
    
    
    a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
    
    //initialize atomsPerProc
    atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
    
    if (worldRank == 0) {
      numerator = info->getNGlobalAtoms();
      denominator = nProcessors;
      precast = numerator / denominator;
      nTarget = (int)(precast + 0.5);
      
      for(i = 0; i < nGlobalMols; i++) {
        done = 0;
        loops = 0;
        
        while (!done) {
          loops++;
          
          // Pick a processor at random
          
          which_proc = (int) (myRandom->rand() * nProcessors);
          
          //get the molecule stamp first
          int stampId = info->getMoleculeStampId(i);
          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
          
          // How many atoms does this processor have so far?
          old_atoms = atomsPerProc[which_proc];
          add_atoms = moleculeStamp->getNAtoms();
          new_atoms = old_atoms + add_atoms;
          
          // If we've been through this loop too many times, we need
          // to just give up and assign the molecule to this processor
          // and be done with it. 
          
          if (loops > 100) {
            sprintf(painCave.errMsg,
                    "I've tried 100 times to assign molecule %d to a "
                    " processor, but can't find a good spot.\n"
                    "I'm assigning it at random to processor %d.\n",
                    i, which_proc);
            
            painCave.isFatal = 0;
            simError();
            
            molToProcMap[i] = which_proc;
            atomsPerProc[which_proc] += add_atoms;
            
            done = 1;
            continue;
          }
          
          // If we can add this molecule to this processor without sending
          // it above nTarget, then go ahead and do it:
          
          if (new_atoms <= nTarget) {
            molToProcMap[i] = which_proc;
            atomsPerProc[which_proc] += add_atoms;
            
            done = 1;
            continue;
          }
          
          // The only situation left is when new_atoms > nTarget.  We
          // want to accept this with some probability that dies off the
          // farther we are from nTarget
          
          // roughly:  x = new_atoms - nTarget
          //           Pacc(x) = exp(- a * x)
          // where a = penalty / (average atoms per molecule)
          
          x = (double)(new_atoms - nTarget);
          y = myRandom->rand();
          
          if (y < exp(- a * x)) {
            molToProcMap[i] = which_proc;
            atomsPerProc[which_proc] += add_atoms;
            
            done = 1;
            continue;
          } else {
            continue;
          }
        }
      }
      
      delete myRandom;
      
      // Spray out this nonsense to all other processors:
      
      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    } else {
      
      // Listen to your marching orders from processor 0:
      
      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    }
    
    info->setMolToProcMap(molToProcMap);
    sprintf(checkPointMsg,
            "Successfully divided the molecules among the processors.\n");
    MPIcheckPoint();
  }
  
#endif
  
  void SimCreator::createMolecules(SimInfo *info) {
    MoleculeCreator molCreator;
    int stampId;
    
    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
      
#ifdef IS_MPI
      
      if (info->getMolToProc(i) == worldRank) {
#endif
        
        stampId = info->getMoleculeStampId(i);
        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
                                                   stampId, i, info->getLocalIndexManager());
        
        info->addMolecule(mol);
        
#ifdef IS_MPI
        
      }
      
#endif
      
    } //end for(int i=0)   
  }
  
  void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
    int i;
    char * id;
    MoleculeStamp * currentStamp;
    Component** the_components = simParams->getComponents();
    int n_components = simParams->getNComponents();
    
    if (!simParams->haveNMol()) {
      // we don't have the total number of molecules, so we assume it is
      // given in each component
      
      for(i = 0; i < n_components; i++) {
        if (!the_components[i]->haveNMol()) {
          // we have a problem
          sprintf(painCave.errMsg,
                  "SimCreator Error. No global NMol or component NMol given.\n"
                  "\tCannot calculate the number of atoms.\n");
          
          painCave.isFatal = 1;
          simError();
        }
        
        id = the_components[i]->getType();

        currentStamp = stamps->getMolStamp(id);
        if (currentStamp == NULL) {
          sprintf(painCave.errMsg,
                  "SimCreator error: Component \"%s\" was not found in the "
                  "list of declared molecules\n", id);
          
          painCave.isFatal = 1;
          simError();
        }
        
        moleculeStampPairs.push_back(
                                     std::make_pair(currentStamp, the_components[i]->getNMol()));
      } //end for (i = 0; i < n_components; i++)
    } else {
      sprintf(painCave.errMsg, "SimSetup error.\n"
              "\tSorry, the ability to specify total"
              " nMols and then give molfractions in the components\n"
              "\tis not currently supported."
              " Please give nMol in the components.\n");
      
      painCave.isFatal = 1;
      simError();
    }
    
#ifdef IS_MPI
    
    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
    MPIcheckPoint();
    
#endif // is_mpi
    
  }
  
  void SimCreator::setGlobalIndex(SimInfo *info) {
    SimInfo::MoleculeIterator mi;
    Molecule::AtomIterator ai;
    Molecule::RigidBodyIterator ri;
    Molecule::CutoffGroupIterator ci;
    Molecule * mol;
    Atom * atom;
    RigidBody * rb;
    CutoffGroup * cg;
    int beginAtomIndex;
    int beginRigidBodyIndex;
    int beginCutoffGroupIndex;
    int nGlobalAtoms = info->getNGlobalAtoms();
    
#ifndef IS_MPI
    
    beginAtomIndex = 0;
    beginRigidBodyIndex = 0;
    beginCutoffGroupIndex = 0;
    
#else
    
    int nproc;
    int myNode;
    
    myNode = worldRank;
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
    
    std::vector < int > tmpAtomsInProc(nproc, 0);
    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
    std::vector < int > NumAtomsInProc(nproc, 0);
    std::vector < int > NumRigidBodiesInProc(nproc, 0);
    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
    
    tmpAtomsInProc[myNode] = info->getNAtoms();
    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
    
    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
                  MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    
    beginAtomIndex = 0;
    beginRigidBodyIndex = 0;
    beginCutoffGroupIndex = 0;
    
    for(int i = 0; i < myNode; i++) {
      beginAtomIndex += NumAtomsInProc[i];
      beginRigidBodyIndex += NumRigidBodiesInProc[i];
      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
    }
    
#endif
    
    //rigidbody's index begins right after atom's
    beginRigidBodyIndex += info->getNGlobalAtoms();
    
    for(mol = info->beginMolecule(mi); mol != NULL;
        mol = info->nextMolecule(mi)) {
      
      //local index(index in DataStorge) of atom is important
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
        atom->setGlobalIndex(beginAtomIndex++);
      }
      
      for(rb = mol->beginRigidBody(ri); rb != NULL;
          rb = mol->nextRigidBody(ri)) {
        rb->setGlobalIndex(beginRigidBodyIndex++);
      }
      
      //local index of cutoff group is trivial, it only depends on the order of travesing
      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
          cg = mol->nextCutoffGroup(ci)) {
        cg->setGlobalIndex(beginCutoffGroupIndex++);
      }
    }
    
    //fill globalGroupMembership
    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
        
        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
        }
        
      }       
    }
    
#ifdef IS_MPI    
    // Since the globalGroupMembership has been zero filled and we've only
    // poked values into the atoms we know, we can do an Allreduce
    // to get the full globalGroupMembership array (We think).
    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
    // docs said we could.
    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    info->setGlobalGroupMembership(tmpGroupMembership);
#else
    info->setGlobalGroupMembership(globalGroupMembership);
#endif
    
    //fill molMembership
    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
    
    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
      
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
      }
    }
    
#ifdef IS_MPI
    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
    
    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    
    info->setGlobalMolMembership(tmpMolMembership);
#else
    info->setGlobalMolMembership(globalMolMembership);
#endif
    
  }
  
  void SimCreator::loadCoordinates(SimInfo* info) {
    Globals* simParams;
    simParams = info->getSimParams();
    
    if (!simParams->haveInitialConfig()) {
      sprintf(painCave.errMsg,
              "Cannot intialize a simulation without an initial configuration file.\n");
      painCave.isFatal = 1;;
      simError();
    }
    
    DumpReader reader(info, simParams->getInitialConfig());
    int nframes = reader.getNFrames();
    
    if (nframes > 0) {
      reader.readFrame(nframes - 1);
    } else {
      //invalid initial coordinate file
      sprintf(painCave.errMsg, 
              "Initial configuration file %s should at least contain one frame\n",
              simParams->getInitialConfig().c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    //copy the current snapshot to previous snapshot
    info->getSnapshotManager()->advance();
  }
  
} //end namespace oopse


Revision:	2448
Committed:	Wed Nov 16 23:10:02 2005 UTC (18 years, 7 months ago) by tim
File size:	19974 byte(s)
Log Message:	OptionSectionParser get compiled
#	Content
1	/*
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
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
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41
42	/**
43	* @file SimCreator.cpp
44	* @author tlin
45	* @date 11/03/2004
46	* @time 13:51am
47	* @version 1.0
48	*/
49
50	#include "brains/MoleculeCreator.hpp"
51	#include "brains/SimCreator.hpp"
52	#include "brains/SimSnapshotManager.hpp"
53	#include "io/DumpReader.hpp"
54	#include "io/parse_me.h"
55	#include "UseTheForce/ForceFieldFactory.hpp"
56	#include "utils/simError.h"
57	#include "utils/StringUtils.hpp"
58	#include "math/SeqRandNumGen.hpp"
59	#ifdef IS_MPI
60	#include "io/mpiBASS.h"
61	#include "math/ParallelRandNumGen.hpp"
62	#endif
63
64	namespace oopse {
65
66	void SimCreator::parseFile(const std::string mdFileName, MakeStamps* stamps,
67	Globals* simParams){
68
69	#ifdef IS_MPI
70
71	if (worldRank == 0) {
72	#endif // is_mpi
73
74	set_interface_stamps(stamps, simParams);
75
76	#ifdef IS_MPI
77
78	mpiEventInit();
79
80	#endif
81
82	yacc_BASS(mdFileName.c_str());
83
84	#ifdef IS_MPI
85
86	throwMPIEvent(NULL);
87	} else {
88	set_interface_stamps(stamps, simParams);
89	mpiEventInit();
90	MPIcheckPoint();
91	mpiEventLoop();
92	}
93
94	#endif
95
96	}
97
98	SimInfo* SimCreator::createSim(const std::string & mdFileName,
99	bool loadInitCoords) {
100
101	MakeStamps * stamps = new MakeStamps();
102
103	Globals * simParams = new Globals();
104
105	//parse meta-data file
106	parseFile(mdFileName, stamps, simParams);
107
108	//create the force field
109	ForceField * ff = ForceFieldFactory::getInstance()
110	->createForceField(simParams->getForceField());
111
112	if (ff == NULL) {
113	sprintf(painCave.errMsg,
114	"ForceField Factory can not create %s force field\n",
115	simParams->getForceField().c_str());
116	painCave.isFatal = 1;
117	simError();
118	}
119
120	if (simParams->haveForceFieldFileName()) {
121	ff->setForceFieldFileName(simParams->getForceFieldFileName());
122	}
123
124	std::string forcefieldFileName;
125	forcefieldFileName = ff->getForceFieldFileName();
126
127	if (simParams->haveForceFieldVariant()) {
128	//If the force field has variant, the variant force field name will be
129	//Base.variant.frc. For exampel EAM.u6.frc
130
131	std::string variant = simParams->getForceFieldVariant();
132
133	std::string::size_type pos = forcefieldFileName.rfind(".frc");
134	variant = "." + variant;
135	if (pos != std::string::npos) {
136	forcefieldFileName.insert(pos, variant);
137	} else {
138	//If the default force field file name does not containt .frc suffix, just append the .variant
139	forcefieldFileName.append(variant);
140	}
141	}
142
143	ff->parse(forcefieldFileName);
144
145	//extract the molecule stamps
146	std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
147	compList(stamps, simParams, moleculeStampPairs);
148
149	//create SimInfo
150	SimInfo * info = new SimInfo(stamps, moleculeStampPairs, ff, simParams);
151
152	//gather parameters (SimCreator only retrieves part of the parameters)
153	gatherParameters(info, mdFileName);
154
155	//divide the molecules and determine the global index of molecules
156	#ifdef IS_MPI
157	divideMolecules(info);
158	#endif
159
160	//create the molecules
161	createMolecules(info);
162
163
164	//allocate memory for DataStorage(circular reference, need to break it)
165	info->setSnapshotManager(new SimSnapshotManager(info));
166
167	//set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
168	//global index will never change again). Local indices of atoms and rigidbodies are already set by
169	//MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
170	setGlobalIndex(info);
171
172	//Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
173	//atoms don't have the global index yet (their global index are all initialized to -1).
174	//Therefore we have to call addExcludePairs explicitly here. A way to work around is that
175	//we can determine the beginning global indices of atoms before they get created.
176	SimInfo::MoleculeIterator mi;
177	Molecule* mol;
178	for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
179	info->addExcludePairs(mol);
180	}
181
182	if (loadInitCoords)
183	loadCoordinates(info);
184
185	return info;
186	}
187
188	void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
189
190	//figure out the output file names
191	std::string prefix;
192
193	#ifdef IS_MPI
194
195	if (worldRank == 0) {
196	#endif // is_mpi
197	Globals * simParams = info->getSimParams();
198	if (simParams->haveFinalConfig()) {
199	prefix = getPrefix(simParams->getFinalConfig());
200	} else {
201	prefix = getPrefix(mdfile);
202	}
203
204	info->setFinalConfigFileName(prefix + ".eor");
205	info->setDumpFileName(prefix + ".dump");
206	info->setStatFileName(prefix + ".stat");
207	info->setRestFileName(prefix + ".zang");
208
209	#ifdef IS_MPI
210
211	}
212
213	#endif
214
215	}
216
217	#ifdef IS_MPI
218	void SimCreator::divideMolecules(SimInfo *info) {
219	double numerator;
220	double denominator;
221	double precast;
222	double x;
223	double y;
224	double a;
225	int old_atoms;
226	int add_atoms;
227	int new_atoms;
228	int nTarget;
229	int done;
230	int i;
231	int j;
232	int loops;
233	int which_proc;
234	int nProcessors;
235	std::vector<int> atomsPerProc;
236	int nGlobalMols = info->getNGlobalMolecules();
237	std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
238
239	MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
240
241	if (nProcessors > nGlobalMols) {
242	sprintf(painCave.errMsg,
243	"nProcessors (%d) > nMol (%d)\n"
244	"\tThe number of processors is larger than\n"
245	"\tthe number of molecules. This will not result in a \n"
246	"\tusable division of atoms for force decomposition.\n"
247	"\tEither try a smaller number of processors, or run the\n"
248	"\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
249
250	painCave.isFatal = 1;
251	simError();
252	}
253
254	int seedValue;
255	Globals * simParams = info->getSimParams();
256	SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
257	if (simParams->haveSeed()) {
258	seedValue = simParams->getSeed();
259	myRandom = new SeqRandNumGen(seedValue);
260	}else {
261	myRandom = new SeqRandNumGen();
262	}
263
264
265	a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
266
267	//initialize atomsPerProc
268	atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
269
270	if (worldRank == 0) {
271	numerator = info->getNGlobalAtoms();
272	denominator = nProcessors;
273	precast = numerator / denominator;
274	nTarget = (int)(precast + 0.5);
275
276	for(i = 0; i < nGlobalMols; i++) {
277	done = 0;
278	loops = 0;
279
280	while (!done) {
281	loops++;
282
283	// Pick a processor at random
284
285	which_proc = (int) (myRandom->rand() * nProcessors);
286
287	//get the molecule stamp first
288	int stampId = info->getMoleculeStampId(i);
289	MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
290
291	// How many atoms does this processor have so far?
292	old_atoms = atomsPerProc[which_proc];
293	add_atoms = moleculeStamp->getNAtoms();
294	new_atoms = old_atoms + add_atoms;
295
296	// If we've been through this loop too many times, we need
297	// to just give up and assign the molecule to this processor
298	// and be done with it.
299
300	if (loops > 100) {
301	sprintf(painCave.errMsg,
302	"I've tried 100 times to assign molecule %d to a "
303	" processor, but can't find a good spot.\n"
304	"I'm assigning it at random to processor %d.\n",
305	i, which_proc);
306
307	painCave.isFatal = 0;
308	simError();
309
310	molToProcMap[i] = which_proc;
311	atomsPerProc[which_proc] += add_atoms;
312
313	done = 1;
314	continue;
315	}
316
317	// If we can add this molecule to this processor without sending
318	// it above nTarget, then go ahead and do it:
319
320	if (new_atoms <= nTarget) {
321	molToProcMap[i] = which_proc;
322	atomsPerProc[which_proc] += add_atoms;
323
324	done = 1;
325	continue;
326	}
327
328	// The only situation left is when new_atoms > nTarget. We
329	// want to accept this with some probability that dies off the
330	// farther we are from nTarget
331
332	// roughly: x = new_atoms - nTarget
333	// Pacc(x) = exp(- a * x)
334	// where a = penalty / (average atoms per molecule)
335
336	x = (double)(new_atoms - nTarget);
337	y = myRandom->rand();
338
339	if (y < exp(- a * x)) {
340	molToProcMap[i] = which_proc;
341	atomsPerProc[which_proc] += add_atoms;
342
343	done = 1;
344	continue;
345	} else {
346	continue;
347	}
348	}
349	}
350
351	delete myRandom;
352
353	// Spray out this nonsense to all other processors:
354
355	MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
356	} else {
357
358	// Listen to your marching orders from processor 0:
359
360	MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
361	}
362
363	info->setMolToProcMap(molToProcMap);
364	sprintf(checkPointMsg,
365	"Successfully divided the molecules among the processors.\n");
366	MPIcheckPoint();
367	}
368
369	#endif
370
371	void SimCreator::createMolecules(SimInfo *info) {
372	MoleculeCreator molCreator;
373	int stampId;
374
375	for(int i = 0; i < info->getNGlobalMolecules(); i++) {
376
377	#ifdef IS_MPI
378
379	if (info->getMolToProc(i) == worldRank) {
380	#endif
381
382	stampId = info->getMoleculeStampId(i);
383	Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
384	stampId, i, info->getLocalIndexManager());
385
386	info->addMolecule(mol);
387
388	#ifdef IS_MPI
389
390	}
391
392	#endif
393
394	} //end for(int i=0)
395	}
396
397	void SimCreator::compList(MakeStamps stamps, Globals simParams,
398	std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
399	int i;
400	char * id;
401	MoleculeStamp * currentStamp;
402	Component** the_components = simParams->getComponents();
403	int n_components = simParams->getNComponents();
404
405	if (!simParams->haveNMol()) {
406	// we don't have the total number of molecules, so we assume it is
407	// given in each component
408
409	for(i = 0; i < n_components; i++) {
410	if (!the_components[i]->haveNMol()) {
411	// we have a problem
412	sprintf(painCave.errMsg,
413	"SimCreator Error. No global NMol or component NMol given.\n"
414	"\tCannot calculate the number of atoms.\n");
415
416	painCave.isFatal = 1;
417	simError();
418	}
419
420	id = the_components[i]->getType();
421
422	currentStamp = stamps->getMolStamp(id);
423	if (currentStamp == NULL) {
424	sprintf(painCave.errMsg,
425	"SimCreator error: Component \"%s\" was not found in the "
426	"list of declared molecules\n", id);
427
428	painCave.isFatal = 1;
429	simError();
430	}
431
432	moleculeStampPairs.push_back(
433	std::make_pair(currentStamp, the_components[i]->getNMol()));
434	} //end for (i = 0; i < n_components; i++)
435	} else {
436	sprintf(painCave.errMsg, "SimSetup error.\n"
437	"\tSorry, the ability to specify total"
438	" nMols and then give molfractions in the components\n"
439	"\tis not currently supported."
440	" Please give nMol in the components.\n");
441
442	painCave.isFatal = 1;
443	simError();
444	}
445
446	#ifdef IS_MPI
447
448	strcpy(checkPointMsg, "Component stamps successfully extracted\n");
449	MPIcheckPoint();
450
451	#endif // is_mpi
452
453	}
454
455	void SimCreator::setGlobalIndex(SimInfo *info) {
456	SimInfo::MoleculeIterator mi;
457	Molecule::AtomIterator ai;
458	Molecule::RigidBodyIterator ri;
459	Molecule::CutoffGroupIterator ci;
460	Molecule * mol;
461	Atom * atom;
462	RigidBody * rb;
463	CutoffGroup * cg;
464	int beginAtomIndex;
465	int beginRigidBodyIndex;
466	int beginCutoffGroupIndex;
467	int nGlobalAtoms = info->getNGlobalAtoms();
468
469	#ifndef IS_MPI
470
471	beginAtomIndex = 0;
472	beginRigidBodyIndex = 0;
473	beginCutoffGroupIndex = 0;
474
475	#else
476
477	int nproc;
478	int myNode;
479
480	myNode = worldRank;
481	MPI_Comm_size(MPI_COMM_WORLD, &nproc);
482
483	std::vector < int > tmpAtomsInProc(nproc, 0);
484	std::vector < int > tmpRigidBodiesInProc(nproc, 0);
485	std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
486	std::vector < int > NumAtomsInProc(nproc, 0);
487	std::vector < int > NumRigidBodiesInProc(nproc, 0);
488	std::vector < int > NumCutoffGroupsInProc(nproc, 0);
489
490	tmpAtomsInProc[myNode] = info->getNAtoms();
491	tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
492	tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
493
494	//do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
495	MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
496	MPI_SUM, MPI_COMM_WORLD);
497	MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
498	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
499	MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
500	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
501
502	beginAtomIndex = 0;
503	beginRigidBodyIndex = 0;
504	beginCutoffGroupIndex = 0;
505
506	for(int i = 0; i < myNode; i++) {
507	beginAtomIndex += NumAtomsInProc[i];
508	beginRigidBodyIndex += NumRigidBodiesInProc[i];
509	beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
510	}
511
512	#endif
513
514	//rigidbody's index begins right after atom's
515	beginRigidBodyIndex += info->getNGlobalAtoms();
516
517	for(mol = info->beginMolecule(mi); mol != NULL;
518	mol = info->nextMolecule(mi)) {
519
520	//local index(index in DataStorge) of atom is important
521	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
522	atom->setGlobalIndex(beginAtomIndex++);
523	}
524
525	for(rb = mol->beginRigidBody(ri); rb != NULL;
526	rb = mol->nextRigidBody(ri)) {
527	rb->setGlobalIndex(beginRigidBodyIndex++);
528	}
529
530	//local index of cutoff group is trivial, it only depends on the order of travesing
531	for(cg = mol->beginCutoffGroup(ci); cg != NULL;
532	cg = mol->nextCutoffGroup(ci)) {
533	cg->setGlobalIndex(beginCutoffGroupIndex++);
534	}
535	}
536
537	//fill globalGroupMembership
538	std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
539	for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
540	for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
541
542	for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
543	globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
544	}
545
546	}
547	}
548
549	#ifdef IS_MPI
550	// Since the globalGroupMembership has been zero filled and we've only
551	// poked values into the atoms we know, we can do an Allreduce
552	// to get the full globalGroupMembership array (We think).
553	// This would be prettier if we could use MPI_IN_PLACE like the MPI-2
554	// docs said we could.
555	std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
556	MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
557	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
558	info->setGlobalGroupMembership(tmpGroupMembership);
559	#else
560	info->setGlobalGroupMembership(globalGroupMembership);
561	#endif
562
563	//fill molMembership
564	std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
565
566	for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
567
568	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
569	globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
570	}
571	}
572
573	#ifdef IS_MPI
574	std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
575
576	MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
577	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
578
579	info->setGlobalMolMembership(tmpMolMembership);
580	#else
581	info->setGlobalMolMembership(globalMolMembership);
582	#endif
583
584	}
585
586	void SimCreator::loadCoordinates(SimInfo* info) {
587	Globals* simParams;
588	simParams = info->getSimParams();
589
590	if (!simParams->haveInitialConfig()) {
591	sprintf(painCave.errMsg,
592	"Cannot intialize a simulation without an initial configuration file.\n");
593	painCave.isFatal = 1;;
594	simError();
595	}
596
597	DumpReader reader(info, simParams->getInitialConfig());
598	int nframes = reader.getNFrames();
599
600	if (nframes > 0) {
601	reader.readFrame(nframes - 1);
602	} else {
603	//invalid initial coordinate file
604	sprintf(painCave.errMsg,
605	"Initial configuration file %s should at least contain one frame\n",
606	simParams->getInitialConfig().c_str());
607	painCave.isFatal = 1;
608	simError();
609	}
610
611	//copy the current snapshot to previous snapshot
612	info->getSnapshotManager()->advance();
613	}
614
615	} //end namespace oopse
616
617