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

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

Diff Legend

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