ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/brains/SimCreator.cpp
(Generate patch)

Comparing trunk/OOPSE-3.0/src/brains/SimCreator.cpp (file contents):
Revision 2087 by gezelter, Tue Mar 8 21:06:49 2005 UTC vs.
Revision 2211 by chrisfen, Thu Apr 21 14:12:19 2005 UTC

# Line 62 | Line 62 | namespace oopse {
62   #endif
63  
64   namespace oopse {
65 <
66 <  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, Globals* simParams){
67 <
65 >  
66 >  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps,
67 >                             Globals* simParams){
68 >    
69   #ifdef IS_MPI
70 <
70 >    
71      if (worldRank == 0) {
72   #endif // is_mpi
73 <
73 >      
74        simParams->initalize();
75        set_interface_stamps(stamps, simParams);
76 <
76 >      
77   #ifdef IS_MPI
78 <
78 >      
79        mpiEventInit();
80 <
80 >      
81   #endif
82 <
82 >      
83        yacc_BASS(mdFileName.c_str());
84 <
84 >      
85   #ifdef IS_MPI
86 <
86 >      
87        throwMPIEvent(NULL);
88      } else {
89        set_interface_stamps(stamps, simParams);
# Line 90 | Line 91 | namespace oopse {
91        MPIcheckPoint();
92        mpiEventLoop();
93      }
94 <
94 >    
95   #endif
96 <
96 >    
97    }
98 <
99 <  SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
98 >  
99 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
100 >                                  bool loadInitCoords) {
101      
102      MakeStamps * stamps = new MakeStamps();
103 <
103 >    
104      Globals * simParams = new Globals();
105 <
105 >    
106      //parse meta-data file
107      parseFile(mdFileName, stamps, simParams);
108 <
108 >    
109      //create the force field
110 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
111 <                                                                         simParams->getForceField());
110 >    ForceField * ff = ForceFieldFactory::getInstance()
111 >      ->createForceField(simParams->getForceField());
112      
113      if (ff == NULL) {
114 <      sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
114 >      sprintf(painCave.errMsg,
115 >              "ForceField Factory can not create %s force field\n",
116                simParams->getForceField());
117        painCave.isFatal = 1;
118        simError();
119      }
120 <
120 >    
121      if (simParams->haveForceFieldFileName()) {
122        ff->setForceFieldFileName(simParams->getForceFieldFileName());
123      }
124      
125      std::string forcefieldFileName;
126      forcefieldFileName = ff->getForceFieldFileName();
127 <
127 >    
128      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 <        
131 >      
132        std::string variant = simParams->getForceFieldVariant();
133 <
133 >      
134        std::string::size_type pos = forcefieldFileName.rfind(".frc");
135        variant = "." + variant;
136        if (pos != std::string::npos) {
# Line 143 | Line 146 | namespace oopse {
146      //extract the molecule stamps
147      std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
148      compList(stamps, simParams, moleculeStampPairs);
149 <
149 >    
150      //create SimInfo
151 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
152 <
151 >    SimInfo * info = new SimInfo(stamps, moleculeStampPairs, ff, simParams);
152 >    
153      //gather parameters (SimCreator only retrieves part of the parameters)
154      gatherParameters(info, mdFileName);
155 <
155 >    
156      //divide the molecules and determine the global index of molecules
157   #ifdef IS_MPI
158      divideMolecules(info);
159   #endif
160 <
160 >    
161      //create the molecules
162      createMolecules(info);
163 <
164 <
163 >    
164 >    
165      //allocate memory for DataStorage(circular reference, need to break it)
166      info->setSnapshotManager(new SimSnapshotManager(info));
167      
# Line 166 | Line 169 | namespace oopse {
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 <
172 >    
173      //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
# Line 177 | Line 180 | namespace oopse {
180        info->addExcludePairs(mol);
181      }
182      
180
181    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
182    //eta, chi for NPT integrator)
183      if (loadInitCoords)
184        loadCoordinates(info);    
185      
186      return info;
187    }
188 <
188 >  
189    void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
190 <
190 >    
191      //figure out the ouput file names
192      std::string prefix;
193 <
193 >    
194   #ifdef IS_MPI
195 <
195 >    
196      if (worldRank == 0) {
197   #endif // is_mpi
198        Globals * simParams = info->getSimParams();
# Line 201 | Line 201 | namespace oopse {
201        } else {
202          prefix = getPrefix(mdfile);
203        }
204 <
204 >      
205        info->setFinalConfigFileName(prefix + ".eor");
206        info->setDumpFileName(prefix + ".dump");
207        info->setStatFileName(prefix + ".stat");
208 <
208 >      info->setRestFileName(prefix + ".zang");
209 >      
210   #ifdef IS_MPI
211 <
211 >      
212      }
213 <
213 >    
214   #endif
215 <
215 >    
216    }
217 <
217 >  
218   #ifdef IS_MPI
219    void SimCreator::divideMolecules(SimInfo *info) {
220      double numerator;
# Line 237 | Line 238 | namespace oopse {
238      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
239      
240      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
241 <
241 >    
242      if (nProcessors > nGlobalMols) {
243        sprintf(painCave.errMsg,
244                "nProcessors (%d) > nMol (%d)\n"
# Line 246 | Line 247 | namespace oopse {
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 <
250 >      
251        painCave.isFatal = 1;
252        simError();
253      }
254 <
254 >    
255      int seedValue;
256      Globals * simParams = info->getSimParams();
257      SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
# Line 260 | Line 261 | namespace oopse {
261      }else {
262        myRandom = new SeqRandNumGen();
263      }  
264 <
265 <
264 >    
265 >    
266      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
267 <
267 >    
268      //initialize atomsPerProc
269      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
270 <
270 >    
271      if (worldRank == 0) {
272        numerator = info->getNGlobalAtoms();
273        denominator = nProcessors;
274        precast = numerator / denominator;
275        nTarget = (int)(precast + 0.5);
276 <
276 >      
277        for(i = 0; i < nGlobalMols; i++) {
278          done = 0;
279          loops = 0;
280 <
280 >        
281          while (!done) {
282            loops++;
283 <
283 >          
284            // Pick a processor at random
285 <
285 >          
286            which_proc = (int) (myRandom->rand() * nProcessors);
287 <
287 >          
288            //get the molecule stamp first
289            int stampId = info->getMoleculeStampId(i);
290            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
291 <
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 <
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 <
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 <
307 >            
308              painCave.isFatal = 0;
309              simError();
310 <
310 >            
311              molToProcMap[i] = which_proc;
312              atomsPerProc[which_proc] += add_atoms;
313 <
313 >            
314              done = 1;
315              continue;
316            }
317 <
317 >          
318            // If we can add this molecule to this processor without sending
319            // it above nTarget, then go ahead and do it:
320 <
320 >          
321            if (new_atoms <= nTarget) {
322              molToProcMap[i] = which_proc;
323              atomsPerProc[which_proc] += add_atoms;
324 <
324 >            
325              done = 1;
326              continue;
327            }
328 <
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 <
332 >          
333            // roughly:  x = new_atoms - nTarget
334            //           Pacc(x) = exp(- a * x)
335            // where a = penalty / (average atoms per molecule)
336 <
336 >          
337            x = (double)(new_atoms - nTarget);
338            y = myRandom->rand();
339 <
339 >          
340            if (y < exp(- a * x)) {
341              molToProcMap[i] = which_proc;
342              atomsPerProc[which_proc] += add_atoms;
343 <
343 >            
344              done = 1;
345              continue;
346            } else {
# Line 347 | Line 348 | namespace oopse {
348            }
349          }
350        }
351 <
351 >      
352        delete myRandom;
353 <        
353 >      
354        // Spray out this nonsense to all other processors:
355 <
355 >      
356        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
357      } else {
358 <
358 >      
359        // Listen to your marching orders from processor 0:
360 <
360 >      
361        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
362      }
363 <
363 >    
364      info->setMolToProcMap(molToProcMap);
365      sprintf(checkPointMsg,
366              "Successfully divided the molecules among the processors.\n");
367      MPIcheckPoint();
368    }
369 <
369 >  
370   #endif
371 <
371 >  
372    void SimCreator::createMolecules(SimInfo *info) {
373      MoleculeCreator molCreator;
374      int stampId;
375 <
375 >    
376      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
377 <
377 >      
378   #ifdef IS_MPI
379 <
379 >      
380        if (info->getMolToProc(i) == worldRank) {
381   #endif
382 <
382 >        
383          stampId = info->getMoleculeStampId(i);
384          Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
385                                                     stampId, i, info->getLocalIndexManager());
386 <
386 >        
387          info->addMolecule(mol);
388 <
388 >        
389   #ifdef IS_MPI
390 <
390 >        
391        }
392 <
392 >      
393   #endif
394 <
394 >      
395      } //end for(int i=0)  
396    }
397 <
397 >  
398    void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
399                              std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
400      int i;
401      char * id;
401    LinkedMolStamp* extractedStamp = NULL;
402      MoleculeStamp * currentStamp;
403      Component** the_components = simParams->getComponents();
404      int n_components = simParams->getNComponents();
405 <
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 <
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 <
416 >          
417            painCave.isFatal = 1;
418            simError();
419          }
420 <      
420 >        
421          id = the_components[i]->getType();
422  
423 <        extractedStamp = stamps->extractMolStamp(id);
424 <        if (extractedStamp == NULL) {
423 >        currentStamp = stamps->getMolStamp(id);
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 <
428 >          
429            painCave.isFatal = 1;
430            simError();
431          }
432 <
433 <        currentStamp = extractedStamp->getStamp();
434 <
435 <
432 >        
433          moleculeStampPairs.push_back(
434                                       std::make_pair(currentStamp, the_components[i]->getNMol()));
435        } //end for (i = 0; i < n_components; i++)
# Line 442 | Line 439 | namespace oopse {
439                " nMols and then give molfractions in the components\n"
440                "\tis not currently supported."
441                " Please give nMol in the components.\n");
442 <
442 >      
443        painCave.isFatal = 1;
444        simError();
445      }
446 <
446 >    
447   #ifdef IS_MPI
448 <
448 >    
449      strcpy(checkPointMsg, "Component stamps successfully extracted\n");
450      MPIcheckPoint();
451 <
451 >    
452   #endif // is_mpi
453 <
453 >    
454    }
455 <
455 >  
456    void SimCreator::setGlobalIndex(SimInfo *info) {
457      SimInfo::MoleculeIterator mi;
458      Molecule::AtomIterator ai;
# Line 471 | Line 468 | namespace oopse {
468      int nGlobalAtoms = info->getNGlobalAtoms();
469      
470   #ifndef IS_MPI
471 <
471 >    
472      beginAtomIndex = 0;
473      beginRigidBodyIndex = 0;
474      beginCutoffGroupIndex = 0;
475 <
475 >    
476   #else
477 <
477 >    
478      int nproc;
479      int myNode;
480 <
480 >    
481      myNode = worldRank;
482      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
483 <
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 <
490 >    
491      tmpAtomsInProc[myNode] = info->getNAtoms();
492      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
493      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
494 <
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);
# Line 502 | Line 499 | namespace oopse {
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 <
502 >    
503      beginAtomIndex = 0;
504      beginRigidBodyIndex = 0;
505      beginCutoffGroupIndex = 0;
506 <
506 >    
507      for(int i = 0; i < myNode; i++) {
508        beginAtomIndex += NumAtomsInProc[i];
509        beginRigidBodyIndex += NumRigidBodiesInProc[i];
510        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
511      }
512 <
512 >    
513   #endif
514 <
514 >    
515      //rigidbody's index begins right after atom's
516      beginRigidBodyIndex += info->getNGlobalAtoms();
517 <
517 >    
518      for(mol = info->beginMolecule(mi); mol != NULL;
519          mol = info->nextMolecule(mi)) {
520 <
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 <
525 >      
526        for(rb = mol->beginRigidBody(ri); rb != NULL;
527            rb = mol->nextRigidBody(ri)) {
528          rb->setGlobalIndex(beginRigidBodyIndex++);
529        }
530 <
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 <
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 <
542 >        
543          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
544            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
545          }
546 <
546 >        
547        }      
548      }
549 <
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
# Line 563 | Line 560 | namespace oopse {
560   #else
561      info->setGlobalGroupMembership(globalGroupMembership);
562   #endif
563 <
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 <
568 >      
569        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
570          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
571        }
572      }
573 <
573 >    
574   #ifdef IS_MPI
575      std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
576 <
576 >    
577      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
578                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
579      
# Line 584 | Line 581 | namespace oopse {
581   #else
582      info->setGlobalMolMembership(globalMolMembership);
583   #endif
584 <
584 >    
585    }
586 <
586 >  
587    void SimCreator::loadCoordinates(SimInfo* info) {
588      Globals* simParams;
589      simParams = info->getSimParams();
# Line 597 | Line 594 | namespace oopse {
594        painCave.isFatal = 1;;
595        simError();
596      }
597 <        
597 >    
598      DumpReader reader(info, simParams->getInitialConfig());
599      int nframes = reader.getNFrames();
600 <
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",
605 >      sprintf(painCave.errMsg,
606 >              "Initial configuration file %s should at least contain one frame\n",
607                simParams->getInitialConfig());
608        painCave.isFatal = 1;
609        simError();
610      }
611 <
611 >    
612      //copy the current snapshot to previous snapshot
613      info->getSnapshotManager()->advance();
614    }
615 <
615 >  
616   } //end namespace oopse
617  
618  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines