ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/brains/SimCreator.cpp
Revision: 2072
Committed: Tue Mar 1 23:02:33 2005 UTC (19 years, 4 months ago) by tim
File size: 20729 byte(s)
Log Message:
replace SPRNG by MTRand

File Contents

# User Rev Content
1 tim 2065 /*
2     * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9     * 1. Acknowledgement of the program authors must be made in any
10     * publication of scientific results based in part on use of the
11     * program. An acceptable form of acknowledgement is citation of
12     * the article in which the program was described (Matthew
13     * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     * Parallel Simulation Engine for Molecular Dynamics,"
16     * J. Comput. Chem. 26, pp. 252-271 (2005))
17     *
18     * 2. Redistributions of source code must retain the above copyright
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41    
42     /**
43     * @file SimCreator.cpp
44     * @author tlin
45     * @date 11/03/2004
46     * @time 13:51am
47     * @version 1.0
48     */
49    
50     #include <sprng.h>
51    
52     #include "brains/MoleculeCreator.hpp"
53     #include "brains/SimCreator.hpp"
54     #include "brains/SimSnapshotManager.hpp"
55     #include "io/DumpReader.hpp"
56     #include "io/parse_me.h"
57     #include "UseTheForce/ForceFieldFactory.hpp"
58     #include "utils/simError.h"
59     #include "utils/StringUtils.hpp"
60 tim 2072 #include "math/MersenneTwister.hpp"
61 tim 2065 #ifdef IS_MPI
62     #include "io/mpiBASS.h"
63     #endif
64    
65     namespace oopse {
66    
67     void SimCreator::parseFile(const std::string mdFileName, MakeStamps* stamps, Globals* simParams){
68    
69     #ifdef IS_MPI
70    
71     if (worldRank == 0) {
72     #endif // is_mpi
73    
74     simParams->initalize();
75     set_interface_stamps(stamps, simParams);
76    
77     #ifdef IS_MPI
78    
79     mpiEventInit();
80    
81     #endif
82    
83     yacc_BASS(mdFileName.c_str());
84    
85     #ifdef IS_MPI
86    
87     throwMPIEvent(NULL);
88     } else {
89     set_interface_stamps(stamps, simParams);
90     mpiEventInit();
91     MPIcheckPoint();
92     mpiEventLoop();
93     }
94    
95     #endif
96    
97     }
98    
99     SimInfo* SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
100    
101     MakeStamps * stamps = new MakeStamps();
102    
103     Globals * simParams = new Globals();
104    
105     //parse meta-data file
106     parseFile(mdFileName, stamps, simParams);
107    
108     //create the force field
109     ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
110     simParams->getForceField());
111    
112     if (ff == NULL) {
113     sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
114     simParams->getForceField());
115     painCave.isFatal = 1;
116     simError();
117     }
118    
119     if (simParams->haveForceFieldFileName()) {
120     ff->setForceFieldFileName(simParams->getForceFieldFileName());
121     }
122    
123     std::string forcefieldFileName;
124     forcefieldFileName = ff->getForceFieldFileName();
125    
126     if (simParams->haveForceFieldVariant()) {
127     //If the force field has variant, the variant force field name will be
128     //Base.variant.frc. For exampel EAM.u6.frc
129    
130     std::string variant = simParams->getForceFieldVariant();
131    
132     std::string::size_type pos = forcefieldFileName.rfind(".frc");
133     variant = "." + variant;
134     if (pos != std::string::npos) {
135     forcefieldFileName.insert(pos, variant);
136     } else {
137     //If the default force field file name does not containt .frc suffix, just append the .variant
138     forcefieldFileName.append(variant);
139     }
140     }
141    
142     ff->parse(forcefieldFileName);
143    
144     //extract the molecule stamps
145     std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
146     compList(stamps, simParams, moleculeStampPairs);
147    
148     //create SimInfo
149     SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
150    
151     //gather parameters (SimCreator only retrieves part of the parameters)
152     gatherParameters(info, mdFileName);
153    
154     //divide the molecules and determine the global index of molecules
155     #ifdef IS_MPI
156     divideMolecules(info);
157     #endif
158    
159     //create the molecules
160     createMolecules(info);
161    
162    
163     //allocate memory for DataStorage(circular reference, need to break it)
164     info->setSnapshotManager(new SimSnapshotManager(info));
165    
166     //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
167     //global index will never change again). Local indices of atoms and rigidbodies are already set by
168     //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
169     setGlobalIndex(info);
170    
171     //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
172     //atoms don't have the global index yet (their global index are all initialized to -1).
173     //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
174     //we can determine the beginning global indices of atoms before they get created.
175     SimInfo::MoleculeIterator mi;
176     Molecule* mol;
177     for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
178     info->addExcludePairs(mol);
179     }
180    
181    
182     //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
183     //eta, chi for NPT integrator)
184     if (loadInitCoords)
185     loadCoordinates(info);
186    
187     return info;
188     }
189    
190     void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
191    
192     //figure out the ouput file names
193     std::string prefix;
194    
195     #ifdef IS_MPI
196    
197     if (worldRank == 0) {
198     #endif // is_mpi
199 tim 2068 Globals * simParams = info->getSimParams();
200 tim 2065 if (simParams->haveFinalConfig()) {
201     prefix = getPrefix(simParams->getFinalConfig());
202     } else {
203     prefix = getPrefix(mdfile);
204     }
205    
206     info->setFinalConfigFileName(prefix + ".eor");
207     info->setDumpFileName(prefix + ".dump");
208     info->setStatFileName(prefix + ".stat");
209    
210     #ifdef IS_MPI
211    
212     }
213    
214     #endif
215    
216     }
217    
218     #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    
242     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    
251     painCave.isFatal = 1;
252     simError();
253     }
254    
255 tim 2068 int seedValue;
256     Globals * simParams = info->getSimParams();
257 tim 2072 MTRand* myRandom; //divide labor does not need Parallel random number generator
258 tim 2068 if (simParams->haveSeed()) {
259     seedValue = simParams->getSeed();
260 tim 2072 myRandom = new MTRand(seedValue);
261 tim 2068 }else {
262 tim 2072 myRandom = new MTRand();
263 tim 2068 }
264 tim 2065
265    
266     a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
267    
268     //initialize atomsPerProc
269     atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
270    
271     if (worldRank == 0) {
272     numerator = info->getNGlobalAtoms();
273     denominator = nProcessors;
274     precast = numerator / denominator;
275     nTarget = (int)(precast + 0.5);
276    
277     for(i = 0; i < nGlobalMols; i++) {
278     done = 0;
279     loops = 0;
280    
281     while (!done) {
282     loops++;
283    
284     // Pick a processor at random
285    
286 tim 2068 which_proc = (int) (myRandom->rand() * nProcessors);
287 tim 2065
288     //get the molecule stamp first
289     int stampId = info->getMoleculeStampId(i);
290     MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
291    
292     // 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    
297     // 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    
301     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    
308     painCave.isFatal = 0;
309     simError();
310    
311     molToProcMap[i] = which_proc;
312     atomsPerProc[which_proc] += add_atoms;
313    
314     done = 1;
315     continue;
316     }
317    
318     // If we can add this molecule to this processor without sending
319     // it above nTarget, then go ahead and do it:
320    
321     if (new_atoms <= nTarget) {
322     molToProcMap[i] = which_proc;
323     atomsPerProc[which_proc] += add_atoms;
324    
325     done = 1;
326     continue;
327     }
328    
329     // 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    
333     // roughly: x = new_atoms - nTarget
334     // Pacc(x) = exp(- a * x)
335     // where a = penalty / (average atoms per molecule)
336    
337     x = (double)(new_atoms - nTarget);
338 tim 2068 y = myRandom->rand();
339 tim 2065
340     if (y < exp(- a * x)) {
341     molToProcMap[i] = which_proc;
342     atomsPerProc[which_proc] += add_atoms;
343    
344     done = 1;
345     continue;
346     } else {
347     continue;
348     }
349     }
350     }
351    
352 tim 2068 delete myRandom;
353    
354 tim 2065 // Spray out this nonsense to all other processors:
355    
356     MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
357     } else {
358    
359     // Listen to your marching orders from processor 0:
360    
361     MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
362     }
363    
364     info->setMolToProcMap(molToProcMap);
365     sprintf(checkPointMsg,
366     "Successfully divided the molecules among the processors.\n");
367     MPIcheckPoint();
368     }
369    
370     #endif
371    
372     void SimCreator::createMolecules(SimInfo *info) {
373     MoleculeCreator molCreator;
374     int stampId;
375    
376     for(int i = 0; i < info->getNGlobalMolecules(); i++) {
377    
378     #ifdef IS_MPI
379    
380     if (info->getMolToProc(i) == worldRank) {
381     #endif
382    
383     stampId = info->getMoleculeStampId(i);
384     Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
385     stampId, i, info->getLocalIndexManager());
386    
387     info->addMolecule(mol);
388    
389     #ifdef IS_MPI
390    
391     }
392    
393     #endif
394    
395     } //end for(int i=0)
396     }
397    
398     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    
406     if (!simParams->haveNMol()) {
407     // we don't have the total number of molecules, so we assume it is
408     // given in each component
409    
410     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    
417     painCave.isFatal = 1;
418     simError();
419     }
420    
421     id = the_components[i]->getType();
422     currentStamp = (stamps->extractMolStamp(id))->getStamp();
423    
424     if (currentStamp == NULL) {
425     sprintf(painCave.errMsg,
426     "SimCreator error: Component \"%s\" was not found in the "
427     "list of declared molecules\n", id);
428    
429     painCave.isFatal = 1;
430     simError();
431     }
432    
433     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    
443     painCave.isFatal = 1;
444     simError();
445     }
446    
447     #ifdef IS_MPI
448    
449     strcpy(checkPointMsg, "Component stamps successfully extracted\n");
450     MPIcheckPoint();
451    
452     #endif // is_mpi
453    
454     }
455    
456     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    
472     beginAtomIndex = 0;
473     beginRigidBodyIndex = 0;
474     beginCutoffGroupIndex = 0;
475    
476     #else
477    
478     int nproc;
479     int myNode;
480    
481     myNode = worldRank;
482     MPI_Comm_size(MPI_COMM_WORLD, &nproc);
483    
484     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    
491     tmpAtomsInProc[myNode] = info->getNAtoms();
492     tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
493     tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
494    
495     //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    
503     beginAtomIndex = 0;
504     beginRigidBodyIndex = 0;
505     beginCutoffGroupIndex = 0;
506    
507     for(int i = 0; i < myNode; i++) {
508     beginAtomIndex += NumAtomsInProc[i];
509     beginRigidBodyIndex += NumRigidBodiesInProc[i];
510     beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
511     }
512    
513     #endif
514    
515     //rigidbody's index begins right after atom's
516     beginRigidBodyIndex += info->getNGlobalAtoms();
517    
518     for(mol = info->beginMolecule(mi); mol != NULL;
519     mol = info->nextMolecule(mi)) {
520    
521     //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    
526     for(rb = mol->beginRigidBody(ri); rb != NULL;
527     rb = mol->nextRigidBody(ri)) {
528     rb->setGlobalIndex(beginRigidBodyIndex++);
529     }
530    
531     //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    
538     //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    
543     for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
544     globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
545     }
546    
547     }
548     }
549    
550     #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    
564     //fill molMembership
565     std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
566    
567     for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
568    
569     for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
570     globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
571     }
572     }
573    
574     #ifdef IS_MPI
575     std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
576    
577     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    
585     }
586    
587     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    
598     DumpReader reader(info, simParams->getInitialConfig());
599     int nframes = reader.getNFrames();
600    
601     if (nframes > 0) {
602     reader.readFrame(nframes - 1);
603     } else {
604     //invalid initial coordinate file
605     sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
606     simParams->getInitialConfig());
607     painCave.isFatal = 1;
608     simError();
609     }
610    
611     //copy the current snapshot to previous snapshot
612     info->getSnapshotManager()->advance();
613     }
614    
615     } //end namespace oopse
616    
617