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 SimSetup::parseFile(char* mdfile,  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(mdfile);

#ifdef IS_MPI
    throwMPIEvent(NULL);
  }
  else{
  set_interface_stamps(stamps, globals);
  mpiEventInit();
  MPIcheckPoint();
  mpiEventLoop();
  }
#endif

}


SimInfo* SimCreator::createSim(const std::string& mdfile) {
    MakeStamps* stamps = new MakeStamps();
    
    Globals* globals = new Globals();

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

    //create the force field
    ForceFiled* ff = ForceFieldFactory::getInstance()->createObject(globals->getForceField());
    
    //create SimInfo
    SimInfo* info = new SimInfo();
    info->setGlobals(globals);
    info->setForceField(ff);

    //extract the molecule stamps and add them into SimInfo
    compList(stamps, info);
    
    //gather parameters (SimCreator only retrieves the parameters which will be used to create
    // the simulation)
    gatherParameters(info);
    
    //divide the molecules and determine the global index of molecules
    
    //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);
    
    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
    //eta, chi for NPT integrator)
    DumpReader reader(info->getInitFilename());
    int nframes = reader->getNframes();

    if (nframes > 0) {
        reader.readFrame(info, nframes - 1);
    } else {
        //invalid initial coordinate file
        
    }

    
    //initialize fortran    
    
    return info;
}    

void SimCreator::gatherParameters(SimInfo* info) {
    //model->addProperty(new StringGenericData("Ensemble", globals->getForceFiled()));
    //model->addProperty(new DoubleGenericData("dt"), globals->getDt());

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

    //
    std::string prefix;
#ifdef IS_MPI
    if (worldRank == 0){
#endif // is_mpi
      
        if(globals->haveFinalConfig()) {
            prefix = getPrefix(globals->getFinalConfig());   
        } else {
            prefix = getPrefix(mdfile_);
        }
        
        info->setFinalConfFileName(prefix + ".eor");
        info->setDumpFileName(prefix + ".dump");
        info->setStatFileName(prefix + ".stat");

#ifdef IS_MPI
    }
#endif    


}

#ifdef IS_MPI
void SimCreator::divideMolecules(SimInfo* info){
    int nComponents;
    MoleculeStamp ** compStamps;
    randomSPRNG * myRandom;
    int * componentsNmol;
    int * AtomsPerProc;

    double numerator;
    double denominator;
    double precast;
    double x;
    double y;
    double a;
    int old_atoms;
    int add_atoms;
    int new_atoms;

    int nTarget;
    int molIndex;
    int atomIndex;
    int done;
    int i;
    int j;
    int loops;
    int which_proc;
    int nmol_global,
    nmol_local;
    int natoms_global,

    int baseSeed = info->getSeed();
    CutoffGroupStamp * cg;

    nComponents = info->nComponents;
    compStamps = info->compStamps;
    componentsNmol = info->componentsNmol;
    AtomsPerProc = new int[parallelData->nProcessors];

    parallelData->nMolGlobal = info->n_mol;

    if (parallelData->nProcessors > parallelData->nMolGlobal) {
        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",
                parallelData->nProcessors,
                parallelData->nMolGlobal);

        painCave.isFatal = 1;
        simError();
    }

    myRandom = new randomSPRNG(baseSeed);

    a = 3.0 * parallelData->nMolGlobal / parallelData->nAtomsGlobal;

    // Initialize things that we'll send out later:
    for( i = 0; i < parallelData->nProcessors; i++ ) {
        AtomsPerProc[i] = 0;
    }

    for( i = 0; i < parallelData->nMolGlobal; i++ ) {
        // default to an error condition:
        MolToProcMap[i] = -1;
    }

    if (parallelData->myNode == 0) {
        numerator = info->n_atoms;
        denominator = parallelData->nProcessors;
        precast = numerator / denominator;
        nTarget = (int)(precast + 0.5);

        // Build the array of molecule component types first
        molIndex = 0;

        for( i = 0; i < nComponents; i++ ) {
            for( j = 0; j < componentsNmol[i]; j++ ) {
                molIndex++;
            }
        }

        for( i = 0; i < molIndex; i++ ) {
            done = 0;
            loops = 0;

            while (!done) {
                loops++;

                // Pick a processor at random

                which_proc

                = (int) (myRandom->getRandom() * parallelData->nProcessors);

                // How many atoms does this processor have?

                old_atoms = AtomsPerProc[which_proc];
                add_atoms = compStamps[MolComponentType[i]]->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, parallelData->nMolGlobal, MPI_INT, 0,
                  MPI_COMM_WORLD);
    } else {

    // Listen to your marching orders from processor 0:

        MPI_Bcast(MolToProcMap, parallelData->nMolGlobal, MPI_INT, 0,
                  MPI_COMM_WORLD); }

    // Let's all check for sanity:

    nmol_local = 0;

    for( i = 0; i < parallelData->nMolGlobal; i++ ) {
        if (MolToProcMap[i] == parallelData->myNode) {
            nmol_local++;
        }
    }

    MPI_Allreduce(&nmol_local, &nmol_global, 1, MPI_INT, MPI_SUM,
                  MPI_COMM_WORLD);

    if (nmol_global != info->n_mol) {
        sprintf(painCave.errMsg,
                "The sum of all nmol_local, %d, did not equal the "
                    "total number of molecules, %d.\n",
                nmol_global,
                info->n_mol);

        painCave.isFatal = 1;
        simError();
    }

    sprintf(checkPointMsg,
            "Successfully divided the molecules among the processors.\n");
    MPIcheckPoint();
 
}
#endif

Molecule* SimCreator::createMolecules(SimInfo* info) {
    MoleculeCreator molCreator;
    int stampId;

    
    for (int i = 0; i < info->getNGlobalMolecules(); i++){
        if (info->getMolToProc(i) == worldRank){

            stampId = info->getMoleculeStampId(i);
            Molecule* mol = molCreator.createMolecule(info->getForceField(), 
                                                                                info->getMoleculeStamp(stampId), stampId, i);
                        
            info->addMolecule(mol);
        }
    }        
       
}

void SimSetup::compList(MakeStamps* stamps,SimInfo* info) {
    int i;
    char* id;
    MoleculeStamp* currentStamp;
    Component* the_components = info->getGlobals()->getComponents();
    int n_components = info->getGlobals()->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();
            }      

            info->addMoleculeStamp(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::setGlobalIndices(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;

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

    tmpAtomInProc[myNode] = info->getNAtoms();
    tmpRigidBodiesInProc[myNode] = info->getNRigidBodiess();
    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroupss();
    
    //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++);
        }        
        
    }
}

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