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

Comparing trunk/OOPSE-2.0/src/brains/SimCreator.cpp (file contents):
Revision 1981 by tim, Mon Feb 7 19:14:26 2005 UTC vs.
Revision 2087 by gezelter, Tue Mar 8 21:06:49 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 47 | Line 47
47   * @version 1.0
48   */
49  
50 – #include <sprng.h>
51 –
50   #include "brains/MoleculeCreator.hpp"
51   #include "brains/SimCreator.hpp"
52   #include "brains/SimSnapshotManager.hpp"
# Line 57 | Line 55
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/randomSPRNG.hpp"
61 > #include "math/ParallelRandNumGen.hpp"
62   #endif
63  
64   namespace oopse {
65  
66 < void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, Globals* simParams){
66 >  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, Globals* simParams){
67  
68   #ifdef IS_MPI
69  
70      if (worldRank == 0) {
71   #endif // is_mpi
72  
73 <        simParams->initalize();
74 <        set_interface_stamps(stamps, simParams);
73 >      simParams->initalize();
74 >      set_interface_stamps(stamps, simParams);
75  
76   #ifdef IS_MPI
77  
78 <        mpiEventInit();
78 >      mpiEventInit();
79  
80   #endif
81  
82 <        yacc_BASS(mdFileName.c_str());
82 >      yacc_BASS(mdFileName.c_str());
83  
84   #ifdef IS_MPI
85  
86 <        throwMPIEvent(NULL);
86 >      throwMPIEvent(NULL);
87      } else {
88 <        set_interface_stamps(stamps, simParams);
89 <        mpiEventInit();
90 <        MPIcheckPoint();
91 <        mpiEventLoop();
88 >      set_interface_stamps(stamps, simParams);
89 >      mpiEventInit();
90 >      MPIcheckPoint();
91 >      mpiEventLoop();
92      }
93  
94   #endif
95  
96 < }
96 >  }
97  
98 < SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
98 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
99      
100      MakeStamps * stamps = new MakeStamps();
101  
# Line 107 | Line 106 | SimInfo*  SimCreator::createSim(const std::string & md
106  
107      //create the force field
108      ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
109 <                          simParams->getForceField());
109 >                                                                         simParams->getForceField());
110      
111      if (ff == NULL) {
112 <        sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
113 <                simParams->getForceField());
114 <        painCave.isFatal = 1;
115 <        simError();
112 >      sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
113 >              simParams->getForceField());
114 >      painCave.isFatal = 1;
115 >      simError();
116      }
117  
118      if (simParams->haveForceFieldFileName()) {
119 <        ff->setForceFieldFileName(simParams->getForceFieldFileName());
119 >      ff->setForceFieldFileName(simParams->getForceFieldFileName());
120      }
121      
122      std::string forcefieldFileName;
123      forcefieldFileName = ff->getForceFieldFileName();
124  
125      if (simParams->haveForceFieldVariant()) {
126 <        //If the force field has variant, the variant force field name will be
127 <        //Base.variant.frc. For exampel EAM.u6.frc
126 >      //If the force field has variant, the variant force field name will be
127 >      //Base.variant.frc. For exampel EAM.u6.frc
128          
129 <        std::string variant = simParams->getForceFieldVariant();
129 >      std::string variant = simParams->getForceFieldVariant();
130  
131 <        std::string::size_type pos = forcefieldFileName.rfind(".frc");
132 <        variant = "." + variant;
133 <        if (pos != std::string::npos) {
134 <            forcefieldFileName.insert(pos, variant);
135 <        } else {
136 <            //If the default force field file name does not containt .frc suffix, just append the .variant
137 <            forcefieldFileName.append(variant);
138 <        }
131 >      std::string::size_type pos = forcefieldFileName.rfind(".frc");
132 >      variant = "." + variant;
133 >      if (pos != std::string::npos) {
134 >        forcefieldFileName.insert(pos, variant);
135 >      } else {
136 >        //If the default force field file name does not containt .frc suffix, just append the .variant
137 >        forcefieldFileName.append(variant);
138 >      }
139      }
140      
141      ff->parse(forcefieldFileName);
# Line 175 | Line 174 | SimInfo*  SimCreator::createSim(const std::string & md
174      SimInfo::MoleculeIterator mi;
175      Molecule* mol;
176      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
177 <        info->addExcludePairs(mol);
177 >      info->addExcludePairs(mol);
178      }
179      
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);    
184 >      loadCoordinates(info);    
185      
186      return info;
187 < }
187 >  }
188  
189 < void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
189 >  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
190  
192 –    //setup seed for random number generator
193 –    int seedValue;
194 –    Globals * simParams = info->getSimParams();
195 –
196 –    if (simParams->haveSeed()) {
197 –        seedValue = simParams->getSeed();
198 –
199 –        if (seedValue < 100000000 ) {
200 –            sprintf(painCave.errMsg,
201 –                    "Seed for sprng library should contain at least 9 digits\n"
202 –                        "OOPSE will generate a seed for user\n");
203 –
204 –            painCave.isFatal = 0;
205 –            simError();
206 –
207 –            //using seed generated by system instead of invalid seed set by user
208 –
209 – #ifndef IS_MPI
210 –
211 –            seedValue = make_sprng_seed();
212 –
213 – #else
214 –
215 –            if (worldRank == 0) {
216 –                seedValue = make_sprng_seed();
217 –            }
218 –
219 –            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
220 –
221 – #endif
222 –
223 –        } //end if (seedValue /1000000000 == 0)
224 –    } else {
225 –
226 – #ifndef IS_MPI
227 –
228 –        seedValue = make_sprng_seed();
229 –
230 – #else
231 –
232 –        if (worldRank == 0) {
233 –            seedValue = make_sprng_seed();
234 –        }
235 –
236 –        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
237 –
238 – #endif
239 –
240 –    } //end of simParams->haveSeed()
241 –
242 –    info->setSeed(seedValue);
243 –
244 –
191      //figure out the ouput file names
192      std::string prefix;
193  
# Line 249 | Line 195 | void SimCreator::gatherParameters(SimInfo *info, const
195  
196      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  
205 <        if (simParams->haveFinalConfig()) {
206 <            prefix = getPrefix(simParams->getFinalConfig());
207 <        } else {
256 <            prefix = getPrefix(mdfile);
257 <        }
205 >      info->setFinalConfigFileName(prefix + ".eor");
206 >      info->setDumpFileName(prefix + ".dump");
207 >      info->setStatFileName(prefix + ".stat");
208  
259 –        info->setFinalConfigFileName(prefix + ".eor");
260 –        info->setDumpFileName(prefix + ".dump");
261 –        info->setStatFileName(prefix + ".stat");
262 –
209   #ifdef IS_MPI
210  
211      }
212  
213   #endif
214  
215 < }
215 >  }
216  
217   #ifdef IS_MPI
218 < void SimCreator::divideMolecules(SimInfo *info) {
218 >  void SimCreator::divideMolecules(SimInfo *info) {
219      double numerator;
220      double denominator;
221      double precast;
# Line 287 | Line 233 | void SimCreator::divideMolecules(SimInfo *info) {
233      int which_proc;
234      int nProcessors;
235      std::vector<int> atomsPerProc;
290 –    randomSPRNG myRandom(info->getSeed());
236      int nGlobalMols = info->getNGlobalMolecules();
237      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
238      
239      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
240  
241      if (nProcessors > nGlobalMols) {
242 <        sprintf(painCave.errMsg,
243 <                "nProcessors (%d) > nMol (%d)\n"
244 <                    "\tThe number of processors is larger than\n"
245 <                    "\tthe number of molecules.  This will not result in a \n"
246 <                    "\tusable division of atoms for force decomposition.\n"
247 <                    "\tEither try a smaller number of processors, or run the\n"
248 <                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
242 >      sprintf(painCave.errMsg,
243 >              "nProcessors (%d) > nMol (%d)\n"
244 >              "\tThe number of processors is larger than\n"
245 >              "\tthe number of molecules.  This will not result in a \n"
246 >              "\tusable division of atoms for force decomposition.\n"
247 >              "\tEither try a smaller number of processors, or run the\n"
248 >              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
249  
250 <        painCave.isFatal = 1;
251 <        simError();
250 >      painCave.isFatal = 1;
251 >      simError();
252      }
253  
254 +    int seedValue;
255 +    Globals * simParams = info->getSimParams();
256 +    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
257 +    if (simParams->haveSeed()) {
258 +      seedValue = simParams->getSeed();
259 +      myRandom = new SeqRandNumGen(seedValue);
260 +    }else {
261 +      myRandom = new SeqRandNumGen();
262 +    }  
263 +
264 +
265      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
266  
267      //initialize atomsPerProc
268      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
269  
270      if (worldRank == 0) {
271 <        numerator = info->getNGlobalAtoms();
272 <        denominator = nProcessors;
273 <        precast = numerator / denominator;
274 <        nTarget = (int)(precast + 0.5);
271 >      numerator = info->getNGlobalAtoms();
272 >      denominator = nProcessors;
273 >      precast = numerator / denominator;
274 >      nTarget = (int)(precast + 0.5);
275  
276 <        for(i = 0; i < nGlobalMols; i++) {
277 <            done = 0;
278 <            loops = 0;
276 >      for(i = 0; i < nGlobalMols; i++) {
277 >        done = 0;
278 >        loops = 0;
279  
280 <            while (!done) {
281 <                loops++;
280 >        while (!done) {
281 >          loops++;
282  
283 <                // Pick a processor at random
283 >          // Pick a processor at random
284  
285 <                which_proc = (int) (myRandom.getRandom() * nProcessors);
285 >          which_proc = (int) (myRandom->rand() * nProcessors);
286  
287 <                //get the molecule stamp first
288 <                int stampId = info->getMoleculeStampId(i);
289 <                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
287 >          //get the molecule stamp first
288 >          int stampId = info->getMoleculeStampId(i);
289 >          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
290  
291 <                // How many atoms does this processor have so far?
292 <                old_atoms = atomsPerProc[which_proc];
293 <                add_atoms = moleculeStamp->getNAtoms();
294 <                new_atoms = old_atoms + add_atoms;
291 >          // How many atoms does this processor have so far?
292 >          old_atoms = atomsPerProc[which_proc];
293 >          add_atoms = moleculeStamp->getNAtoms();
294 >          new_atoms = old_atoms + add_atoms;
295  
296 <                // If we've been through this loop too many times, we need
297 <                // to just give up and assign the molecule to this processor
298 <                // and be done with it.
296 >          // If we've been through this loop too many times, we need
297 >          // to just give up and assign the molecule to this processor
298 >          // and be done with it.
299  
300 <                if (loops > 100) {
301 <                    sprintf(painCave.errMsg,
302 <                            "I've tried 100 times to assign molecule %d to a "
303 <                                " processor, but can't find a good spot.\n"
304 <                                "I'm assigning it at random to processor %d.\n",
305 <                            i, which_proc);
300 >          if (loops > 100) {
301 >            sprintf(painCave.errMsg,
302 >                    "I've tried 100 times to assign molecule %d to a "
303 >                    " processor, but can't find a good spot.\n"
304 >                    "I'm assigning it at random to processor %d.\n",
305 >                    i, which_proc);
306  
307 <                    painCave.isFatal = 0;
308 <                    simError();
307 >            painCave.isFatal = 0;
308 >            simError();
309  
310 <                    molToProcMap[i] = which_proc;
311 <                    atomsPerProc[which_proc] += add_atoms;
310 >            molToProcMap[i] = which_proc;
311 >            atomsPerProc[which_proc] += add_atoms;
312  
313 <                    done = 1;
314 <                    continue;
315 <                }
313 >            done = 1;
314 >            continue;
315 >          }
316  
317 <                // If we can add this molecule to this processor without sending
318 <                // it above nTarget, then go ahead and do it:
317 >          // If we can add this molecule to this processor without sending
318 >          // it above nTarget, then go ahead and do it:
319  
320 <                if (new_atoms <= nTarget) {
321 <                    molToProcMap[i] = which_proc;
322 <                    atomsPerProc[which_proc] += add_atoms;
320 >          if (new_atoms <= nTarget) {
321 >            molToProcMap[i] = which_proc;
322 >            atomsPerProc[which_proc] += add_atoms;
323  
324 <                    done = 1;
325 <                    continue;
326 <                }
371 <
372 <                // The only situation left is when new_atoms > nTarget.  We
373 <                // want to accept this with some probability that dies off the
374 <                // farther we are from nTarget
324 >            done = 1;
325 >            continue;
326 >          }
327  
328 <                // roughly:  x = new_atoms - nTarget
329 <                //           Pacc(x) = exp(- a * x)
330 <                // where a = penalty / (average atoms per molecule)
328 >          // The only situation left is when new_atoms > nTarget.  We
329 >          // want to accept this with some probability that dies off the
330 >          // farther we are from nTarget
331  
332 <                x = (double)(new_atoms - nTarget);
333 <                y = myRandom.getRandom();
332 >          // roughly:  x = new_atoms - nTarget
333 >          //           Pacc(x) = exp(- a * x)
334 >          // where a = penalty / (average atoms per molecule)
335  
336 <                if (y < exp(- a * x)) {
337 <                    molToProcMap[i] = which_proc;
385 <                    atomsPerProc[which_proc] += add_atoms;
336 >          x = (double)(new_atoms - nTarget);
337 >          y = myRandom->rand();
338  
339 <                    done = 1;
340 <                    continue;
341 <                } else {
342 <                    continue;
343 <                }
344 <            }
339 >          if (y < exp(- a * x)) {
340 >            molToProcMap[i] = which_proc;
341 >            atomsPerProc[which_proc] += add_atoms;
342 >
343 >            done = 1;
344 >            continue;
345 >          } else {
346 >            continue;
347 >          }
348          }
349 +      }
350  
351 <        // Spray out this nonsense to all other processors:
351 >      delete myRandom;
352 >        
353 >      // Spray out this nonsense to all other processors:
354  
355 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
355 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
356      } else {
357  
358 <        // Listen to your marching orders from processor 0:
358 >      // Listen to your marching orders from processor 0:
359  
360 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
360 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
361      }
362  
363      info->setMolToProcMap(molToProcMap);
364      sprintf(checkPointMsg,
365              "Successfully divided the molecules among the processors.\n");
366      MPIcheckPoint();
367 < }
367 >  }
368  
369   #endif
370  
371 < void SimCreator::createMolecules(SimInfo *info) {
371 >  void SimCreator::createMolecules(SimInfo *info) {
372      MoleculeCreator molCreator;
373      int stampId;
374  
# Line 418 | Line 376 | void SimCreator::createMolecules(SimInfo *info) {
376  
377   #ifdef IS_MPI
378  
379 <        if (info->getMolToProc(i) == worldRank) {
379 >      if (info->getMolToProc(i) == worldRank) {
380   #endif
381  
382 <            stampId = info->getMoleculeStampId(i);
383 <            Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
384 <                                                                                    stampId, i, info->getLocalIndexManager());
382 >        stampId = info->getMoleculeStampId(i);
383 >        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
384 >                                                   stampId, i, info->getLocalIndexManager());
385  
386 <            info->addMolecule(mol);
386 >        info->addMolecule(mol);
387  
388   #ifdef IS_MPI
389  
390 <        }
390 >      }
391  
392   #endif
393  
394      } //end for(int i=0)  
395 < }
395 >  }
396  
397 < void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
398 <                        std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
397 >  void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
398 >                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
399      int i;
400      char * id;
401 +    LinkedMolStamp* extractedStamp = NULL;
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
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");
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 <            }
417 >          painCave.isFatal = 1;
418 >          simError();
419 >        }
420 >      
421 >        id = the_components[i]->getType();
422  
423 <            id = the_components[i]->getType();
424 <            currentStamp = (stamps->extractMolStamp(id))->getStamp();
423 >        extractedStamp = stamps->extractMolStamp(id);
424 >        if (extractedStamp == NULL) {
425 >          sprintf(painCave.errMsg,
426 >                  "SimCreator error: Component \"%s\" was not found in the "
427 >                  "list of declared molecules\n", id);
428  
429 <            if (currentStamp == NULL) {
430 <                sprintf(painCave.errMsg,
431 <                        "SimCreator error: Component \"%s\" was not found in the "
468 <                            "list of declared molecules\n", id);
429 >          painCave.isFatal = 1;
430 >          simError();
431 >        }
432  
433 <                painCave.isFatal = 1;
471 <                simError();
472 <            }
433 >        currentStamp = extractedStamp->getStamp();
434  
435 <            moleculeStampPairs.push_back(
436 <                std::make_pair(currentStamp, the_components[i]->getNMol()));
437 <        } //end for (i = 0; i < n_components; i++)
435 >
436 >        moleculeStampPairs.push_back(
437 >                                     std::make_pair(currentStamp, the_components[i]->getNMol()));
438 >      } //end for (i = 0; i < n_components; i++)
439      } else {
440 <        sprintf(painCave.errMsg, "SimSetup error.\n"
441 <                                     "\tSorry, the ability to specify total"
442 <                                     " nMols and then give molfractions in the components\n"
443 <                                     "\tis not currently supported."
444 <                                     " Please give nMol in the components.\n");
440 >      sprintf(painCave.errMsg, "SimSetup error.\n"
441 >              "\tSorry, the ability to specify total"
442 >              " nMols and then give molfractions in the components\n"
443 >              "\tis not currently supported."
444 >              " Please give nMol in the components.\n");
445  
446 <        painCave.isFatal = 1;
447 <        simError();
446 >      painCave.isFatal = 1;
447 >      simError();
448      }
449  
450   #ifdef IS_MPI
# Line 492 | Line 454 | void SimCreator::compList(MakeStamps *stamps, Globals*
454  
455   #endif // is_mpi
456  
457 < }
457 >  }
458  
459 < void SimCreator::setGlobalIndex(SimInfo *info) {
459 >  void SimCreator::setGlobalIndex(SimInfo *info) {
460      SimInfo::MoleculeIterator mi;
461      Molecule::AtomIterator ai;
462      Molecule::RigidBodyIterator ri;
# Line 546 | Line 508 | void SimCreator::setGlobalIndex(SimInfo *info) {
508      beginCutoffGroupIndex = 0;
509  
510      for(int i = 0; i < myNode; i++) {
511 <        beginAtomIndex += NumAtomsInProc[i];
512 <        beginRigidBodyIndex += NumRigidBodiesInProc[i];
513 <        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
511 >      beginAtomIndex += NumAtomsInProc[i];
512 >      beginRigidBodyIndex += NumRigidBodiesInProc[i];
513 >      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
514      }
515  
516   #endif
# Line 559 | Line 521 | void SimCreator::setGlobalIndex(SimInfo *info) {
521      for(mol = info->beginMolecule(mi); mol != NULL;
522          mol = info->nextMolecule(mi)) {
523  
524 <        //local index(index in DataStorge) of atom is important
525 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
526 <            atom->setGlobalIndex(beginAtomIndex++);
527 <        }
524 >      //local index(index in DataStorge) of atom is important
525 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
526 >        atom->setGlobalIndex(beginAtomIndex++);
527 >      }
528  
529 <        for(rb = mol->beginRigidBody(ri); rb != NULL;
530 <            rb = mol->nextRigidBody(ri)) {
531 <            rb->setGlobalIndex(beginRigidBodyIndex++);
532 <        }
529 >      for(rb = mol->beginRigidBody(ri); rb != NULL;
530 >          rb = mol->nextRigidBody(ri)) {
531 >        rb->setGlobalIndex(beginRigidBodyIndex++);
532 >      }
533  
534 <        //local index of cutoff group is trivial, it only depends on the order of travesing
535 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
536 <            cg = mol->nextCutoffGroup(ci)) {
537 <            cg->setGlobalIndex(beginCutoffGroupIndex++);
538 <        }
534 >      //local index of cutoff group is trivial, it only depends on the order of travesing
535 >      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
536 >          cg = mol->nextCutoffGroup(ci)) {
537 >        cg->setGlobalIndex(beginCutoffGroupIndex++);
538 >      }
539      }
540  
541      //fill globalGroupMembership
542      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
543      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
544 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
544 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
545  
546 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
547 <                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
548 <            }
546 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
547 >          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
548 >        }
549  
550 <        }      
550 >      }      
551      }
552  
553   #ifdef IS_MPI    
# Line 597 | Line 559 | void SimCreator::setGlobalIndex(SimInfo *info) {
559      std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
560      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
561                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
562 <     info->setGlobalGroupMembership(tmpGroupMembership);
562 >    info->setGlobalGroupMembership(tmpGroupMembership);
563   #else
564      info->setGlobalGroupMembership(globalGroupMembership);
565   #endif
# Line 607 | Line 569 | void SimCreator::setGlobalIndex(SimInfo *info) {
569      
570      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
571  
572 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
573 <            globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
574 <        }
572 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
573 >        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
574 >      }
575      }
576  
577   #ifdef IS_MPI
# Line 623 | Line 585 | void SimCreator::setGlobalIndex(SimInfo *info) {
585      info->setGlobalMolMembership(globalMolMembership);
586   #endif
587  
588 < }
588 >  }
589  
590 < void SimCreator::loadCoordinates(SimInfo* info) {
590 >  void SimCreator::loadCoordinates(SimInfo* info) {
591      Globals* simParams;
592      simParams = info->getSimParams();
593      
594      if (!simParams->haveInitialConfig()) {
595 <        sprintf(painCave.errMsg,
596 <                "Cannot intialize a simulation without an initial configuration file.\n");
597 <        painCave.isFatal = 1;;
598 <        simError();
595 >      sprintf(painCave.errMsg,
596 >              "Cannot intialize a simulation without an initial configuration file.\n");
597 >      painCave.isFatal = 1;;
598 >      simError();
599      }
600          
601      DumpReader reader(info, simParams->getInitialConfig());
602      int nframes = reader.getNFrames();
603  
604      if (nframes > 0) {
605 <        reader.readFrame(nframes - 1);
605 >      reader.readFrame(nframes - 1);
606      } else {
607 <        //invalid initial coordinate file
608 <        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
609 <                simParams->getInitialConfig());
610 <        painCave.isFatal = 1;
611 <        simError();
607 >      //invalid initial coordinate file
608 >      sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
609 >              simParams->getInitialConfig());
610 >      painCave.isFatal = 1;
611 >      simError();
612      }
613  
614      //copy the current snapshot to previous snapshot
615      info->getSnapshotManager()->advance();
616 < }
616 >  }
617  
618   } //end namespace oopse
619  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines