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 "brains/SimCreator.hpp"

namespace oopse {

void SimCreator::parseFile(const std::string& mdFileName, MakeStamps *stamps, Globals *globals) {

#ifdef IS_MPI

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

        globals->initalize();
        set_interface_stamps(stamps, globals);

#ifdef IS_MPI

        mpiEventInit();

#endif

        yacc_BASS(mdFileName.c_str());

#ifdef IS_MPI

        throwMPIEvent(NULL);
    } else {
        set_interface_stamps(stamps, globals);
        mpiEventInit();
        MPIcheckPoint();
        mpiEventLoop();
    }

#endif

}

SimInfo*  SimCreator::createSim(const std::string & mdFileName) {

    mdFileName_ = mdFileName;
    
    MakeStamps * stamps = new MakeStamps();

    Globals * globals = new Globals();

    //parse meta-data file
    parseFile(mdFileName_, stamps, globals);

    //create the force field
    ForceFiled * ff = ForceFieldFactory::getInstance()->createObject(
                          globals->getForceField());
    
    if (ff == NULL) {
        sprintf(painCave.errMsg, "ForceFiled Factory can not create %s force field\n",
                globals->getForceField());
        painCave.isFatal = 1;
        simError();
    }
    ff->parse();

    //extract the molecule stamps
    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
    compList(stamps, globals, moleculeStampPairs);

    //create SimInfo
    SimInfo * info = new SimInfo(moleculeStampPairs, ff, globals);

    //gather parameters (SimCreator only retrieves part of the parameters)
    gatherParameters(info);

    //divide the molecules and determine the global index of molecules
    divideMolecules(info);

    //create the molecules
    createMolecules(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. 
    setGlobalIndices(info);

    //allocate memory for DataStorage(circular reference, need to break it)
    info->setSnapshotManager(new SimSnapshotManager(info));

    //initialize fortran -- setup the cutoff 
    initFortran(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) {

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

    if (globals->haveSeed()) {
        seedValue = globals->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 globals->haveSeed()

    info->setSeed(seedValue);


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

#ifdef IS_MPI

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

        if (globals->haveFinalConfig()) {
            prefix = StringUtils::getPrefix(globals->getFinalConfig());
        } else {
            prefix = StringUtils::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;
    int nGlobalMols;
    std::vector < int > atomsPerProc;
    randomSPRNG myRandom(info->getSeed());
    int * molToProcMap = info->getMolToProcMapPointer();

    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
    nGlobalMols = info->getNGlobalMolecules();

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

    for(i = 0; i < nGlobalMols; i++) {
        // default to an error condition:
        molToProcMap[i] = -1;
    }

    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, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    } else {

        // Listen to your marching orders from processor 0:

        MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    }

    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->addMolecule(mol);

#ifdef IS_MPI

        }

#endif

    } //end for(int i=0)   
}

void SimSetup::compList(MakeStamps *stamps, Globals* globals,
                        std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
    int i;
    char * id;
    MoleculeStamp * currentStamp;
    Component * the_components = globals->getComponents();
    int n_components = globals->getNComponents();

    if (!globals->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,
                        "SimSetup 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);

            if (currentStamp == NULL) {
                sprintf(painCave.errMsg,
                        "SimSetup error: Component \"%s\" was not found in the "
                            "list of declared molecules\n", id);

                painCave.isFatal = 1;
                simError();
            }

            moleculeStampPairs.push_back(
                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) {
    typename SimInfo::MoleculeIterator mi;
    typename Molecule::AtomIterator ai;
    typename Molecule::RigidBodyIterator ri;
    typename Molecule::CutoffGroupIterator ci;
    Molecule * mol;
    Atom * atom;
    RigidBody * rb;
    CutoffGroup * cg;
    int beginAtomIndex;
    int beginRigidBodyIndex;
    int beginCutoffGroupIndex;

#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 = beginMolecule(mi); mol != NULL; mol = 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.

    MPI_Allreduce(&globalGroupMembership[0],
                  info->getGlobalGroupMembershipPointer(),
                  info->getNGlobalAtoms(),
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#else
    std::copy(globalGroupMembership.begin(), globalGroupMembership.end(), 
        info->getGlobalGroupMembershipPointer());
#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
    MPI_Allreduce(&globalMolMembership[0],
                  info->getGlobalMolMembershipPointer(),
                  info->getNGlobalAtoms(),
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#else
    std::copy(globalMolMembership.begin(), globalMolMembership.end(), 
        info->getGlobalMolMembershipPointer());
#endif

}

void SimCreator::loadCoordinates(SimInfo* info) {
    Globals* globals;
    globals = info->getGlobals();
    
    if (!globals->haveInitialConfig()) {
        sprintf(painCave.errMsg,
                "Cannot intialize a simulation without an initial configuration file.\n");
        painCave.isFatal = 1;;
        simError();
    }
        
    DumpReader reader(globals->getInitialConfig());
    int nframes = reader->getNframes();

    if (nframes > 0) {
        reader.readFrame(info, nframes - 1);
    } else {
        //invalid initial coordinate file
        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
                globals->getInitialConfig());
        painCave.isFatal = 1;
        simError();
    }

}

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