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 <sprng.h>

#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;
    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->extractMolStamp(id))->getStamp();

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

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

} //end namespace oopse


Revision:	393
Committed:	Wed Mar 2 16:28:20 2005 UTC (20 years, 2 months ago) by tim
File size:	20788 byte(s)
Log Message:	info_ in Velocitizer is not initialized which causes a seg fault
#	User	Rev	Content
1	tim	381	/*
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 <sprng.h>
51
52			#include "brains/MoleculeCreator.hpp"
53			#include "brains/SimCreator.hpp"
54			#include "brains/SimSnapshotManager.hpp"
55			#include "io/DumpReader.hpp"
56			#include "io/parse_me.h"
57			#include "UseTheForce/ForceFieldFactory.hpp"
58			#include "utils/simError.h"
59			#include "utils/StringUtils.hpp"
60	tim	393	#include "math/SeqRandNumGen.hpp"
61	tim	381	#ifdef IS_MPI
62			#include "io/mpiBASS.h"
63	tim	392	#include "math/ParallelRandNumGen.hpp"
64	tim	381	#endif
65
66			namespace oopse {
67
68			void SimCreator::parseFile(const std::string mdFileName, MakeStamps* stamps, Globals* simParams){
69
70			#ifdef IS_MPI
71
72			if (worldRank == 0) {
73			#endif // is_mpi
74
75			simParams->initalize();
76			set_interface_stamps(stamps, simParams);
77
78			#ifdef IS_MPI
79
80			mpiEventInit();
81
82			#endif
83
84			yacc_BASS(mdFileName.c_str());
85
86			#ifdef IS_MPI
87
88			throwMPIEvent(NULL);
89			} else {
90			set_interface_stamps(stamps, simParams);
91			mpiEventInit();
92			MPIcheckPoint();
93			mpiEventLoop();
94			}
95
96			#endif
97
98			}
99
100			SimInfo* SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
101
102			MakeStamps * stamps = new MakeStamps();
103
104			Globals * simParams = new Globals();
105
106			//parse meta-data file
107			parseFile(mdFileName, stamps, simParams);
108
109			//create the force field
110			ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
111			simParams->getForceField());
112
113			if (ff == NULL) {
114			sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
115			simParams->getForceField());
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(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
183			//load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
184			//eta, chi for NPT integrator)
185			if (loadInitCoords)
186			loadCoordinates(info);
187
188			return info;
189			}
190
191			void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
192
193			//figure out the ouput file names
194			std::string prefix;
195
196			#ifdef IS_MPI
197
198			if (worldRank == 0) {
199			#endif // is_mpi
200	tim	384	Globals * simParams = info->getSimParams();
201	tim	381	if (simParams->haveFinalConfig()) {
202			prefix = getPrefix(simParams->getFinalConfig());
203			} else {
204			prefix = getPrefix(mdfile);
205			}
206
207			info->setFinalConfigFileName(prefix + ".eor");
208			info->setDumpFileName(prefix + ".dump");
209			info->setStatFileName(prefix + ".stat");
210
211			#ifdef IS_MPI
212
213			}
214
215			#endif
216
217			}
218
219			#ifdef IS_MPI
220			void SimCreator::divideMolecules(SimInfo *info) {
221			double numerator;
222			double denominator;
223			double precast;
224			double x;
225			double y;
226			double a;
227			int old_atoms;
228			int add_atoms;
229			int new_atoms;
230			int nTarget;
231			int done;
232			int i;
233			int j;
234			int loops;
235			int which_proc;
236			int nProcessors;
237			std::vector<int> atomsPerProc;
238			int nGlobalMols = info->getNGlobalMolecules();
239			std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
240
241			MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
242
243			if (nProcessors > nGlobalMols) {
244			sprintf(painCave.errMsg,
245			"nProcessors (%d) > nMol (%d)\n"
246			"\tThe number of processors is larger than\n"
247			"\tthe number of molecules. This will not result in a \n"
248			"\tusable division of atoms for force decomposition.\n"
249			"\tEither try a smaller number of processors, or run the\n"
250			"\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
251
252			painCave.isFatal = 1;
253			simError();
254			}
255
256	tim	384	int seedValue;
257			Globals * simParams = info->getSimParams();
258	tim	392	SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
259	tim	384	if (simParams->haveSeed()) {
260			seedValue = simParams->getSeed();
261	tim	393	myRandom = new SeqRandNumGen(seedValue);
262	tim	384	}else {
263	tim	393	myRandom = new SeqRandNumGen();
264	tim	384	}
265	tim	381
266
267			a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
268
269			//initialize atomsPerProc
270			atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
271
272			if (worldRank == 0) {
273			numerator = info->getNGlobalAtoms();
274			denominator = nProcessors;
275			precast = numerator / denominator;
276			nTarget = (int)(precast + 0.5);
277
278			for(i = 0; i < nGlobalMols; i++) {
279			done = 0;
280			loops = 0;
281
282			while (!done) {
283			loops++;
284
285			// Pick a processor at random
286
287	tim	384	which_proc = (int) (myRandom->rand() * nProcessors);
288	tim	381
289			//get the molecule stamp first
290			int stampId = info->getMoleculeStampId(i);
291			MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
292
293			// How many atoms does this processor have so far?
294			old_atoms = atomsPerProc[which_proc];
295			add_atoms = moleculeStamp->getNAtoms();
296			new_atoms = old_atoms + add_atoms;
297
298			// If we've been through this loop too many times, we need
299			// to just give up and assign the molecule to this processor
300			// and be done with it.
301
302			if (loops > 100) {
303			sprintf(painCave.errMsg,
304			"I've tried 100 times to assign molecule %d to a "
305			" processor, but can't find a good spot.\n"
306			"I'm assigning it at random to processor %d.\n",
307			i, which_proc);
308
309			painCave.isFatal = 0;
310			simError();
311
312			molToProcMap[i] = which_proc;
313			atomsPerProc[which_proc] += add_atoms;
314
315			done = 1;
316			continue;
317			}
318
319			// If we can add this molecule to this processor without sending
320			// it above nTarget, then go ahead and do it:
321
322			if (new_atoms <= nTarget) {
323			molToProcMap[i] = which_proc;
324			atomsPerProc[which_proc] += add_atoms;
325
326			done = 1;
327			continue;
328			}
329
330			// The only situation left is when new_atoms > nTarget. We
331			// want to accept this with some probability that dies off the
332			// farther we are from nTarget
333
334			// roughly: x = new_atoms - nTarget
335			// Pacc(x) = exp(- a * x)
336			// where a = penalty / (average atoms per molecule)
337
338			x = (double)(new_atoms - nTarget);
339	tim	384	y = myRandom->rand();
340	tim	381
341			if (y < exp(- a * x)) {
342			molToProcMap[i] = which_proc;
343			atomsPerProc[which_proc] += add_atoms;
344
345			done = 1;
346			continue;
347			} else {
348			continue;
349			}
350			}
351			}
352
353	tim	384	delete myRandom;
354
355	tim	381	// Spray out this nonsense to all other processors:
356
357			MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
358			} else {
359
360			// Listen to your marching orders from processor 0:
361
362			MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
363			}
364
365			info->setMolToProcMap(molToProcMap);
366			sprintf(checkPointMsg,
367			"Successfully divided the molecules among the processors.\n");
368			MPIcheckPoint();
369			}
370
371			#endif
372
373			void SimCreator::createMolecules(SimInfo *info) {
374			MoleculeCreator molCreator;
375			int stampId;
376
377			for(int i = 0; i < info->getNGlobalMolecules(); i++) {
378
379			#ifdef IS_MPI
380
381			if (info->getMolToProc(i) == worldRank) {
382			#endif
383
384			stampId = info->getMoleculeStampId(i);
385			Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
386			stampId, i, info->getLocalIndexManager());
387
388			info->addMolecule(mol);
389
390			#ifdef IS_MPI
391
392			}
393
394			#endif
395
396			} //end for(int i=0)
397			}
398
399			void SimCreator::compList(MakeStamps stamps, Globals simParams,
400			std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
401			int i;
402			char * id;
403			MoleculeStamp * currentStamp;
404			Component** the_components = simParams->getComponents();
405			int n_components = simParams->getNComponents();
406
407			if (!simParams->haveNMol()) {
408			// we don't have the total number of molecules, so we assume it is
409			// given in each component
410
411			for(i = 0; i < n_components; i++) {
412			if (!the_components[i]->haveNMol()) {
413			// we have a problem
414			sprintf(painCave.errMsg,
415			"SimCreator Error. No global NMol or component NMol given.\n"
416			"\tCannot calculate the number of atoms.\n");
417
418			painCave.isFatal = 1;
419			simError();
420			}
421
422			id = the_components[i]->getType();
423			currentStamp = (stamps->extractMolStamp(id))->getStamp();
424
425			if (currentStamp == NULL) {
426			sprintf(painCave.errMsg,
427			"SimCreator error: Component \"%s\" was not found in the "
428			"list of declared molecules\n", id);
429
430			painCave.isFatal = 1;
431			simError();
432			}
433
434			moleculeStampPairs.push_back(
435			std::make_pair(currentStamp, the_components[i]->getNMol()));
436			} //end for (i = 0; i < n_components; i++)
437			} else {
438			sprintf(painCave.errMsg, "SimSetup error.\n"
439			"\tSorry, the ability to specify total"
440			" nMols and then give molfractions in the components\n"
441			"\tis not currently supported."
442			" Please give nMol in the components.\n");
443
444			painCave.isFatal = 1;
445			simError();
446			}
447
448			#ifdef IS_MPI
449
450			strcpy(checkPointMsg, "Component stamps successfully extracted\n");
451			MPIcheckPoint();
452
453			#endif // is_mpi
454
455			}
456
457			void SimCreator::setGlobalIndex(SimInfo *info) {
458			SimInfo::MoleculeIterator mi;
459			Molecule::AtomIterator ai;
460			Molecule::RigidBodyIterator ri;
461			Molecule::CutoffGroupIterator ci;
462			Molecule * mol;
463			Atom * atom;
464			RigidBody * rb;
465			CutoffGroup * cg;
466			int beginAtomIndex;
467			int beginRigidBodyIndex;
468			int beginCutoffGroupIndex;
469			int nGlobalAtoms = info->getNGlobalAtoms();
470
471			#ifndef IS_MPI
472
473			beginAtomIndex = 0;
474			beginRigidBodyIndex = 0;
475			beginCutoffGroupIndex = 0;
476
477			#else
478
479			int nproc;
480			int myNode;
481
482			myNode = worldRank;
483			MPI_Comm_size(MPI_COMM_WORLD, &nproc);
484
485			std::vector < int > tmpAtomsInProc(nproc, 0);
486			std::vector < int > tmpRigidBodiesInProc(nproc, 0);
487			std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
488			std::vector < int > NumAtomsInProc(nproc, 0);
489			std::vector < int > NumRigidBodiesInProc(nproc, 0);
490			std::vector < int > NumCutoffGroupsInProc(nproc, 0);
491
492			tmpAtomsInProc[myNode] = info->getNAtoms();
493			tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
494			tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
495
496			//do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
497			MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
498			MPI_SUM, MPI_COMM_WORLD);
499			MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
500			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
501			MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
502			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
503
504			beginAtomIndex = 0;
505			beginRigidBodyIndex = 0;
506			beginCutoffGroupIndex = 0;
507
508			for(int i = 0; i < myNode; i++) {
509			beginAtomIndex += NumAtomsInProc[i];
510			beginRigidBodyIndex += NumRigidBodiesInProc[i];
511			beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
512			}
513
514			#endif
515
516			//rigidbody's index begins right after atom's
517			beginRigidBodyIndex += info->getNGlobalAtoms();
518
519			for(mol = info->beginMolecule(mi); mol != NULL;
520			mol = info->nextMolecule(mi)) {
521
522			//local index(index in DataStorge) of atom is important
523			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
524			atom->setGlobalIndex(beginAtomIndex++);
525			}
526
527			for(rb = mol->beginRigidBody(ri); rb != NULL;
528			rb = mol->nextRigidBody(ri)) {
529			rb->setGlobalIndex(beginRigidBodyIndex++);
530			}
531
532			//local index of cutoff group is trivial, it only depends on the order of travesing
533			for(cg = mol->beginCutoffGroup(ci); cg != NULL;
534			cg = mol->nextCutoffGroup(ci)) {
535			cg->setGlobalIndex(beginCutoffGroupIndex++);
536			}
537			}
538
539			//fill globalGroupMembership
540			std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
541			for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
542			for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
543
544			for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
545			globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
546			}
547
548			}
549			}
550
551			#ifdef IS_MPI
552			// Since the globalGroupMembership has been zero filled and we've only
553			// poked values into the atoms we know, we can do an Allreduce
554			// to get the full globalGroupMembership array (We think).
555			// This would be prettier if we could use MPI_IN_PLACE like the MPI-2
556			// docs said we could.
557			std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
558			MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
559			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
560			info->setGlobalGroupMembership(tmpGroupMembership);
561			#else
562			info->setGlobalGroupMembership(globalGroupMembership);
563			#endif
564
565			//fill molMembership
566			std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
567
568			for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
569
570			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
571			globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
572			}
573			}
574
575			#ifdef IS_MPI
576			std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
577
578			MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
579			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
580
581			info->setGlobalMolMembership(tmpMolMembership);
582			#else
583			info->setGlobalMolMembership(globalMolMembership);
584			#endif
585
586			}
587
588			void SimCreator::loadCoordinates(SimInfo* info) {
589			Globals* simParams;
590			simParams = info->getSimParams();
591
592			if (!simParams->haveInitialConfig()) {
593			sprintf(painCave.errMsg,
594			"Cannot intialize a simulation without an initial configuration file.\n");
595			painCave.isFatal = 1;;
596			simError();
597			}
598
599			DumpReader reader(info, simParams->getInitialConfig());
600			int nframes = reader.getNFrames();
601
602			if (nframes > 0) {
603			reader.readFrame(nframes - 1);
604			} else {
605			//invalid initial coordinate file
606			sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
607			simParams->getInitialConfig());
608			painCave.isFatal = 1;
609			simError();
610			}
611
612			//copy the current snapshot to previous snapshot
613			info->getSnapshotManager()->advance();
614			}
615
616			} //end namespace oopse
617
618