ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/brains/SimCreator.cpp
Revision: 2211
Committed: Thu Apr 21 14:12:19 2005 UTC (19 years, 2 months ago) by chrisfen
File size: 19987 byte(s)
Log Message:
Shapes is limping along with a array bounds overwrite (I think...). At least the parser loads the forcefield fine...

File Contents

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