src/brains/SimCreator.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

/**
 * @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"
#ifdef IS_MPI
#include "io/mpiBASS.h"
#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) {
    
    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();
    }

    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);

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

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

    //setup seed for random number generator
    int seedValue;
    Globals * simParams = info->getSimParams();

    if (simParams->haveSeed()) {
        seedValue = simParams->getSeed();

        if (seedValue / 1000000000 == 0) {
            sprintf(painCave.errMsg,
                    "Seed for sprng library should contain at least 9 digits\n"
                        "OOPSE will generate a seed for user\n");

            painCave.isFatal = 0;
            simError();

            //using seed generated by system instead of invalid seed set by user 

#ifndef IS_MPI

            seedValue = make_sprng_seed();

#else

            if (worldRank == 0) {
                seedValue = make_sprng_seed();
            }

            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);

#endif

        } //end if (seedValue /1000000000 == 0)
    } else {

#ifndef IS_MPI

        seedValue = make_sprng_seed();

#else

        if (worldRank == 0) {
            seedValue = make_sprng_seed();
        }

        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);

#endif

    } //end of simParams->haveSeed()

    info->setSeed(seedValue);


    //figure out the ouput file names
    std::string prefix;

#ifdef IS_MPI

    if (worldRank == 0) {
#endif // is_mpi

        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;
    randomSPRNG myRandom(info->getSeed());
    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();
    }

    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.getRandom() * 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.getRandom();

                if (y < exp(- a * x)) {
                    molToProcMap[i] = which_proc;
                    atomsPerProc[which_proc] += add_atoms;

                    done = 1;
                    continue;
                } else {
                    continue;
                }
            }
        }

        // 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 < nproc; i++) {
        beginAtomIndex += NumAtomsInProc[i];
        beginRigidBodyIndex += NumRigidBodiesInProc[i];
        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
    }

#endif

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

}

} //end namespace oopse


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