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

Comparing branches/new_design/OOPSE-3.0/src/brains/SimCreator.cpp (file contents):
Revision 1732 by tim, Thu Nov 11 16:41:58 2004 UTC vs.
Revision 1733 by tim, Fri Nov 12 06:19:04 2004 UTC

# Line 32 | Line 32 | namespace oopse {
32   */
33  
34   #include "brains/SimCreator.hpp"
35 +
36   namespace oopse {
37  
38 < void SimSetup::parseFile(char* mdfile,  MakeStamps* stamps, Globals* globals){
38 > void SimSetup::parseFile(const std::string& mdFileName, MakeStamps *stamps, Globals *globals) {
39 >
40   #ifdef IS_MPI
41 <  if (worldRank == 0){
41 >
42 >    if (worldRank == 0) {
43   #endif // is_mpi
44  
45 +        globals->initalize();
46 +        set_interface_stamps(stamps, globals);
47  
43 –    globals->initalize();
44 –    set_interface_stamps(stamps, globals);
45 –
48   #ifdef IS_MPI
49 <    mpiEventInit();
49 >
50 >        mpiEventInit();
51 >
52   #endif
53  
54 <    yacc_BASS(mdfile);
54 >        yacc_BASS(mdFileName.c_str());
55  
56   #ifdef IS_MPI
57 <    throwMPIEvent(NULL);
58 <  }
59 <  else{
60 <  set_interface_stamps(stamps, globals);
61 <  mpiEventInit();
62 <  MPIcheckPoint();
63 <  mpiEventLoop();
64 <  }
57 >
58 >        throwMPIEvent(NULL);
59 >    } else {
60 >        set_interface_stamps(stamps, globals);
61 >        mpiEventInit();
62 >        MPIcheckPoint();
63 >        mpiEventLoop();
64 >    }
65 >
66   #endif
67  
68   }
69  
70 + SimInfo*  SimCreator::createSim(const std::string & mdFileName) {
71  
72 < SimInfo* SimCreator::createSim(const std::string& mdfile) {
67 <    MakeStamps* stamps = new MakeStamps();
72 >    mdFileName_ = mdFileName;
73      
74 <    Globals* globals = new Globals();
74 >    MakeStamps * stamps = new MakeStamps();
75  
76 +    Globals * globals = new Globals();
77 +
78      //parse meta-data file
79 <    parseFile(mdfile, stamps, globals);
79 >    parseFile(mdFileName_, stamps, globals);
80  
81      //create the force field
82 <    ForceFiled* ff = ForceFieldFactory::getInstance()->createObject(globals->getForceField());
83 <    
82 >    ForceFiled * ff = ForceFieldFactory::getInstance()->createObject(
83 >                          globals->getForceField());
84 >
85 >    if (ff == NULL) {
86 >        sprintf(painCave.errMsg, "ForceFiled Factory can not create %s force field\n",
87 >                globals->getForceField());
88 >        painCave.isFatal = 1;
89 >        simError();
90 >    }
91 >
92 >    //extract the molecule stamps
93 >    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
94 >    compList(stamps, globals, moleculeStampPairs);
95 >
96      //create SimInfo
97 <    SimInfo* info = new SimInfo();
79 <    info->setGlobals(globals);
80 <    info->setForceField(ff);
97 >    SimInfo * info = new SimInfo(moleculeStampPairs, ff, globals);
98  
99 <    //extract the molecule stamps and add them into SimInfo
83 <    compList(stamps, info);
84 <    
85 <    //gather parameters (SimCreator only retrieves the parameters which will be used to create
86 <    // the simulation)
99 >    //gather parameters (SimCreator only retrieves part of the parameters)
100      gatherParameters(info);
101 <    
101 >
102      //divide the molecules and determine the global index of molecules
103 <    
103 >    divideMolecules(info);
104 >
105      //create the molecules
106      createMolecules(info);
107  
108      //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
109      //global index will never change again). Local indices of atoms and rigidbodies are already set by
110 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
110 >    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
111      setGlobalIndices(info);
112 <    
112 >
113      //allocate memory for DataStorage(circular reference, need to break it)
114 <    info->setSnapshotManager(new SimSnapshotManager(info);
101 <    
102 <    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
103 <    //eta, chi for NPT integrator)
104 <    DumpReader reader(info->getInitFilename());
105 <    int nframes = reader->getNframes();
114 >    info->setSnapshotManager(new SimSnapshotManager(info));
115  
116 <    if (nframes > 0) {
117 <        reader.readFrame(info, nframes - 1);
109 <    } else {
110 <        //invalid initial coordinate file
111 <        
112 <    }
116 >    //initialize fortran -- setup the cutoff
117 >    initFortran(info);
118  
119 <    
120 <    //initialize fortran    
121 <    
119 >    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
120 >    //eta, chi for NPT integrator)
121 >    loadCoordinates(info);    
122      return info;
123 < }    
123 > }
124  
125 < void SimCreator::gatherParameters(SimInfo* info) {
121 <    //model->addProperty(new StringGenericData("Ensemble", globals->getForceFiled()));
122 <    //model->addProperty(new DoubleGenericData("dt"), globals->getDt());
125 > void SimCreator::gatherParameters(SimInfo *info) {
126  
127      //setup seed for random number generator
128      int seedValue;
129 <    Globals* globals = info->getGlobals();
130 <    
131 <    if (globals->haveSeed()){
129 >    Globals * globals = info->getGlobals();
130 >
131 >    if (globals->haveSeed()) {
132          seedValue = globals->getSeed();
133  
134 <        if(seedValue / 1000000000 == 0){
134 >        if (seedValue / 1000000000 == 0) {
135              sprintf(painCave.errMsg,
136 <            "Seed for sprng library should contain at least 9 digits\n"
137 <            "OOPSE will generate a seed for user\n");
136 >                    "Seed for sprng library should contain at least 9 digits\n"
137 >                        "OOPSE will generate a seed for user\n");
138 >
139              painCave.isFatal = 0;
140              simError();
141  
142 <    //using seed generated by system instead of invalid seed set by user
142 >            //using seed generated by system instead of invalid seed set by user
143 >
144   #ifndef IS_MPI
145 +
146              seedValue = make_sprng_seed();
147 < #else
148 <            if (worldRank == 0){
147 >
148 > #else
149 >
150 >            if (worldRank == 0) {
151                  seedValue = make_sprng_seed();
152              }
153 <            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
154 < #endif      
153 >
154 >            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
155 >
156 > #endif
157 >
158          } //end if (seedValue /1000000000 == 0)
159 <    } else{
160 <    
159 >    } else {
160 >
161   #ifndef IS_MPI
162 +
163          seedValue = make_sprng_seed();
164 < #else
165 <        if (worldRank == 0){
164 >
165 > #else
166 >
167 >        if (worldRank == 0) {
168              seedValue = make_sprng_seed();
169          }
170 <        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
170 >
171 >        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
172 >
173   #endif
174 <    }//end of globals->haveSeed()
174 >
175 >    } //end of globals->haveSeed()
176 >
177      info->setSeed(seedValue);
178  
179 <    //
179 >
180 >    //figure out the ouput file names
181      std::string prefix;
182 +
183   #ifdef IS_MPI
184 <    if (worldRank == 0){
184 >
185 >    if (worldRank == 0) {
186   #endif // is_mpi
187 <      
188 <        if(globals->haveFinalConfig()) {
189 <            prefix = getPrefix(globals->getFinalConfig());  
187 >
188 >        if (globals->haveFinalConfig()) {
189 >            prefix = StringUtils::getPrefix(globals->getFinalConfig());
190          } else {
191 <            prefix = getPrefix(mdfile_);
191 >            prefix = StringUtils::getPrefix(mdfile_);
192          }
193 <        
194 <        info->setFinalConfFileName(prefix + ".eor");
193 >
194 >        info->setFinalConfigFileName(prefix + ".eor");
195          info->setDumpFileName(prefix + ".dump");
196          info->setStatFileName(prefix + ".stat");
197  
198   #ifdef IS_MPI
199 +
200      }
179 – #endif    
201  
202 + #endif
203  
182 –
204   }
205  
206   #ifdef IS_MPI
186 – void SimCreator::divideMolecules(SimInfo* info){
187 –    int nComponents;
188 –    MoleculeStamp ** compStamps;
189 –    randomSPRNG * myRandom;
190 –    int * componentsNmol;
191 –    int * AtomsPerProc;
207  
208 + void SimCreator::divideMolecules(SimInfo *info) {
209      double numerator;
210      double denominator;
211      double precast;
# Line 199 | Line 215 | void SimCreator::divideMolecules(SimInfo* info){
215      int old_atoms;
216      int add_atoms;
217      int new_atoms;
202 –
218      int nTarget;
204 –    int molIndex;
205 –    int atomIndex;
219      int done;
220      int i;
221      int j;
222      int loops;
223      int which_proc;
224 <    int nmol_global,
225 <    nmol_local;
226 <    int natoms_global,
224 >    int nProcessors;
225 >    int nGlobalMols;
226 >    std::vector < int > atomsPerProc;
227 >    randomSPRNG myRandom(info->getSeed());
228 >    int * molToProcMap = info->getMolToProcMapPointer();
229  
230 <    int baseSeed = info->getSeed();
231 <    CutoffGroupStamp * cg;
230 >    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
231 >    nGlobalMols = info->getNGlobalMolecules();
232  
233 <    nComponents = info->nComponents;
219 <    compStamps = info->compStamps;
220 <    componentsNmol = info->componentsNmol;
221 <    AtomsPerProc = new int[parallelData->nProcessors];
222 <
223 <    parallelData->nMolGlobal = info->n_mol;
224 <
225 <    if (parallelData->nProcessors > parallelData->nMolGlobal) {
233 >    if (nProcessors > nGlobalMols) {
234          sprintf(painCave.errMsg,
235                  "nProcessors (%d) > nMol (%d)\n"
236                      "\tThe number of processors is larger than\n"
237                      "\tthe number of molecules.  This will not result in a \n"
238                      "\tusable division of atoms for force decomposition.\n"
239                      "\tEither try a smaller number of processors, or run the\n"
240 <                    "\tsingle-processor version of OOPSE.\n",
233 <                parallelData->nProcessors,
234 <                parallelData->nMolGlobal);
240 >                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
241  
242          painCave.isFatal = 1;
243          simError();
244      }
245  
246 <    myRandom = new randomSPRNG(baseSeed);
246 >    a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
247  
248 <    a = 3.0 * parallelData->nMolGlobal / parallelData->nAtomsGlobal;
248 >    //initialize atomsPerProc
249 >    atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
250  
251 <    // Initialize things that we'll send out later:
245 <    for( i = 0; i < parallelData->nProcessors; i++ ) {
246 <        AtomsPerProc[i] = 0;
247 <    }
248 <
249 <    for( i = 0; i < parallelData->nMolGlobal; i++ ) {
251 >    for(i = 0; i < nGlobalMols; i++) {
252          // default to an error condition:
253 <        MolToProcMap[i] = -1;
253 >        molToProcMap[i] = -1;
254      }
255  
256 <    if (parallelData->myNode == 0) {
257 <        numerator = info->n_atoms;
258 <        denominator = parallelData->nProcessors;
256 >    if (worldRank == 0) {
257 >        numerator = info->getNGlobalAtoms();
258 >        denominator = nProcessors;
259          precast = numerator / denominator;
260          nTarget = (int)(precast + 0.5);
261  
262 <        // Build the array of molecule component types first
261 <        molIndex = 0;
262 <
263 <        for( i = 0; i < nComponents; i++ ) {
264 <            for( j = 0; j < componentsNmol[i]; j++ ) {
265 <                molIndex++;
266 <            }
267 <        }
268 <
269 <        for( i = 0; i < molIndex; i++ ) {
262 >        for(i = 0; i < nGlobalMols; i++) {
263              done = 0;
264              loops = 0;
265  
# Line 275 | Line 268 | void SimCreator::divideMolecules(SimInfo* info){
268  
269                  // Pick a processor at random
270  
271 <                which_proc
271 >                which_proc = (int) (myRandom.getRandom() * nProcessors);
272  
273 <                = (int) (myRandom->getRandom() * parallelData->nProcessors);
273 >                //get the molecule stamp first
274 >                int stampId = info->getMoleculeStampId(i);
275 >                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
276  
277 <                // How many atoms does this processor have?
278 <
279 <                old_atoms = AtomsPerProc[which_proc];
285 <                add_atoms = compStamps[MolComponentType[i]]->getNAtoms();
277 >                // How many atoms does this processor have so far?
278 >                old_atoms = atomsPerProc[which_proc];
279 >                add_atoms = moleculeStamp->getNAtoms();
280                  new_atoms = old_atoms + add_atoms;
281  
282                  // If we've been through this loop too many times, we need
# Line 290 | Line 284 | void SimCreator::divideMolecules(SimInfo* info){
284                  // and be done with it.
285  
286                  if (loops > 100) {
287 <                    sprintf(
288 <                        painCave.errMsg,
289 <                        "I've tried 100 times to assign molecule %d to a "
290 <                            " processor, but can't find a good spot.\n"
291 <                            "I'm assigning it at random to processor %d.\n",
298 <                        i,
299 <                        which_proc);
287 >                    sprintf(painCave.errMsg,
288 >                            "I've tried 100 times to assign molecule %d to a "
289 >                                " processor, but can't find a good spot.\n"
290 >                                "I'm assigning it at random to processor %d.\n",
291 >                            i, which_proc);
292  
293                      painCave.isFatal = 0;
294                      simError();
295  
296 <                    MolToProcMap[i] = which_proc;
297 <                    AtomsPerProc[which_proc] += add_atoms;
296 >                    molToProcMap[i] = which_proc;
297 >                    atomsPerProc[which_proc] += add_atoms;
298  
299                      done = 1;
300                      continue;
# Line 312 | Line 304 | void SimCreator::divideMolecules(SimInfo* info){
304                  // it above nTarget, then go ahead and do it:
305  
306                  if (new_atoms <= nTarget) {
307 <                    MolToProcMap[i] = which_proc;
308 <                    AtomsPerProc[which_proc] += add_atoms;
307 >                    molToProcMap[i] = which_proc;
308 >                    atomsPerProc[which_proc] += add_atoms;
309  
310                      done = 1;
311                      continue;
# Line 328 | Line 320 | void SimCreator::divideMolecules(SimInfo* info){
320                  // where a = penalty / (average atoms per molecule)
321  
322                  x = (double)(new_atoms - nTarget);
323 <                y = myRandom->getRandom();
323 >                y = myRandom.getRandom();
324  
325                  if (y < exp(- a * x)) {
326 <                    MolToProcMap[i] = which_proc;
327 <                    AtomsPerProc[which_proc] += add_atoms;
326 >                    molToProcMap[i] = which_proc;
327 >                    atomsPerProc[which_proc] += add_atoms;
328  
329                      done = 1;
330                      continue;
331 <                } else { continue; }
331 >                } else {
332 >                    continue;
333 >                }
334              }
335          }
336  
337          // Spray out this nonsense to all other processors:
338  
339 <        MPI_Bcast(MolToProcMap, parallelData->nMolGlobal, MPI_INT, 0,
346 <                  MPI_COMM_WORLD);
339 >        MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
340      } else {
341  
342 <    // Listen to your marching orders from processor 0:
350 <
351 <        MPI_Bcast(MolToProcMap, parallelData->nMolGlobal, MPI_INT, 0,
352 <                  MPI_COMM_WORLD); }
342 >        // Listen to your marching orders from processor 0:
343  
344 <    // Let's all check for sanity:
355 <
356 <    nmol_local = 0;
357 <
358 <    for( i = 0; i < parallelData->nMolGlobal; i++ ) {
359 <        if (MolToProcMap[i] == parallelData->myNode) {
360 <            nmol_local++;
361 <        }
344 >        MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
345      }
346  
364 –    MPI_Allreduce(&nmol_local, &nmol_global, 1, MPI_INT, MPI_SUM,
365 –                  MPI_COMM_WORLD);
366 –
367 –    if (nmol_global != info->n_mol) {
368 –        sprintf(painCave.errMsg,
369 –                "The sum of all nmol_local, %d, did not equal the "
370 –                    "total number of molecules, %d.\n",
371 –                nmol_global,
372 –                info->n_mol);
373 –
374 –        painCave.isFatal = 1;
375 –        simError();
376 –    }
377 –
347      sprintf(checkPointMsg,
348              "Successfully divided the molecules among the processors.\n");
349      MPIcheckPoint();
381 –
350   }
351 +
352   #endif
353  
354 < Molecule* SimCreator::createMolecules(SimInfo* info) {
354 > void SimCreator::createMolecules(SimInfo *info) {
355      MoleculeCreator molCreator;
356      int stampId;
357  
358 <    
390 <    for (int i = 0; i < info->getNGlobalMolecules(); i++){
391 <        if (info->getMolToProc(i) == worldRank){
358 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
359  
360 + #ifdef IS_MPI
361 +
362 +        if (info->getMolToProc(i) == worldRank) {
363 + #endif
364 +
365              stampId = info->getMoleculeStampId(i);
366 <            Molecule* mol = molCreator.createMolecule(info->getForceField(),
367 <                                                                                info->getMoleculeStamp(stampId), stampId, i);
368 <                        
366 >            Molecule * mol = molCreator.createMolecule(info->getForceField(),
367 >                                                info->getMoleculeStamp(stampId), stampId, i);
368 >
369              info->addMolecule(mol);
370 +
371 + #ifdef IS_MPI
372 +
373          }
374 <    }        
375 <      
374 >
375 > #endif
376 >
377 >    } //end for(int i=0)  
378   }
379  
380 < void SimSetup::compList(MakeStamps* stamps,SimInfo* info) {
380 > void SimSetup::compList(MakeStamps *stamps, Globals* globals,
381 >                        std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
382      int i;
383 <    char* id;
384 <    MoleculeStamp* currentStamp;
385 <    Component* the_components = info->getGlobals()->getComponents();
386 <    int n_components = info->getGlobals()->getNComponents();
383 >    char * id;
384 >    MoleculeStamp * currentStamp;
385 >    Component * the_components = globals->getComponents();
386 >    int n_components = globals->getNComponents();
387  
388 <    if (!globals->haveNMol()){
388 >    if (!globals->haveNMol()) {
389          // we don't have the total number of molecules, so we assume it is
390          // given in each component
391  
392 <        for (i = 0; i < n_components; i++){
393 <            if (!the_components[i]->haveNMol()){
392 >        for(i = 0; i < n_components; i++) {
393 >            if (!the_components[i]->haveNMol()) {
394                  // we have a problem
395                  sprintf(painCave.errMsg,
396 <                    "SimSetup Error. No global NMol or component NMol given.\n"
397 <                    "\tCannot calculate the number of atoms.\n");
396 >                        "SimSetup Error. No global NMol or component NMol given.\n"
397 >                            "\tCannot calculate the number of atoms.\n");
398 >
399                  painCave.isFatal = 1;
400                  simError();
401              }
402  
403              id = the_components[i]->getType();
404              currentStamp = stamps->extractMolStamp(id);
405 <            if (currentStamp == NULL){
405 >
406 >            if (currentStamp == NULL) {
407                  sprintf(painCave.errMsg,
408 <                    "SimSetup error: Component \"%s\" was not found in the "
409 <                    "list of declared molecules\n",
410 <                    id);
408 >                        "SimSetup error: Component \"%s\" was not found in the "
409 >                            "list of declared molecules\n", id);
410 >
411                  painCave.isFatal = 1;
412                  simError();
413 <            }      
413 >            }
414  
415 <            info->addMoleculeStamp(currentStamp, the_components[i]->getNMol);
416 <          
415 >            moleculeStampPairs.push_back(
416 >                make_pair(currentStamp, the_components[i]->getNMol));
417          } //end for (i = 0; i < n_components; i++)
418 <        
419 <    } else{
420 <        sprintf(painCave.errMsg,
421 <                "SimSetup error.\n"
422 <                "\tSorry, the ability to specify total"
423 <                " nMols and then give molfractions in the components\n"
424 <                "\tis not currently supported."
445 <                " Please give nMol in the components.\n");
418 >    } else {
419 >        sprintf(painCave.errMsg, "SimSetup error.\n"
420 >                                     "\tSorry, the ability to specify total"
421 >                                     " nMols and then give molfractions in the components\n"
422 >                                     "\tis not currently supported."
423 >                                     " Please give nMol in the components.\n");
424 >
425          painCave.isFatal = 1;
426          simError();
427      }
428 <  
428 >
429   #ifdef IS_MPI
430 <  strcpy(checkPointMsg, "Component stamps successfully extracted\n");
431 <  MPIcheckPoint();
430 >
431 >    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
432 >    MPIcheckPoint();
433 >
434   #endif // is_mpi
435 +
436   }
437  
438 < void SimCreator::setGlobalIndices(SimInfo* info) {
438 > void SimCreator::setGlobalIndex(SimInfo *info) {
439      typename SimInfo::MoleculeIterator mi;
440      typename Molecule::AtomIterator ai;
441      typename Molecule::RigidBodyIterator ri;
442      typename Molecule::CutoffGroupIterator ci;
443 <    Molecule* mol;
444 <    Atom* atom;
445 <    RigidBody* rb;
446 <    CutoffGroup* cg;
443 >    Molecule * mol;
444 >    Atom * atom;
445 >    RigidBody * rb;
446 >    CutoffGroup * cg;
447      int beginAtomIndex;
448      int beginRigidBodyIndex;
449      int beginCutoffGroupIndex;
450  
451 < #ifdef IS_MPI
451 > #ifndef IS_MPI
452 >
453      beginAtomIndex = 0;
454      beginRigidBodyIndex = 0;
455      beginCutoffGroupIndex = 0;
456 +
457   #else
458 +
459      int nproc;
460      int myNode;
461  
462      myNode = worldRank;
463      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
479 –    
480 –    std::vector<int> tmpAtomsInProc(nproc, 0);
481 –    std::vector<int> tmpRigidBodiesInProc(nproc, 0);
482 –    std::vector<int> tmpCutoffGroupsInProc(nproc, 0);    
483 –    std::vector<int> NumAtomsInProc(nproc, 0);
484 –    std::vector<int> NumRigidBodiesInProc(nproc, 0);
485 –    std::vector<int> NumCutoffGroupsInProc(nproc, 0);
464  
465 <    tmpAtomInProc[myNode] = info->getNAtoms();
466 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodiess();
467 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroupss();
468 <    
465 >    std::vector < int > tmpAtomsInProc(nproc, 0);
466 >    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
467 >    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
468 >    std::vector < int > NumAtomsInProc(nproc, 0);
469 >    std::vector < int > NumRigidBodiesInProc(nproc, 0);
470 >    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
471 >
472 >    tmpAtomsInProc[myNode] = info->getNAtoms();
473 >    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
474 >    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
475 >
476      //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
477 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
478 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
479 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
477 >    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
478 >                  MPI_SUM, MPI_COMM_WORLD);
479 >    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
480 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
481 >    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
482 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
483  
484      beginAtomIndex = 0;
485      beginRigidBodyIndex = 0;
486      beginCutoffGroupIndex = 0;
487  
488 <    for (int i = 0; i < nproc; i++) {
488 >    for(int i = 0; i < nproc; i++) {
489          beginAtomIndex += NumAtomsInProc[i];
490          beginRigidBodyIndex += NumRigidBodiesInProc[i];
491 <        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];        
491 >        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
492      }
493 <    
493 >
494   #endif
507 –    
508 –    for (mol = info->beginMolecule(mi); mol != NULL; mol  = info->nextMolecule(mi)) {
495  
496 +    for(mol = info->beginMolecule(mi); mol != NULL;
497 +        mol = info->nextMolecule(mi)) {
498 +
499          //local index(index in DataStorge) of atom is important
500 <        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
500 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
501              atom->setGlobalIndex(beginAtomIndex++);
502          }
503  
504 <        for (rb = mol->beginRigidBody(ri); rb != NULL; rb = mol->nextRigidBody(ri)) {
504 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
505 >            rb = mol->nextRigidBody(ri)) {
506              rb->setGlobalIndex(beginRigidBodyIndex++);
507          }
508 <        
508 >
509          //local index of cutoff group is trivial, it only depends on the order of travesing
510 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
510 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
511 >            cg = mol->nextCutoffGroup(ci)) {
512              cg->setGlobalIndex(beginCutoffGroupIndex++);
513 <        }        
513 >        }
514 >    }
515 >
516 >    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
517 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
518 >        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
519 >
520 >            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
521 >                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
522 >            }
523 >
524 >        }      
525 >    }
526 >    
527 >    // Since the globalGroupMembership has been zero filled and we've only
528 >    // poked values into the atoms we know, we can do an Allreduce
529 >    // to get the full globalGroupMembership array (We think).
530 >    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
531 >    // docs said we could.
532 >
533 >    MPI_Allreduce(&globalGroupMembership[0],
534 >                  info->getGlobalGroupMembershipPointer(),
535 >                  info->getNGlobalAtoms(),
536 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
537 >
538 >    //fill molMembership
539 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
540 >    
541 >    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
542 >
543 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
544 >            globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
545 >        }
546 >
547 >    MPI_Allreduce(&globalMolMembership[0],
548 >                  info->getGlobalMolMembershipPointer(),
549 >                  info->getNGlobalAtoms(),
550 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
551 >        
552 > }
553 >
554 >
555 > void SimCreator::initFortran(SimInfo* info) {
556 >    info->update();
557 >
558 >    std::vector<AtomType*> atomTypes;
559 >    std::vector<AtomType*>::iterator i;
560 >    std::vector<double> cutoffRadius;
561 >    ForceField* ff;
562 >
563 >    //get the unique atom types
564 >    atomTypes = info->getUniqueAtomTypes();
565 >
566 >    ff = info->getForceField();
567 >
568 >    //query the max cutoff radius among these atom types
569 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
570 >        cutoffRadius.push_back(ff->(*i));
571 >    }
572 >
573 >    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
574 >
575 > #ifdef IS_MPI
576 >    //pick the max cutoff radius among the processors
577 > #endif
578 >
579 >    
580 >    
581 > }
582 >
583 > void SimCreator::loadCoordinates(SimInfo* info) {
584 >    Globals* globals;
585 >    globals = info->getGlobals();
586 >    
587 >    if (!globals->haveInitialConfig()) {
588 >        sprintf(painCave.errMsg,
589 >                "Cannot intialize a simulation without an initial configuration file.\n");
590 >        painCave.isFatal = 1;;
591 >        simError();
592 >    }
593          
594 +    DumpReader reader(globals->getInitialConfig());
595 +    int nframes = reader->getNframes();
596 +
597 +    if (nframes > 0) {
598 +        reader.readFrame(info, nframes - 1);
599 +    } else {
600 +    //invalid initial coordinate file
601 +        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
602 +                globals->getInitialConfig());
603 +        painCave.isFatal = 1;
604 +        simError();
605      }
606 +
607   }
608  
609   } //end namespace oopse

Diff Legend

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