ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/SimSetup.cpp
(Generate patch)

Comparing trunk/mdtools/interface_implementation/SimSetup.cpp (file contents):
Revision 11 by mmeineke, Tue Jul 9 18:40:59 2002 UTC vs.
Revision 206 by chuckv, Thu Dec 12 21:21:59 2002 UTC

# Line 6 | Line 6
6   #include "parse_me.h"
7   #include "LRI.hpp"
8   #include "Integrator.hpp"
9 + #include "simError.h"
10  
11 + #ifdef IS_MPI
12 + #include "mpiBASS.h"
13 + #include "bassDiag.hpp"
14 + #endif
15 +
16   SimSetup::SimSetup(){
17    stamps = new MakeStamps();
18    globals = new Globals();
19 +  
20 + #ifdef IS_MPI
21 +  strcpy( checkPointMsg, "SimSetup creation successful" );
22 +  MPIcheckPoint();
23 + #endif // IS_MPI
24   }
25  
26   SimSetup::~SimSetup(){
# Line 19 | Line 30 | void SimSetup::parseFile( char* fileName ){
30  
31   void SimSetup::parseFile( char* fileName ){
32  
33 <  inFileName = fileName;
34 <  set_interface_stamps( stamps, globals );
35 <  yacc_BASS( fileName );
33 > #ifdef IS_MPI
34 >  if( worldRank == 0 ){
35 > #endif // is_mpi
36 >    
37 >    inFileName = fileName;
38 >    set_interface_stamps( stamps, globals );
39 >    
40 > #ifdef IS_MPI
41 >    mpiEventInit();
42 > #endif
43 >
44 >    yacc_BASS( fileName );
45 >
46 > #ifdef IS_MPI
47 >    throwMPIEvent(NULL);
48 >  }
49 >  else receiveParse();
50 > #endif
51 >
52   }
53  
54 + #ifdef IS_MPI
55 + void SimSetup::receiveParse(void){
56 +
57 +    set_interface_stamps( stamps, globals );
58 +    mpiEventInit();
59 +    MPIcheckPoint();
60 +    mpiEventLoop();
61 +
62 + }
63 +
64 +
65 + void SimSetup::testMe(void){
66 +  bassDiag* dumpMe = new bassDiag(globals,stamps);
67 +  dumpMe->dumpStamps();
68 +  delete dumpMe;
69 + }
70 + #endif
71 +
72   void SimSetup::createSim( void ){
73  
74    MakeStamps *the_stamps;
75    Globals* the_globals;
76 <  int i;
76 >  int i, j;
77  
78    // get the stamps and globals;
79    the_stamps = stamps;
# Line 43 | Line 88 | void SimSetup::createSim( void ){
88    n_components = the_globals->getNComponents();
89    strcpy( force_field, the_globals->getForceField() );
90    strcpy( ensemble, the_globals->getEnsemble() );
91 <  
91 >
92    if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
93    else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
94    else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
95    else{
96 <    std::cerr<< "SimSetup Error. Unrecognized force field -> "
97 <             << force_field << "\n";
98 <    exit(8);
96 >    sprintf( painCave.errMsg,
97 >             "SimSetup Error. Unrecognized force field -> %s\n",
98 >             force_field );
99 >    painCave.isFatal = 1;
100 >    simError();
101    }
102  
103 + #ifdef IS_MPI
104 +  strcpy( checkPointMsg, "ForceField creation successful" );
105 +  MPIcheckPoint();
106 + #endif // is_mpi
107 +
108    // get the components and calculate the tot_nMol and indvidual n_mol
109    the_components = the_globals->getComponents();
110    components_nmol = new int[n_components];
111    comp_stamps = new MoleculeStamp*[n_components];
112 <  
112 >
113    if( !the_globals->haveNMol() ){
114 <    // we don't have the total number of molecules, so we assume it is
114 >    // we don't have the total number of molecules, so we assume it is
115      // given in each component
116  
117      tot_nmol = 0;
118      for( i=0; i<n_components; i++ ){
119 <      
119 >
120        if( !the_components[i]->haveNMol() ){
121          // we have a problem
122 <        std::cerr << "SimSetup Error. No global NMol or component NMol"
123 <                  << " given. Cannot calculate the number of atoms.\n";
124 <        exit( 8 );
122 >        sprintf( painCave.errMsg,
123 >                 "SimSetup Error. No global NMol or component NMol"
124 >                 " given. Cannot calculate the number of atoms.\n" );
125 >        painCave.isFatal = 1;
126 >        simError();
127        }
128  
129        tot_nmol += the_components[i]->getNMol();
# Line 77 | Line 131 | void SimSetup::createSim( void ){
131      }
132    }
133    else{
134 <    std::cerr << "NOT A SUPPORTED FEATURE\n";
134 >    sprintf( painCave.errMsg,
135 >             "SimSetup error.\n"
136 >             "\tSorry, the ability to specify total"
137 >             " nMols and then give molfractions in the components\n"
138 >             "\tis not currently supported."
139 >             " Please give nMol in the components.\n" );
140 >    painCave.isFatal = 1;
141 >    simError();
142      
82 //     tot_nmol = the_globals->getNMol();
143      
144 < //     //we have the total number of molecules, now we check for molfractions
145 < //     for( i=0; i<n_components; i++ ){
146 <      
147 < //       if( !the_components[i]->haveMolFraction() ){
148 <        
149 < //      if( !the_components[i]->haveNMol() ){
150 < //        //we have a problem
151 < //        std::cerr << "SimSetup error. Neither molFraction nor "
152 < //                  << " nMol was given in component
153 <
144 >    //     tot_nmol = the_globals->getNMol();
145 >    
146 >    //   //we have the total number of molecules, now we check for molfractions
147 >    //     for( i=0; i<n_components; i++ ){
148 >    
149 >    //       if( !the_components[i]->haveMolFraction() ){
150 >    
151 >    //  if( !the_components[i]->haveNMol() ){
152 >    //    //we have a problem
153 >    //    std::cerr << "SimSetup error. Neither molFraction nor "
154 >    //              << " nMol was given in component
155 >    
156    }
157  
158 + #ifdef IS_MPI
159 +  strcpy( checkPointMsg, "Have the number of components" );
160 +  MPIcheckPoint();
161 + #endif // is_mpi
162 +
163    // make an array of molecule stamps that match the components used.
164 +  // also extract the used stamps out into a separate linked list
165  
166 +  simnfo->nComponents = n_components;
167 +  simnfo->componentsNmol = components_nmol;
168 +  simnfo->compStamps = comp_stamps;
169 +  simnfo->headStamp = new LinkedMolStamp();
170 +  
171 +  char* id;
172 +  LinkedMolStamp* headStamp = simnfo->headStamp;
173 +  LinkedMolStamp* currentStamp = NULL;
174    for( i=0; i<n_components; i++ ){
175  
176 <    comp_stamps[i] =
177 <      the_stamps->getMolecule( the_components[i]->getType() );
176 >    id = the_components[i]->getType();
177 >    comp_stamps[i] = NULL;
178 >    
179 >    // check to make sure the component isn't already in the list
180 >
181 >    comp_stamps[i] = headStamp->match( id );
182 >    if( comp_stamps[i] == NULL ){
183 >      
184 >      // extract the component from the list;
185 >      
186 >      currentStamp = the_stamps->extractMolStamp( id );
187 >      if( currentStamp == NULL ){
188 >        sprintf( painCave.errMsg,
189 >                 "SimSetup error: Component \"%s\" was not found in the "
190 >                 "list of declared molecules\n"
191 >                 id );
192 >        painCave.isFatal = 1;
193 >        simError();
194 >      }
195 >      
196 >      headStamp->add( currentStamp );
197 >      comp_stamps[i] = headStamp->match( id );
198 >    }
199    }
200  
201 + #ifdef IS_MPI
202 +  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
203 +  MPIcheckPoint();
204 + #endif // is_mpi
205    
206  
207 +
208 +
209    // caclulate the number of atoms, bonds, bends and torsions
210  
211    tot_atoms = 0;
# Line 111 | Line 214 | void SimSetup::createSim( void ){
214    tot_torsions = 0;
215    for( i=0; i<n_components; i++ ){
216      
217 <    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
218 <    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
219 <    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
217 >    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
218 >    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
219 >    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
220      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
221    }
222 <  
222 >
223    tot_SRI = tot_bonds + tot_bends + tot_torsions;
224 <  
224 >
225    simnfo->n_atoms = tot_atoms;
226    simnfo->n_bonds = tot_bonds;
227    simnfo->n_bends = tot_bends;
228    simnfo->n_torsions = tot_torsions;
229    simnfo->n_SRI = tot_SRI;
230 +  simnfo->n_mol = tot_nmol;
231  
128  // create the atom and short range interaction arrays
232    
233 <  the_atoms = new Atom*[tot_atoms];
234 <  //  the_molecules = new Molecule[tot_nmol];
233 > #ifdef IS_MPI
234 >
235 >  // divide the molecules among processors here.
236    
237 +  mpiSimulation* mpiSim = new mpiSimulation( simnfo );
238    
239 <  if( tot_SRI ){
240 <    the_sris = new SRI*[tot_SRI];
241 <    the_excludes = new ex_pair[tot_SRI];
239 >  mpiSim->divideLabor();
240 >
241 >  // set up the local variables
242 >  
243 >  int localMol, allMol;
244 >  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
245 >  
246 >  allMol = 0;
247 >  localMol = 0;
248 >  local_atoms = 0;
249 >  local_bonds = 0;
250 >  local_bends = 0;
251 >  local_torsions = 0;
252 >  for( i=0; i<n_components; i++ ){
253 >
254 >    for( j=0; j<components_nmol[i]; j++ ){
255 >      
256 >      if( mpiSim->getMyMolStart() <= allMol &&
257 >          allMol <= mpiSim->getMyMolEnd() ){
258 >        
259 >        local_atoms +=    comp_stamps[i]->getNAtoms();
260 >        local_bonds +=    comp_stamps[i]->getNBonds();
261 >        local_bends +=    comp_stamps[i]->getNBends();
262 >        local_torsions += comp_stamps[i]->getNTorsions();
263 >        localMol++;
264 >      }      
265 >      allMol++;
266 >    }
267    }
268 +  local_SRI = local_bonds + local_bends + local_torsions;
269 +  
270  
271 +  simnfo->n_atoms = mpiSim->getMyNlocal();  
272 +  
273 +  if( local_atoms != simnfo->n_atoms ){
274 +    sprintf( painCave.errMsg,
275 +             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
276 +             " localAtom (%d) are note equal.\n",
277 +             simnfo->n_atoms,
278 +             local_atoms );
279 +    painCave.isFatal = 1;
280 +    simError();
281 +  }
282 +
283 +  simnfo->n_bonds = local_bonds;
284 +  simnfo->n_bends = local_bends;
285 +  simnfo->n_torsions = local_torsions;
286 +  simnfo->n_SRI = local_SRI;
287 +  simnfo->n_mol = localMol;
288 +
289 +  strcpy( checkPointMsg, "Passed nlocal consistency check." );
290 +  MPIcheckPoint();
291 +  
292 +  
293 + #endif // is_mpi
294 +  
295 +
296 +  // create the atom and short range interaction arrays
297 +
298 +  Atom::createArrays(simnfo->n_atoms);
299 +  the_atoms = new Atom*[simnfo->n_atoms];
300 +  the_molecules = new Molecule[simnfo->n_mol];
301 +
302 +
303 +  if( simnfo->n_SRI ){
304 +    the_sris = new SRI*[simnfo->n_SRI];
305 +    the_excludes = new ex_pair[simnfo->n_SRI];
306 +  }
307 +
308    // set the arrays into the SimInfo object
309  
310    simnfo->atoms = the_atoms;
# Line 143 | Line 312 | void SimSetup::createSim( void ){
312    simnfo->n_exclude = tot_SRI;
313    simnfo->excludes = the_excludes;
314  
315 +
316    // initialize the arrays
317 <  
317 >
318    the_ff->setSimInfo( simnfo );
319 <    
319 >
320    makeAtoms();
321  
322    if( tot_bonds ){
# Line 161 | Line 331 | void SimSetup::createSim( void ){
331      makeTorsions();
332    }
333  
164  //  makeMolecules();
334  
335    // get some of the tricky things that may still be in the globals
336  
337    if( simnfo->n_dipoles ){
338  
339      if( !the_globals->haveRRF() ){
340 <      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
341 <      exit(8);
340 >      sprintf( painCave.errMsg,
341 >               "SimSetup Error, system has dipoles, but no rRF was set.\n");
342 >      painCave.isFatal = 1;
343 >      simError();
344      }
345      if( !the_globals->haveDielectric() ){
346 <      std::cerr << "SimSetup Error, system has dipoles, but no"
347 <                << " dielectric was set.\n";
348 <      exit(8);
346 >      sprintf( painCave.errMsg,
347 >               "SimSetup Error, system has dipoles, but no"
348 >               " dielectric was set.\n" );
349 >      painCave.isFatal = 1;
350 >      simError();
351      }
352  
353      simnfo->rRF        = the_globals->getRRF();
354      simnfo->dielectric = the_globals->getDielectric();
355    }
356  
357 + #ifdef IS_MPI
358 +  strcpy( checkPointMsg, "rRf and dielectric check out" );
359 +  MPIcheckPoint();
360 + #endif // is_mpi
361 +  
362    if( the_globals->haveBox() ){
363      simnfo->box_x = the_globals->getBox();
364      simnfo->box_y = the_globals->getBox();
365      simnfo->box_z = the_globals->getBox();
366    }
367    else if( the_globals->haveDensity() ){
368 <    
368 >
369      double vol;
370      vol = (double)tot_nmol / the_globals->getDensity();
371      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 196 | Line 374 | void SimSetup::createSim( void ){
374    }
375    else{
376      if( !the_globals->haveBoxX() ){
377 <      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
378 <      exit(8);
377 >      sprintf( painCave.errMsg,
378 >               "SimSetup error, no periodic BoxX size given.\n" );
379 >      painCave.isFatal = 1;
380 >      simError();
381      }
382      simnfo->box_x = the_globals->getBoxX();
383  
384      if( !the_globals->haveBoxY() ){
385 <      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
386 <      exit(8);
385 >      sprintf( painCave.errMsg,
386 >               "SimSetup error, no periodic BoxY size given.\n" );
387 >      painCave.isFatal = 1;
388 >      simError();
389      }
390      simnfo->box_y = the_globals->getBoxY();
391  
392      if( !the_globals->haveBoxZ() ){
393 <      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
394 <      exit(8);
393 >      sprintf( painCave.errMsg,
394 >               "SimSetup error, no periodic BoxZ size given.\n" );
395 >      painCave.isFatal = 1;
396 >      simError();
397      }
398      simnfo->box_z = the_globals->getBoxZ();
399    }
216    
217  if( the_globals->haveInitialConfig() ){
218    InitializeFromFile* fileInit;
219    fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
220    
221    fileInit->read_xyz( simnfo ); // default velocities on
400  
401 <    delete fileInit;    
402 <  }
403 <  else{
404 <    initFromBass();
227 <  }
401 > #ifdef IS_MPI
402 >  strcpy( checkPointMsg, "Box size set up" );
403 >  MPIcheckPoint();
404 > #endif // is_mpi
405  
406 <  if( the_globals->haveFinalConfig() ){
407 <    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
408 <  }
409 <  else{
410 <    strcpy( simnfo->finalName, inFileName );
411 <    char* endTest;
412 <    int nameLength = strlen( simnfo->finalName );
413 <    endTest = &(simnfo->finalName[nameLength - 5]);
414 <    if( !strcmp( endTest, ".bass" ) ){
415 <      strcpy( endTest, ".eor" );
406 >
407 >
408 > if( the_globals->haveInitialConfig() ){
409 >
410 >     InitializeFromFile* fileInit;
411 > #ifdef IS_MPI // is_mpi
412 >     if( worldRank == 0 ){
413 > #endif //is_mpi
414 >   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
415 > #ifdef IS_MPI
416 >     }else fileInit = new InitializeFromFile( NULL );
417 > #endif
418 >   fileInit->read_xyz( simnfo ); // default velocities on
419 >
420 >   delete fileInit;
421 > }
422 > else{
423 >
424 > #ifdef IS_MPI
425 >
426 >  // no init from bass
427 >  
428 >  sprintf( painCave.errMsg,
429 >           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
430 >  painCave.isFatal;
431 >  simError();
432 >  
433 > #else
434 >
435 >  initFromBass();
436 >
437 >
438 > #endif
439 > }
440 >
441 > #ifdef IS_MPI
442 >  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
443 >  MPIcheckPoint();
444 > #endif // is_mpi
445 >
446 >
447 >  
448 >
449 >  
450 >
451 >  
452 > #ifdef IS_MPI
453 >  if( worldRank == 0 ){
454 > #endif // is_mpi
455 >    
456 >    if( the_globals->haveFinalConfig() ){
457 >      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
458      }
240    else if( !strcmp( endTest, ".BASS" ) ){
241      strcpy( endTest, ".eor" );
242    }
459      else{
460 <      endTest = &(simnfo->finalName[nameLength - 4]);
461 <      if( !strcmp( endTest, ".bss" ) ){
460 >      strcpy( simnfo->finalName, inFileName );
461 >      char* endTest;
462 >      int nameLength = strlen( simnfo->finalName );
463 >      endTest = &(simnfo->finalName[nameLength - 5]);
464 >      if( !strcmp( endTest, ".bass" ) ){
465          strcpy( endTest, ".eor" );
466        }
467 <      else if( !strcmp( endTest, ".mdl" ) ){
467 >      else if( !strcmp( endTest, ".BASS" ) ){
468          strcpy( endTest, ".eor" );
469        }
470        else{
471 <        strcat( simnfo->finalName, ".eor" );
471 >        endTest = &(simnfo->finalName[nameLength - 4]);
472 >        if( !strcmp( endTest, ".bss" ) ){
473 >          strcpy( endTest, ".eor" );
474 >        }
475 >        else if( !strcmp( endTest, ".mdl" ) ){
476 >          strcpy( endTest, ".eor" );
477 >        }
478 >        else{
479 >          strcat( simnfo->finalName, ".eor" );
480 >        }
481        }
482      }
255  }
483      
484 <  // make the sample and status out names
485 <
486 <  strcpy( simnfo->sampleName, inFileName );
487 <  char* endTest;
488 <  int nameLength = strlen( simnfo->sampleName );
489 <  endTest = &(simnfo->sampleName[nameLength - 5]);
490 <  if( !strcmp( endTest, ".bass" ) ){
264 <    strcpy( endTest, ".dump" );
265 <  }
266 <  else if( !strcmp( endTest, ".BASS" ) ){
267 <    strcpy( endTest, ".dump" );
268 <  }
269 <  else{
270 <    endTest = &(simnfo->sampleName[nameLength - 4]);
271 <    if( !strcmp( endTest, ".bss" ) ){
484 >    // make the sample and status out names
485 >    
486 >    strcpy( simnfo->sampleName, inFileName );
487 >    char* endTest;
488 >    int nameLength = strlen( simnfo->sampleName );
489 >    endTest = &(simnfo->sampleName[nameLength - 5]);
490 >    if( !strcmp( endTest, ".bass" ) ){
491        strcpy( endTest, ".dump" );
492      }
493 <    else if( !strcmp( endTest, ".mdl" ) ){
493 >    else if( !strcmp( endTest, ".BASS" ) ){
494        strcpy( endTest, ".dump" );
495      }
496      else{
497 <      strcat( simnfo->sampleName, ".dump" );
497 >      endTest = &(simnfo->sampleName[nameLength - 4]);
498 >      if( !strcmp( endTest, ".bss" ) ){
499 >        strcpy( endTest, ".dump" );
500 >      }
501 >      else if( !strcmp( endTest, ".mdl" ) ){
502 >        strcpy( endTest, ".dump" );
503 >      }
504 >      else{
505 >        strcat( simnfo->sampleName, ".dump" );
506 >      }
507      }
508 <  }
509 <  
510 <  strcpy( simnfo->statusName, inFileName );
511 <  nameLength = strlen( simnfo->statusName );
512 <  endTest = &(simnfo->statusName[nameLength - 5]);
285 <  if( !strcmp( endTest, ".bass" ) ){
286 <    strcpy( endTest, ".stat" );
287 <  }
288 <  else if( !strcmp( endTest, ".BASS" ) ){
289 <    strcpy( endTest, ".stat" );
290 <  }
291 <  else{
292 <    endTest = &(simnfo->statusName[nameLength - 4]);
293 <    if( !strcmp( endTest, ".bss" ) ){
508 >    
509 >    strcpy( simnfo->statusName, inFileName );
510 >    nameLength = strlen( simnfo->statusName );
511 >    endTest = &(simnfo->statusName[nameLength - 5]);
512 >    if( !strcmp( endTest, ".bass" ) ){
513        strcpy( endTest, ".stat" );
514      }
515 <    else if( !strcmp( endTest, ".mdl" ) ){
515 >    else if( !strcmp( endTest, ".BASS" ) ){
516        strcpy( endTest, ".stat" );
517      }
518      else{
519 <      strcat( simnfo->statusName, ".stat" );
519 >      endTest = &(simnfo->statusName[nameLength - 4]);
520 >      if( !strcmp( endTest, ".bss" ) ){
521 >        strcpy( endTest, ".stat" );
522 >      }
523 >      else if( !strcmp( endTest, ".mdl" ) ){
524 >        strcpy( endTest, ".stat" );
525 >      }
526 >      else{
527 >        strcat( simnfo->statusName, ".stat" );
528 >      }
529      }
530 +    
531 + #ifdef IS_MPI
532    }
533 + #endif // is_mpi
534    
535    // set the status, sample, and themal kick times
536 <
536 >  
537    if( the_globals->haveSampleTime() ){
538 <    simnfo->sampleTime = the_globals->getSampleTime();
538 >    simnfo->sampleTime = the_globals->getSampleTime();
539      simnfo->statusTime = simnfo->sampleTime;
540      simnfo->thermalTime = simnfo->sampleTime;
541    }
542    else{
543 <    simnfo->sampleTime = the_globals->getRunTime();
543 >    simnfo->sampleTime = the_globals->getRunTime();
544      simnfo->statusTime = simnfo->sampleTime;
545      simnfo->thermalTime = simnfo->sampleTime;
546    }
# Line 325 | Line 556 | void SimSetup::createSim( void ){
556    // check for the temperature set flag
557  
558    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
559 <  
560 <    
559 >
560 >
561    // make the longe range forces and the integrator
562 <  
562 >
563    new AllLong( simnfo );
564 <      
564 >
565    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
566    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
567    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
568   }
569  
570   void SimSetup::makeAtoms( void ){
571 <  
571 >
572    int i, j, k, index;
573    double ux, uy, uz, uSqr, u;
574    AtomStamp* current_atom;
575    DirectionalAtom* dAtom;
576 +  int molIndex, molStart, molEnd, nMemb, lMolIndex;
577  
578 +  lMolIndex = 0;
579 +  molIndex = 0;
580    index = 0;
581    for( i=0; i<n_components; i++ ){
582 <    
582 >
583      for( j=0; j<components_nmol[i]; j++ ){
350      
351      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
352        
353        current_atom = comp_stamps[i]->getAtom( k );
354        if( current_atom->haveOrientation() ){
584  
585 <          dAtom = new DirectionalAtom;
586 <          simnfo->n_oriented++;
587 <          the_atoms[index] = dAtom;
588 <      
589 <          ux = current_atom->getOrntX();
590 <          uy = current_atom->getOrntY();
591 <          uz = current_atom->getOrntZ();
585 > #ifdef IS_MPI
586 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
587 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
588 > #endif // is_mpi        
589 >
590 >        molStart = index;
591 >        nMemb = comp_stamps[i]->getNAtoms();
592 >        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
593            
594 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
594 >          current_atom = comp_stamps[i]->getAtom( k );
595 >          if( current_atom->haveOrientation() ){
596 >            
597 >            dAtom = new DirectionalAtom(index);
598 >            simnfo->n_oriented++;
599 >            the_atoms[index] = dAtom;
600 >            
601 >            ux = current_atom->getOrntX();
602 >            uy = current_atom->getOrntY();
603 >            uz = current_atom->getOrntZ();
604 >            
605 >            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
606 >            
607 >            u = sqrt( uSqr );
608 >            ux = ux / u;
609 >            uy = uy / u;
610 >            uz = uz / u;
611 >            
612 >            dAtom->setSUx( ux );
613 >            dAtom->setSUy( uy );
614 >            dAtom->setSUz( uz );
615 >          }
616 >          else{
617 >            the_atoms[index] = new GeneralAtom(index);
618 >          }
619 >          the_atoms[index]->setType( current_atom->getType() );
620 >          the_atoms[index]->setIndex( index );
621            
622 <          u = sqrt( uSqr );
623 <          ux = ux / u;
368 <          uy = uy / u;
369 <          uz = uz / u;
370 <          
371 <          dAtom->setSUx( ux );
372 <          dAtom->setSUy( uy );
373 <          dAtom->setSUz( uz );
622 >          // increment the index and repeat;
623 >          index++;
624          }
375        else{
376          the_atoms[index] = new GeneralAtom;
377        }
378        the_atoms[index]->setType( current_atom->getType() );
379        the_atoms[index]->setIndex( index );
625          
626 <        // increment the index and repeat;
627 <        index++;
626 >        molEnd = index -1;
627 >        the_molecules[lMolIndex].setNMembers( nMemb );
628 >        the_molecules[lMolIndex].setStartAtom( molStart );
629 >        the_molecules[lMolIndex].setEndAtom( molEnd );
630 >        the_molecules[lMolIndex].setStampID( i );
631 >        lMolIndex++;
632 >
633 > #ifdef IS_MPI
634        }
635 + #endif //is_mpi
636 +      
637 +      molIndex++;
638      }
639    }
640 <  
640 >
641    the_ff->initializeAtoms();
642   }
643  
644   void SimSetup::makeBonds( void ){
645  
646 <  int i, j, k, index, offset;
646 >  int i, j, k, index, offset, molIndex;
647    bond_pair* the_bonds;
648    BondStamp* current_bond;
649  
650    the_bonds = new bond_pair[tot_bonds];
651    index = 0;
652    offset = 0;
653 +  molIndex = 0;
654    for( i=0; i<n_components; i++ ){
655 <    
655 >
656      for( j=0; j<components_nmol[i]; j++ ){
402      
403      for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
404        
405        current_bond = comp_stamps[i]->getBond( k );
406        the_bonds[index].a = current_bond->getA() + offset;
407        the_bonds[index].b = current_bond->getB() + offset;
657  
658 <        the_excludes[index].i = the_bonds[index].a;
659 <        the_excludes[index].j = the_bonds[index].b;
660 <
661 <        // increment the index and repeat;
662 <        index++;
663 <      }
664 <      offset += comp_stamps[i]->getNAtoms();
665 <    }
658 > #ifdef IS_MPI
659 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
660 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
661 > #endif // is_mpi        
662 >        
663 >        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
664 >          
665 >          current_bond = comp_stamps[i]->getBond( k );
666 >          the_bonds[index].a = current_bond->getA() + offset;
667 >          the_bonds[index].b = current_bond->getB() + offset;
668 >          
669 >          the_excludes[index].i = the_bonds[index].a;
670 >          the_excludes[index].j = the_bonds[index].b;
671 >          
672 >          // increment the index and repeat;
673 >          index++;
674 >        }
675 >        offset += comp_stamps[i]->getNAtoms();
676 >        
677 > #ifdef IS_MPI
678 >      }
679 > #endif is_mpi
680 >      
681 >      molIndex++;
682 >    }      
683    }
684 <  
684 >
685    the_ff->initializeBonds( the_bonds );
686   }
687  
688   void SimSetup::makeBends( void ){
689  
690 <  int i, j, k, index, offset;
690 >  int i, j, k, index, offset, molIndex;
691    bend_set* the_bends;
692    BendStamp* current_bend;
693  
694    the_bends = new bend_set[tot_bends];
695    index = 0;
696    offset = 0;
697 +  molIndex = 0;
698    for( i=0; i<n_components; i++ ){
699 <    
699 >
700      for( j=0; j<components_nmol[i]; j++ ){
434      
435      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
436        
437        current_bend = comp_stamps[i]->getBend( k );
438        the_bends[index].a = current_bend->getA() + offset;
439        the_bends[index].b = current_bend->getB() + offset;
440        the_bends[index].c = current_bend->getC() + offset;
701  
702 <        the_excludes[index + tot_bonds].i = the_bends[index].a;
703 <        the_excludes[index + tot_bonds].j = the_bends[index].c;
702 > #ifdef IS_MPI
703 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
704 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
705 > #endif // is_mpi        
706  
707 <        // increment the index and repeat;
708 <        index++;
707 >        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
708 >          
709 >          current_bend = comp_stamps[i]->getBend( k );
710 >          the_bends[index].a = current_bend->getA() + offset;
711 >          the_bends[index].b = current_bend->getB() + offset;
712 >          the_bends[index].c = current_bend->getC() + offset;
713 >          
714 >          the_excludes[index + tot_bonds].i = the_bends[index].a;
715 >          the_excludes[index + tot_bonds].j = the_bends[index].c;
716 >          
717 >          // increment the index and repeat;
718 >          index++;
719 >        }
720 >        offset += comp_stamps[i]->getNAtoms();
721 >        
722 > #ifdef IS_MPI
723        }
724 <      offset += comp_stamps[i]->getNAtoms();
724 > #endif //is_mpi
725 >
726 >      molIndex++;
727      }
728    }
729 <  
729 >
730    the_ff->initializeBends( the_bends );
731   }
732  
733   void SimSetup::makeTorsions( void ){
734  
735 <  int i, j, k, index, offset;
735 >  int i, j, k, index, offset, molIndex;
736    torsion_set* the_torsions;
737    TorsionStamp* current_torsion;
738  
739    the_torsions = new torsion_set[tot_torsions];
740    index = 0;
741    offset = 0;
742 +  molIndex = 0;
743    for( i=0; i<n_components; i++ ){
744 <    
744 >
745      for( j=0; j<components_nmol[i]; j++ ){
746 <      
746 >
747 > #ifdef IS_MPI
748 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
749 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
750 > #endif // is_mpi        
751 >
752        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
753 <        
753 >
754          current_torsion = comp_stamps[i]->getTorsion( k );
755          the_torsions[index].a = current_torsion->getA() + offset;
756          the_torsions[index].b = current_torsion->getB() + offset;
# Line 480 | Line 764 | void SimSetup::makeTorsions( void ){
764          index++;
765        }
766        offset += comp_stamps[i]->getNAtoms();
767 +
768 + #ifdef IS_MPI
769 +      }
770 + #endif //is_mpi      
771 +
772 +      molIndex++;
773      }
774    }
775 <  
775 >
776    the_ff->initializeTorsions( the_torsions );
777   }
778  
489 void SimSetup::makeMolecules( void ){
490
491  //empy for now
492 }
493
779   void SimSetup::initFromBass( void ){
780  
781    int i, j, k;
# Line 518 | Line 803 | void SimSetup::initFromBass( void ){
803      n_per_extra = (int)ceil( temp1 );
804  
805      if( n_per_extra > 4){
806 <      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
807 <      exit(8);
806 >      sprintf( painCave.errMsg,
807 >               "SimSetup error. There has been an error in constructing"
808 >               " the non-complete lattice.\n" );
809 >      painCave.isFatal = 1;
810 >      simError();
811      }
812    }
813    else{
# Line 528 | Line 816 | void SimSetup::initFromBass( void ){
816      celly = simnfo->box_y / temp3;
817      cellz = simnfo->box_z / temp3;
818    }
819 <  
819 >
820    current_mol = 0;
821    current_comp_mol = 0;
822    current_comp = 0;
823    current_atom_ndx = 0;
824 <  
824 >
825    for( i=0; i < n_cells ; i++ ){
826      for( j=0; j < n_cells; j++ ){
827        for( k=0; k < n_cells; k++ ){
828 <        
828 >
829          makeElement( i * cellx,
830                       j * celly,
831                       k * cellz );
832 <        
832 >
833          makeElement( i * cellx + 0.5 * cellx,
834                       j * celly + 0.5 * celly,
835                       k * cellz );
836 <        
836 >
837          makeElement( i * cellx,
838                       j * celly + 0.5 * celly,
839                       k * cellz + 0.5 * cellz );
840 <        
840 >
841          makeElement( i * cellx + 0.5 * cellx,
842                       j * celly,
843                       k * cellz + 0.5 * cellz );
# Line 559 | Line 847 | void SimSetup::initFromBass( void ){
847  
848    if( have_extra ){
849      done = 0;
850 <    
850 >
851      int start_ndx;
852      for( i=0; i < (n_cells+1) && !done; i++ ){
853        for( j=0; j < (n_cells+1) && !done; j++ ){
854 <        
854 >
855          if( i < n_cells ){
856 <          
856 >
857            if( j < n_cells ){
858              start_ndx = n_cells;
859            }
860            else start_ndx = 0;
861          }
862          else start_ndx = 0;
863 <        
863 >
864          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
865 <          
865 >
866            makeElement( i * cellx,
867                         j * celly,
868                         k * cellz );
869            done = ( current_mol >= tot_nmol );
870 <          
870 >
871            if( !done && n_per_extra > 1 ){
872              makeElement( i * cellx + 0.5 * cellx,
873                           j * celly + 0.5 * celly,
874                           k * cellz );
875              done = ( current_mol >= tot_nmol );
876            }
877 <          
877 >
878            if( !done && n_per_extra > 2){
879              makeElement( i * cellx,
880                           j * celly + 0.5 * celly,
881                           k * cellz + 0.5 * cellz );
882              done = ( current_mol >= tot_nmol );
883            }
884 <          
884 >
885            if( !done && n_per_extra > 3){
886              makeElement( i * cellx + 0.5 * cellx,
887                           j * celly,
# Line 604 | Line 892 | void SimSetup::initFromBass( void ){
892        }
893      }
894    }
895 <  
896 <  
895 >
896 >
897    for( i=0; i<simnfo->n_atoms; i++ ){
898      simnfo->atoms[i]->set_vx( 0.0 );
899      simnfo->atoms[i]->set_vy( 0.0 );
# Line 621 | Line 909 | void SimSetup::makeElement( double x, double y, double
909    double rotMat[3][3];
910  
911    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
912 <    
912 >
913      current_atom = comp_stamps[current_comp]->getAtom( k );
914      if( !current_atom->havePosition() ){
915 <      std::cerr << "Component " << comp_stamps[current_comp]->getID()
916 <                << ", atom " << current_atom->getType()
917 <                << " does not have a position specified.\n"
918 <                << "The initialization routine is unable to give a start"
919 <                << " position.\n";
920 <      exit(8);
915 >      sprintf( painCave.errMsg,
916 >               "SimSetup:initFromBass error.\n"
917 >               "\tComponent %s, atom %s does not have a position specified.\n"
918 >               "\tThe initialization routine is unable to give a start"
919 >               " position.\n",
920 >               comp_stamps[current_comp]->getID(),
921 >               current_atom->getType() );
922 >      painCave.isFatal = 1;
923 >      simError();
924      }
925 <    
925 >
926      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
927      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
928      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
929 <    
929 >
930      if( the_atoms[current_atom_ndx]->isDirectional() ){
931 <      
931 >
932        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
933 <      
933 >
934        rotMat[0][0] = 1.0;
935        rotMat[0][1] = 0.0;
936        rotMat[0][2] = 0.0;
# Line 657 | Line 948 | void SimSetup::makeElement( double x, double y, double
948  
949      current_atom_ndx++;
950    }
951 <  
951 >
952    current_mol++;
953    current_comp_mol++;
954  
955    if( current_comp_mol >= components_nmol[current_comp] ){
956 <    
956 >
957      current_comp_mol = 0;
958      current_comp++;
959    }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines