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

      simParams->initalize();
      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());
      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(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);
    }
    

    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
    //eta, chi for NPT integrator)
    if (loadInitCoords)
      loadCoordinates(info);    
    
    return info;
  }

  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {

    //figure out the ouput 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");

#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;
    LinkedMolStamp* extractedStamp = NULL;
    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();

        extractedStamp = stamps->extractMolStamp(id);
        if (extractedStamp == NULL) {
          sprintf(painCave.errMsg,
                  "SimCreator error: Component \"%s\" was not found in the "
                  "list of declared molecules\n", id);

          painCave.isFatal = 1;
          simError();
        }

        currentStamp = extractedStamp->getStamp();


        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());
      painCave.isFatal = 1;
      simError();
    }

    //copy the current snapshot to previous snapshot
    info->getSnapshotManager()->advance();
  }

} //end namespace oopse


Revision:	2087
Committed:	Tue Mar 8 21:06:49 2005 UTC (19 years, 4 months ago) by gezelter
File size:	19491 byte(s)
Log Message:	Using OOPSE_itoa and checking a NULL left over from the LinkedList days
#	User	Rev	Content
1	gezelter	2087	/*
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, Globals* simParams){
67
68			#ifdef IS_MPI
69
70			if (worldRank == 0) {
71			#endif // is_mpi
72
73			simParams->initalize();
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, bool loadInitCoords) {
99
100			MakeStamps * stamps = new MakeStamps();
101
102			Globals * simParams = new Globals();
103
104			//parse meta-data file
105			parseFile(mdFileName, stamps, simParams);
106
107			//create the force field
108			ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
109			simParams->getForceField());
110
111			if (ff == NULL) {
112			sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
113			simParams->getForceField());
114			painCave.isFatal = 1;
115			simError();
116			}
117
118			if (simParams->haveForceFieldFileName()) {
119			ff->setForceFieldFileName(simParams->getForceFieldFileName());
120			}
121
122			std::string forcefieldFileName;
123			forcefieldFileName = ff->getForceFieldFileName();
124
125			if (simParams->haveForceFieldVariant()) {
126			//If the force field has variant, the variant force field name will be
127			//Base.variant.frc. For exampel EAM.u6.frc
128
129			std::string variant = simParams->getForceFieldVariant();
130
131			std::string::size_type pos = forcefieldFileName.rfind(".frc");
132			variant = "." + variant;
133			if (pos != std::string::npos) {
134			forcefieldFileName.insert(pos, variant);
135			} else {
136			//If the default force field file name does not containt .frc suffix, just append the .variant
137			forcefieldFileName.append(variant);
138			}
139			}
140
141			ff->parse(forcefieldFileName);
142
143			//extract the molecule stamps
144			std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
145			compList(stamps, simParams, moleculeStampPairs);
146
147			//create SimInfo
148			SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
149
150			//gather parameters (SimCreator only retrieves part of the parameters)
151			gatherParameters(info, mdFileName);
152
153			//divide the molecules and determine the global index of molecules
154			#ifdef IS_MPI
155			divideMolecules(info);
156			#endif
157
158			//create the molecules
159			createMolecules(info);
160
161
162			//allocate memory for DataStorage(circular reference, need to break it)
163			info->setSnapshotManager(new SimSnapshotManager(info));
164
165			//set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
166			//global index will never change again). Local indices of atoms and rigidbodies are already set by
167			//MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
168			setGlobalIndex(info);
169
170			//Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
171			//atoms don't have the global index yet (their global index are all initialized to -1).
172			//Therefore we have to call addExcludePairs explicitly here. A way to work around is that
173			//we can determine the beginning global indices of atoms before they get created.
174			SimInfo::MoleculeIterator mi;
175			Molecule* mol;
176			for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
177			info->addExcludePairs(mol);
178			}
179
180
181			//load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
182			//eta, chi for NPT integrator)
183			if (loadInitCoords)
184			loadCoordinates(info);
185
186			return info;
187			}
188
189			void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
190
191			//figure out the ouput file names
192			std::string prefix;
193
194			#ifdef IS_MPI
195
196			if (worldRank == 0) {
197			#endif // is_mpi
198			Globals * simParams = info->getSimParams();
199			if (simParams->haveFinalConfig()) {
200			prefix = getPrefix(simParams->getFinalConfig());
201			} else {
202			prefix = getPrefix(mdfile);
203			}
204
205			info->setFinalConfigFileName(prefix + ".eor");
206			info->setDumpFileName(prefix + ".dump");
207			info->setStatFileName(prefix + ".stat");
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			LinkedMolStamp* extractedStamp = NULL;
402			MoleculeStamp * currentStamp;
403			Component** the_components = simParams->getComponents();
404			int n_components = simParams->getNComponents();
405
406			if (!simParams->haveNMol()) {
407			// we don't have the total number of molecules, so we assume it is
408			// given in each component
409
410			for(i = 0; i < n_components; i++) {
411			if (!the_components[i]->haveNMol()) {
412			// we have a problem
413			sprintf(painCave.errMsg,
414			"SimCreator Error. No global NMol or component NMol given.\n"
415			"\tCannot calculate the number of atoms.\n");
416
417			painCave.isFatal = 1;
418			simError();
419			}
420
421			id = the_components[i]->getType();
422
423			extractedStamp = stamps->extractMolStamp(id);
424			if (extractedStamp == NULL) {
425			sprintf(painCave.errMsg,
426			"SimCreator error: Component \"%s\" was not found in the "
427			"list of declared molecules\n", id);
428
429			painCave.isFatal = 1;
430			simError();
431			}
432
433			currentStamp = extractedStamp->getStamp();
434
435
436			moleculeStampPairs.push_back(
437			std::make_pair(currentStamp, the_components[i]->getNMol()));
438			} //end for (i = 0; i < n_components; i++)
439			} else {
440			sprintf(painCave.errMsg, "SimSetup error.\n"
441			"\tSorry, the ability to specify total"
442			" nMols and then give molfractions in the components\n"
443			"\tis not currently supported."
444			" Please give nMol in the components.\n");
445
446			painCave.isFatal = 1;
447			simError();
448			}
449
450			#ifdef IS_MPI
451
452			strcpy(checkPointMsg, "Component stamps successfully extracted\n");
453			MPIcheckPoint();
454
455			#endif // is_mpi
456
457			}
458
459			void SimCreator::setGlobalIndex(SimInfo *info) {
460			SimInfo::MoleculeIterator mi;
461			Molecule::AtomIterator ai;
462			Molecule::RigidBodyIterator ri;
463			Molecule::CutoffGroupIterator ci;
464			Molecule * mol;
465			Atom * atom;
466			RigidBody * rb;
467			CutoffGroup * cg;
468			int beginAtomIndex;
469			int beginRigidBodyIndex;
470			int beginCutoffGroupIndex;
471			int nGlobalAtoms = info->getNGlobalAtoms();
472
473			#ifndef IS_MPI
474
475			beginAtomIndex = 0;
476			beginRigidBodyIndex = 0;
477			beginCutoffGroupIndex = 0;
478
479			#else
480
481			int nproc;
482			int myNode;
483
484			myNode = worldRank;
485			MPI_Comm_size(MPI_COMM_WORLD, &nproc);
486
487			std::vector < int > tmpAtomsInProc(nproc, 0);
488			std::vector < int > tmpRigidBodiesInProc(nproc, 0);
489			std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
490			std::vector < int > NumAtomsInProc(nproc, 0);
491			std::vector < int > NumRigidBodiesInProc(nproc, 0);
492			std::vector < int > NumCutoffGroupsInProc(nproc, 0);
493
494			tmpAtomsInProc[myNode] = info->getNAtoms();
495			tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
496			tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
497
498			//do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
499			MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
500			MPI_SUM, MPI_COMM_WORLD);
501			MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
502			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
503			MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
504			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
505
506			beginAtomIndex = 0;
507			beginRigidBodyIndex = 0;
508			beginCutoffGroupIndex = 0;
509
510			for(int i = 0; i < myNode; i++) {
511			beginAtomIndex += NumAtomsInProc[i];
512			beginRigidBodyIndex += NumRigidBodiesInProc[i];
513			beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
514			}
515
516			#endif
517
518			//rigidbody's index begins right after atom's
519			beginRigidBodyIndex += info->getNGlobalAtoms();
520
521			for(mol = info->beginMolecule(mi); mol != NULL;
522			mol = info->nextMolecule(mi)) {
523
524			//local index(index in DataStorge) of atom is important
525			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
526			atom->setGlobalIndex(beginAtomIndex++);
527			}
528
529			for(rb = mol->beginRigidBody(ri); rb != NULL;
530			rb = mol->nextRigidBody(ri)) {
531			rb->setGlobalIndex(beginRigidBodyIndex++);
532			}
533
534			//local index of cutoff group is trivial, it only depends on the order of travesing
535			for(cg = mol->beginCutoffGroup(ci); cg != NULL;
536			cg = mol->nextCutoffGroup(ci)) {
537			cg->setGlobalIndex(beginCutoffGroupIndex++);
538			}
539			}
540
541			//fill globalGroupMembership
542			std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
543			for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
544			for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
545
546			for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
547			globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
548			}
549
550			}
551			}
552
553			#ifdef IS_MPI
554			// Since the globalGroupMembership has been zero filled and we've only
555			// poked values into the atoms we know, we can do an Allreduce
556			// to get the full globalGroupMembership array (We think).
557			// This would be prettier if we could use MPI_IN_PLACE like the MPI-2
558			// docs said we could.
559			std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
560			MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
561			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
562			info->setGlobalGroupMembership(tmpGroupMembership);
563			#else
564			info->setGlobalGroupMembership(globalGroupMembership);
565			#endif
566
567			//fill molMembership
568			std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
569
570			for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
571
572			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
573			globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
574			}
575			}
576
577			#ifdef IS_MPI
578			std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
579
580			MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
581			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
582
583			info->setGlobalMolMembership(tmpMolMembership);
584			#else
585			info->setGlobalMolMembership(globalMolMembership);
586			#endif
587
588			}
589
590			void SimCreator::loadCoordinates(SimInfo* info) {
591			Globals* simParams;
592			simParams = info->getSimParams();
593
594			if (!simParams->haveInitialConfig()) {
595			sprintf(painCave.errMsg,
596			"Cannot intialize a simulation without an initial configuration file.\n");
597			painCave.isFatal = 1;;
598			simError();
599			}
600
601			DumpReader reader(info, simParams->getInitialConfig());
602			int nframes = reader.getNFrames();
603
604			if (nframes > 0) {
605			reader.readFrame(nframes - 1);
606			} else {
607			//invalid initial coordinate file
608			sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
609			simParams->getInitialConfig());
610			painCave.isFatal = 1;
611			simError();
612			}
613
614			//copy the current snapshot to previous snapshot
615			info->getSnapshotManager()->advance();
616			}
617
618			} //end namespace oopse
619
620