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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC vs.
Revision 621 by gezelter, Wed Jul 16 02:11:02 2003 UTC

# Line 14 | Line 14
14  
15   // some defines for ensemble and Forcefield  cases
16  
17 < #define NVE_ENS  0
18 < #define NVT_ENS  1
19 < #define NPTi_ENS 2
20 < #define NPTf_ENS 3
17 > #define NVE_ENS   0
18 > #define NVT_ENS   1
19 > #define NPTi_ENS  2
20 > #define NPTf_ENS  3
21 > #define NPTim_ENS 4
22 > #define NPTfm_ENS 5
23  
24  
25   #define FF_DUFF 0
# Line 76 | Line 78 | void SimSetup::createSim( void ){
78  
79   void SimSetup::createSim( void ){
80  
79 –  MakeStamps *the_stamps;
80 –  Globals* the_globals;
81    int i, j, k, globalAtomIndex;
82    
83 <  int ensembleCase;
84 <  int ffCase;
83 >  // gather all of the information from the Bass file
84    
85 <  ensembleCase = -1;
87 <  ffCase = -1;
85 >  gatherInfo();
86  
87 <  // get the stamps and globals;
90 <  the_stamps = stamps;
91 <  the_globals = globals;
87 >  // creation of complex system objects
88  
89 <  // set the easy ones first
94 <  simnfo->target_temp = the_globals->getTargetTemp();
95 <  simnfo->dt = the_globals->getDt();
96 <  simnfo->run_time = the_globals->getRunTime();
89 >  sysObjectsCreation();
90  
91 <  // get the ones we know are there, yet still may need some work.
99 <  n_components = the_globals->getNComponents();
100 <  strcpy( force_field, the_globals->getForceField() );
101 <
102 <  if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
103 <  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
104 <  else{
105 <    sprintf( painCave.errMsg,
106 <             "SimSetup Error. Unrecognized force field -> %s\n",
107 <             force_field );
108 <    painCave.isFatal = 1;
109 <    simError();
110 <  }
111 <
112 <  // get the ensemble:
113 <  strcpy( ensemble, the_globals->getEnsemble() );
114 <
115 <  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
116 <  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
117 <  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") )
118 <    ensembleCase = NPTi_ENS;
119 <  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
120 <  else{
121 <    sprintf( painCave.errMsg,
122 <             "SimSetup Warning. Unrecognized Ensemble -> %s, "
123 <             "reverting to NVE for this simulation.\n",
124 <             ensemble );
125 <    painCave.isFatal = 0;
126 <    simError();
127 <    strcpy( ensemble, "NVE" );
128 <    ensembleCase = NVE_ENS;
129 <  }  
130 <  strcpy( simnfo->ensemble, ensemble );
131 <
132 <
133 < //   if( !strcasecmp( ensemble, "NPT" ) ) {
134 < //     the_extendedsystem = new ExtendedSystem( simnfo );
135 < //     the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
136 < //     if (the_globals->haveTargetPressure())
137 < //       the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
138 < //     else {
139 < //       sprintf( painCave.errMsg,
140 < //                "SimSetup error: If you use the constant pressure\n"
141 < //                "    ensemble, you must set targetPressure.\n"
142 < //                "    This was found in the BASS file.\n");
143 < //       painCave.isFatal = 1;
144 < //       simError();
145 < //     }
146 <
147 < //     if (the_globals->haveTauThermostat())
148 < //       the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
149 < //     else if (the_globals->haveQmass())
150 < //       the_extendedsystem->setQmass(the_globals->getQmass());
151 < //     else {
152 < //       sprintf( painCave.errMsg,
153 < //                "SimSetup error: If you use one of the constant temperature\n"
154 < //                "    ensembles, you must set either tauThermostat or qMass.\n"
155 < //                "    Neither of these was found in the BASS file.\n");
156 < //       painCave.isFatal = 1;
157 < //       simError();
158 < //     }
159 <
160 < //     if (the_globals->haveTauBarostat())
161 < //       the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
162 < //     else {
163 < //       sprintf( painCave.errMsg,
164 < //                "SimSetup error: If you use the constant pressure\n"
165 < //                "    ensemble, you must set tauBarostat.\n"
166 < //                "    This was found in the BASS file.\n");
167 < //       painCave.isFatal = 1;
168 < //       simError();
169 < //     }
170 <
171 < //   } else if ( !strcasecmp( ensemble, "NVT") ) {
172 < //     the_extendedsystem = new ExtendedSystem( simnfo );
173 < //     the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
174 <
175 < //     if (the_globals->haveTauThermostat())
176 < //       the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
177 < //     else if (the_globals->haveQmass())
178 < //       the_extendedsystem->setQmass(the_globals->getQmass());
179 < //     else {
180 < //       sprintf( painCave.errMsg,
181 < //                "SimSetup error: If you use one of the constant temperature\n"
182 < //                "    ensembles, you must set either tauThermostat or qMass.\n"
183 < //                "    Neither of these was found in the BASS file.\n");
184 < //       painCave.isFatal = 1;
185 < //       simError();
186 < //     }
187 <
188 <  strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
189 <  simnfo->usePBC = the_globals->getPBC();
190 <          
191 <  int usesDipoles = 0;
192 <  switch( ffCase ){
193 <
194 <  case FF_DUFF:
195 <    the_ff = new DUFF();
196 <    usesDipoles = 1;
197 <    break;
198 <
199 <  case FF_LJ:
200 <    the_ff = new LJFF();
201 <    break;
202 <
203 <  default:
204 <    sprintf( painCave.errMsg,
205 <             "SimSetup Error. Unrecognized force field in case statement.\n");
206 <    painCave.isFatal = 1;
207 <    simError();
208 <  }
209 <
210 < #ifdef IS_MPI
211 <  strcpy( checkPointMsg, "ForceField creation successful" );
212 <  MPIcheckPoint();
213 < #endif // is_mpi
214 <
215 <  // get the components and calculate the tot_nMol and indvidual n_mol
216 <  the_components = the_globals->getComponents();
217 <  components_nmol = new int[n_components];
218 <  comp_stamps = new MoleculeStamp*[n_components];
219 <
220 <  if( !the_globals->haveNMol() ){
221 <    // we don't have the total number of molecules, so we assume it is
222 <    // given in each component
223 <
224 <    tot_nmol = 0;
225 <    for( i=0; i<n_components; i++ ){
226 <
227 <      if( !the_components[i]->haveNMol() ){
228 <        // we have a problem
229 <        sprintf( painCave.errMsg,
230 <                 "SimSetup Error. No global NMol or component NMol"
231 <                 " given. Cannot calculate the number of atoms.\n" );
232 <        painCave.isFatal = 1;
233 <        simError();
234 <      }
235 <
236 <      tot_nmol += the_components[i]->getNMol();
237 <      components_nmol[i] = the_components[i]->getNMol();
238 <    }
239 <  }
240 <  else{
241 <    sprintf( painCave.errMsg,
242 <             "SimSetup error.\n"
243 <             "\tSorry, the ability to specify total"
244 <             " nMols and then give molfractions in the components\n"
245 <             "\tis not currently supported."
246 <             " Please give nMol in the components.\n" );
247 <    painCave.isFatal = 1;
248 <    simError();
249 <    
250 <    
251 <    //     tot_nmol = the_globals->getNMol();
252 <    
253 <    //   //we have the total number of molecules, now we check for molfractions
254 <    //     for( i=0; i<n_components; i++ ){
255 <    
256 <    //       if( !the_components[i]->haveMolFraction() ){
257 <    
258 <    //  if( !the_components[i]->haveNMol() ){
259 <    //    //we have a problem
260 <    //    std::cerr << "SimSetup error. Neither molFraction nor "
261 <    //              << " nMol was given in component
262 <    
263 <  }
264 <
265 < #ifdef IS_MPI
266 <  strcpy( checkPointMsg, "Have the number of components" );
267 <  MPIcheckPoint();
268 < #endif // is_mpi
269 <
270 <  // make an array of molecule stamps that match the components used.
271 <  // also extract the used stamps out into a separate linked list
272 <
273 <  simnfo->nComponents = n_components;
274 <  simnfo->componentsNmol = components_nmol;
275 <  simnfo->compStamps = comp_stamps;
276 <  simnfo->headStamp = new LinkedMolStamp();
91 >  // check on the post processing info
92    
93 <  char* id;
279 <  LinkedMolStamp* headStamp = simnfo->headStamp;
280 <  LinkedMolStamp* currentStamp = NULL;
281 <  for( i=0; i<n_components; i++ ){
93 >  finalInfoCheck();
94  
95 <    id = the_components[i]->getType();
284 <    comp_stamps[i] = NULL;
285 <    
286 <    // check to make sure the component isn't already in the list
95 >  // initialize the system coordinates
96  
97 <    comp_stamps[i] = headStamp->match( id );
289 <    if( comp_stamps[i] == NULL ){
290 <      
291 <      // extract the component from the list;
292 <      
293 <      currentStamp = the_stamps->extractMolStamp( id );
294 <      if( currentStamp == NULL ){
295 <        sprintf( painCave.errMsg,
296 <                 "SimSetup error: Component \"%s\" was not found in the "
297 <                 "list of declared molecules\n",
298 <                 id );
299 <        painCave.isFatal = 1;
300 <        simError();
301 <      }
302 <      
303 <      headStamp->add( currentStamp );
304 <      comp_stamps[i] = headStamp->match( id );
305 <    }
306 <  }
307 <
308 < #ifdef IS_MPI
309 <  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
310 <  MPIcheckPoint();
311 < #endif // is_mpi
97 >  initSystemCoords();
98    
99  
100 <
315 <
316 <  // caclulate the number of atoms, bonds, bends and torsions
317 <
318 <  tot_atoms = 0;
319 <  tot_bonds = 0;
320 <  tot_bends = 0;
321 <  tot_torsions = 0;
322 <  for( i=0; i<n_components; i++ ){
323 <    
324 <    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
325 <    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
326 <    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
327 <    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
328 <  }
329 <
330 <  tot_SRI = tot_bonds + tot_bends + tot_torsions;
331 <
332 <  simnfo->n_atoms = tot_atoms;
333 <  simnfo->n_bonds = tot_bonds;
334 <  simnfo->n_bends = tot_bends;
335 <  simnfo->n_torsions = tot_torsions;
336 <  simnfo->n_SRI = tot_SRI;
337 <  simnfo->n_mol = tot_nmol;
338 <  
339 <  simnfo->molMembershipArray = new int[tot_atoms];
340 <
341 < #ifdef IS_MPI
342 <
343 <  // divide the molecules among processors here.
344 <  
345 <  mpiSim = new mpiSimulation( simnfo );
346 <  
347 <  globalIndex = mpiSim->divideLabor();
348 <
349 <  // set up the local variables
350 <  
351 <  int localMol, allMol;
352 <  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
353 <
354 <  int* mol2proc = mpiSim->getMolToProcMap();
355 <  int* molCompType = mpiSim->getMolComponentType();
356 <  
357 <  allMol = 0;
358 <  localMol = 0;
359 <  local_atoms = 0;
360 <  local_bonds = 0;
361 <  local_bends = 0;
362 <  local_torsions = 0;
363 <  globalAtomIndex = 0;
364 <
365 <
366 <  for( i=0; i<n_components; i++ ){
367 <
368 <    for( j=0; j<components_nmol[i]; j++ ){
369 <      
370 <      if( mol2proc[allMol] == worldRank ){
371 <        
372 <        local_atoms +=    comp_stamps[i]->getNAtoms();
373 <        local_bonds +=    comp_stamps[i]->getNBonds();
374 <        local_bends +=    comp_stamps[i]->getNBends();
375 <        local_torsions += comp_stamps[i]->getNTorsions();
376 <        localMol++;
377 <      }      
378 <      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
379 <        simnfo->molMembershipArray[globalAtomIndex] = allMol;
380 <        globalAtomIndex++;
381 <      }
382 <
383 <      allMol++;      
384 <    }
385 <  }
386 <  local_SRI = local_bonds + local_bends + local_torsions;
387 <  
388 <  simnfo->n_atoms = mpiSim->getMyNlocal();  
389 <  
390 <  if( local_atoms != simnfo->n_atoms ){
391 <    sprintf( painCave.errMsg,
392 <             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
393 <             " localAtom (%d) are not equal.\n",
394 <             simnfo->n_atoms,
395 <             local_atoms );
396 <    painCave.isFatal = 1;
397 <    simError();
398 <  }
399 <
400 <  simnfo->n_bonds = local_bonds;
401 <  simnfo->n_bends = local_bends;
402 <  simnfo->n_torsions = local_torsions;
403 <  simnfo->n_SRI = local_SRI;
404 <  simnfo->n_mol = localMol;
405 <
406 <  strcpy( checkPointMsg, "Passed nlocal consistency check." );
407 <  MPIcheckPoint();
408 <  
409 <  
410 < #endif // is_mpi
411 <  
412 <
413 <  // create the atom and short range interaction arrays
414 <
415 <  Atom::createArrays(simnfo->n_atoms);
416 <  the_atoms = new Atom*[simnfo->n_atoms];
417 <  the_molecules = new Molecule[simnfo->n_mol];
418 <  int molIndex;
419 <
420 <  // initialize the molecule's stampID's
421 <
422 < #ifdef IS_MPI
423 <  
424 <
425 <  molIndex = 0;
426 <  for(i=0; i<mpiSim->getTotNmol(); i++){
427 <    
428 <    if(mol2proc[i] == worldRank ){
429 <      the_molecules[molIndex].setStampID( molCompType[i] );
430 <      the_molecules[molIndex].setMyIndex( molIndex );
431 <      the_molecules[molIndex].setGlobalIndex( i );
432 <      molIndex++;
433 <    }
434 <  }
435 <
436 < #else // is_mpi
437 <  
438 <  molIndex = 0;
439 <  globalAtomIndex = 0;
440 <  for(i=0; i<n_components; i++){
441 <    for(j=0; j<components_nmol[i]; j++ ){
442 <      the_molecules[molIndex].setStampID( i );
443 <      the_molecules[molIndex].setMyIndex( molIndex );
444 <      the_molecules[molIndex].setGlobalIndex( molIndex );
445 <      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
446 <        simnfo->molMembershipArray[globalAtomIndex] = molIndex;
447 <        globalAtomIndex++;
448 <      }
449 <      molIndex++;
450 <    }
451 <  }
452 <    
453 <
454 < #endif // is_mpi
455 <
456 <
457 <  if( simnfo->n_SRI ){
458 <    
459 <    Exclude::createArray(simnfo->n_SRI);
460 <    the_excludes = new Exclude*[simnfo->n_SRI];
461 <    for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
462 <    simnfo->globalExcludes = new int;
463 <    simnfo->n_exclude = simnfo->n_SRI;
464 <  }
465 <  else{
466 <    
467 <    Exclude::createArray( 1 );
468 <    the_excludes = new Exclude*;
469 <    the_excludes[0] = new Exclude(0);
470 <    the_excludes[0]->setPair( 0,0 );
471 <    simnfo->globalExcludes = new int;
472 <    simnfo->globalExcludes[0] = 0;
473 <    simnfo->n_exclude = 0;
474 <  }
475 <
476 <  // set the arrays into the SimInfo object
477 <
478 <  simnfo->atoms = the_atoms;
479 <  simnfo->molecules = the_molecules;
480 <  simnfo->nGlobalExcludes = 0;
481 <  simnfo->excludes = the_excludes;
482 <
483 <
484 <  // get some of the tricky things that may still be in the globals
485 <
486 <  double boxVector[3];
487 <  if( the_globals->haveBox() ){
488 <    boxVector[0] = the_globals->getBox();
489 <    boxVector[1] = the_globals->getBox();
490 <    boxVector[2] = the_globals->getBox();
491 <    
492 <    simnfo->setBox( boxVector );
493 <  }
494 <  else if( the_globals->haveDensity() ){
495 <
496 <    double vol;
497 <    vol = (double)tot_nmol / the_globals->getDensity();
498 <     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
499 <     boxVector[1] = boxVector[0];
500 <     boxVector[2] = boxVector[0];
501 <
502 <    simnfo->setBox( boxVector );
503 <  }
504 <  else{
505 <    if( !the_globals->haveBoxX() ){
506 <      sprintf( painCave.errMsg,
507 <               "SimSetup error, no periodic BoxX size given.\n" );
508 <      painCave.isFatal = 1;
509 <      simError();
510 <    }
511 <    boxVector[0] = the_globals->getBoxX();
512 <
513 <    if( !the_globals->haveBoxY() ){
514 <      sprintf( painCave.errMsg,
515 <               "SimSetup error, no periodic BoxY size given.\n" );
516 <      painCave.isFatal = 1;
517 <      simError();
518 <    }
519 <    boxVector[1] = the_globals->getBoxY();
520 <
521 <    if( !the_globals->haveBoxZ() ){
522 <      sprintf( painCave.errMsg,
523 <               "SimSetup error, no periodic BoxZ size given.\n" );
524 <      painCave.isFatal = 1;
525 <      simError();
526 <    }
527 <    boxVector[2] = the_globals->getBoxZ();
528 <
529 <    simnfo->setBox( boxVector );
530 <  }
531 <
532 < #ifdef IS_MPI
533 <  strcpy( checkPointMsg, "Box size set up" );
534 <  MPIcheckPoint();
535 < #endif // is_mpi
100 >  // make the output filenames
101  
102 <
538 <  // initialize the arrays
539 <
540 <  the_ff->setSimInfo( simnfo );
541 <
542 <  makeMolecules();
543 <  simnfo->identArray = new int[simnfo->n_atoms];
544 <  for(i=0; i<simnfo->n_atoms; i++){
545 <    simnfo->identArray[i] = the_atoms[i]->getIdent();
546 <  }
547 <  
548 <  if (the_globals->getUseRF() ) {
549 <    simnfo->useReactionField = 1;
550 <  
551 <    if( !the_globals->haveECR() ){
552 <      sprintf( painCave.errMsg,
553 <               "SimSetup Warning: using default value of 1/2 the smallest "
554 <               "box length for the electrostaticCutoffRadius.\n"
555 <               "I hope you have a very fast processor!\n");
556 <      painCave.isFatal = 0;
557 <      simError();
558 <      double smallest;
559 <      smallest = simnfo->boxLx;
560 <      if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
561 <      if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
562 <      simnfo->ecr = 0.5 * smallest;
563 <    } else {
564 <      simnfo->ecr        = the_globals->getECR();
565 <    }
566 <
567 <    if( !the_globals->haveEST() ){
568 <      sprintf( painCave.errMsg,
569 <               "SimSetup Warning: using default value of 0.05 * the "
570 <               "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
571 <               );
572 <      painCave.isFatal = 0;
573 <      simError();
574 <      simnfo->est = 0.05 * simnfo->ecr;
575 <    } else {
576 <      simnfo->est        = the_globals->getEST();
577 <    }
578 <    
579 <    if(!the_globals->haveDielectric() ){
580 <      sprintf( painCave.errMsg,
581 <               "SimSetup Error: You are trying to use Reaction Field without"
582 <               "setting a dielectric constant!\n"
583 <               );
584 <      painCave.isFatal = 1;
585 <      simError();
586 <    }
587 <    simnfo->dielectric = the_globals->getDielectric();  
588 <  } else {
589 <    if (usesDipoles) {
590 <      
591 <      if( !the_globals->haveECR() ){
592 <        sprintf( painCave.errMsg,
593 <                 "SimSetup Warning: using default value of 1/2 the smallest "
594 <                 "box length for the electrostaticCutoffRadius.\n"
595 <                 "I hope you have a very fast processor!\n");
596 <        painCave.isFatal = 0;
597 <        simError();
598 <        double smallest;
599 <        smallest = simnfo->boxLx;
600 <        if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
601 <        if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
602 <        simnfo->ecr = 0.5 * smallest;
603 <      } else {
604 <        simnfo->ecr        = the_globals->getECR();
605 <      }
606 <      
607 <      if( !the_globals->haveEST() ){
608 <        sprintf( painCave.errMsg,
609 <                 "SimSetup Warning: using default value of 5%% of the "
610 <                 "electrostaticCutoffRadius for the "
611 <                 "electrostaticSkinThickness\n"
612 <                 );
613 <        painCave.isFatal = 0;
614 <        simError();
615 <        simnfo->est = 0.05 * simnfo->ecr;
616 <      } else {
617 <        simnfo->est        = the_globals->getEST();
618 <      }
619 <    }
620 <  }  
621 <
622 < #ifdef IS_MPI
623 <  strcpy( checkPointMsg, "electrostatic parameters check out" );
624 <  MPIcheckPoint();
625 < #endif // is_mpi
626 <
627 < if( the_globals->haveInitialConfig() ){
628 <
629 <     InitializeFromFile* fileInit;
630 < #ifdef IS_MPI // is_mpi
631 <     if( worldRank == 0 ){
632 < #endif //is_mpi
633 <   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
634 < #ifdef IS_MPI
635 <     }else fileInit = new InitializeFromFile( NULL );
636 < #endif
637 <   fileInit->read_xyz( simnfo ); // default velocities on
638 <
639 <   delete fileInit;
640 < }
641 < else{
642 <
643 < #ifdef IS_MPI
644 <
645 <  // no init from bass
646 <  
647 <  sprintf( painCave.errMsg,
648 <           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
649 <  painCave.isFatal;
650 <  simError();
651 <  
652 < #else
653 <
654 <  initFromBass();
655 <
656 <
657 < #endif
658 < }
659 <
660 < #ifdef IS_MPI
661 <  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
662 <  MPIcheckPoint();
663 < #endif // is_mpi
664 <
665 <
102 >  makeOutNames();
103    
667 –
668 –  
669 –
670 –  
671 – #ifdef IS_MPI
672 –  if( worldRank == 0 ){
673 – #endif // is_mpi
674 –    
675 –    if( the_globals->haveFinalConfig() ){
676 –      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
677 –    }
678 –    else{
679 –      strcpy( simnfo->finalName, inFileName );
680 –      char* endTest;
681 –      int nameLength = strlen( simnfo->finalName );
682 –      endTest = &(simnfo->finalName[nameLength - 5]);
683 –      if( !strcmp( endTest, ".bass" ) ){
684 –        strcpy( endTest, ".eor" );
685 –      }
686 –      else if( !strcmp( endTest, ".BASS" ) ){
687 –        strcpy( endTest, ".eor" );
688 –      }
689 –      else{
690 –        endTest = &(simnfo->finalName[nameLength - 4]);
691 –        if( !strcmp( endTest, ".bss" ) ){
692 –          strcpy( endTest, ".eor" );
693 –        }
694 –        else if( !strcmp( endTest, ".mdl" ) ){
695 –          strcpy( endTest, ".eor" );
696 –        }
697 –        else{
698 –          strcat( simnfo->finalName, ".eor" );
699 –        }
700 –      }
701 –    }
702 –    
703 –    // make the sample and status out names
704 –    
705 –    strcpy( simnfo->sampleName, inFileName );
706 –    char* endTest;
707 –    int nameLength = strlen( simnfo->sampleName );
708 –    endTest = &(simnfo->sampleName[nameLength - 5]);
709 –    if( !strcmp( endTest, ".bass" ) ){
710 –      strcpy( endTest, ".dump" );
711 –    }
712 –    else if( !strcmp( endTest, ".BASS" ) ){
713 –      strcpy( endTest, ".dump" );
714 –    }
715 –    else{
716 –      endTest = &(simnfo->sampleName[nameLength - 4]);
717 –      if( !strcmp( endTest, ".bss" ) ){
718 –        strcpy( endTest, ".dump" );
719 –      }
720 –      else if( !strcmp( endTest, ".mdl" ) ){
721 –        strcpy( endTest, ".dump" );
722 –      }
723 –      else{
724 –        strcat( simnfo->sampleName, ".dump" );
725 –      }
726 –    }
727 –    
728 –    strcpy( simnfo->statusName, inFileName );
729 –    nameLength = strlen( simnfo->statusName );
730 –    endTest = &(simnfo->statusName[nameLength - 5]);
731 –    if( !strcmp( endTest, ".bass" ) ){
732 –      strcpy( endTest, ".stat" );
733 –    }
734 –    else if( !strcmp( endTest, ".BASS" ) ){
735 –      strcpy( endTest, ".stat" );
736 –    }
737 –    else{
738 –      endTest = &(simnfo->statusName[nameLength - 4]);
739 –      if( !strcmp( endTest, ".bss" ) ){
740 –        strcpy( endTest, ".stat" );
741 –      }
742 –      else if( !strcmp( endTest, ".mdl" ) ){
743 –        strcpy( endTest, ".stat" );
744 –      }
745 –      else{
746 –        strcat( simnfo->statusName, ".stat" );
747 –      }
748 –    }
749 –    
750 – #ifdef IS_MPI
751 –  }
752 – #endif // is_mpi
753 –  
754 –  // set the status, sample, and themal kick times
755 –  
756 –  if( the_globals->haveSampleTime() ){
757 –    simnfo->sampleTime = the_globals->getSampleTime();
758 –    simnfo->statusTime = simnfo->sampleTime;
759 –    simnfo->thermalTime = simnfo->sampleTime;
760 –  }
761 –  else{
762 –    simnfo->sampleTime = the_globals->getRunTime();
763 –    simnfo->statusTime = simnfo->sampleTime;
764 –    simnfo->thermalTime = simnfo->sampleTime;
765 –  }
766 –
767 –  if( the_globals->haveStatusTime() ){
768 –    simnfo->statusTime = the_globals->getStatusTime();
769 –  }
770 –
771 –  if( the_globals->haveThermalTime() ){
772 –    simnfo->thermalTime = the_globals->getThermalTime();
773 –  }
774 –
775 –  // check for the temperature set flag
776 –
777 –  if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
778 –
779 –
104    // make the integrator
105    
106 +  makeIntegrator();
107    
783 –  NVT*  myNVT = NULL;
784 –  NPTi* myNPTi = NULL;
785 –  NPTf* myNPTf = NULL;
786 –  switch( ensembleCase ){
787 –
788 –  case NVE_ENS:
789 –    new NVE( simnfo, the_ff );
790 –    break;
791 –
792 –  case NVT_ENS:
793 –    myNVT = new NVT( simnfo, the_ff );
794 –    myNVT->setTargetTemp(the_globals->getTargetTemp());
795 –
796 –    if (the_globals->haveTauThermostat())
797 –      myNVT->setTauThermostat(the_globals->getTauThermostat());
798 –
799 –    else {
800 –      sprintf( painCave.errMsg,
801 –               "SimSetup error: If you use the NVT\n"
802 –               "    ensemble, you must set tauThermostat.\n");
803 –      painCave.isFatal = 1;
804 –      simError();
805 –    }
806 –    break;
807 –
808 –  case NPTi_ENS:
809 –    myNPTi = new NPTi( simnfo, the_ff );
810 –    myNPTi->setTargetTemp( the_globals->getTargetTemp());
811 –
812 –    if (the_globals->haveTargetPressure())
813 –      myNPTi->setTargetPressure(the_globals->getTargetPressure());
814 –    else {
815 –      sprintf( painCave.errMsg,
816 –               "SimSetup error: If you use a constant pressure\n"
817 –               "    ensemble, you must set targetPressure in the BASS file.\n");
818 –      painCave.isFatal = 1;
819 –      simError();
820 –    }
821 –    
822 –    if( the_globals->haveTauThermostat() )
823 –      myNPTi->setTauThermostat( the_globals->getTauThermostat() );
824 –    else{
825 –      sprintf( painCave.errMsg,
826 –               "SimSetup error: If you use an NPT\n"
827 –               "    ensemble, you must set tauThermostat.\n");
828 –      painCave.isFatal = 1;
829 –      simError();
830 –    }
831 –
832 –    if( the_globals->haveTauBarostat() )
833 –      myNPTi->setTauBarostat( the_globals->getTauBarostat() );
834 –    else{
835 –      sprintf( painCave.errMsg,
836 –               "SimSetup error: If you use an NPT\n"
837 –               "    ensemble, you must set tauBarostat.\n");
838 –      painCave.isFatal = 1;
839 –      simError();
840 –    }
841 –    break;
842 –
843 –  case NPTf_ENS:
844 –    myNPTf = new NPTf( simnfo, the_ff );
845 –    myNPTf->setTargetTemp( the_globals->getTargetTemp());
846 –
847 –    if (the_globals->haveTargetPressure())
848 –      myNPTf->setTargetPressure(the_globals->getTargetPressure());
849 –    else {
850 –      sprintf( painCave.errMsg,
851 –               "SimSetup error: If you use a constant pressure\n"
852 –               "    ensemble, you must set targetPressure in the BASS file.\n");
853 –      painCave.isFatal = 1;
854 –      simError();
855 –    }    
856 –
857 –    if( the_globals->haveTauThermostat() )
858 –      myNPTf->setTauThermostat( the_globals->getTauThermostat() );
859 –    else{
860 –      sprintf( painCave.errMsg,
861 –               "SimSetup error: If you use an NPT\n"
862 –               "    ensemble, you must set tauThermostat.\n");
863 –      painCave.isFatal = 1;
864 –      simError();
865 –    }
866 –
867 –    if( the_globals->haveTauBarostat() )
868 –      myNPTf->setTauBarostat( the_globals->getTauBarostat() );
869 –    else{
870 –      sprintf( painCave.errMsg,
871 –               "SimSetup error: If you use an NPT\n"
872 –               "    ensemble, you must set tauBarostat.\n");
873 –      painCave.isFatal = 1;
874 –      simError();
875 –    }
876 –    break;
877 –    
878 –    
879 –
880 –  default:
881 –    sprintf( painCave.errMsg,
882 –             "SimSetup Error. Unrecognized ensemble in case statement.\n");
883 –    painCave.isFatal = 1;
884 –    simError();
885 –  }
886 –
887 –
108   #ifdef IS_MPI
109    mpiSim->mpiRefresh();
110   #endif
111  
112    // initialize the Fortran
113  
114 +  initFortran();
115  
895 –  simnfo->refreshSim();
896 –  
897 –  if( !strcmp( simnfo->mixingRule, "standard") ){
898 –    the_ff->initForceField( LB_MIXING_RULE );
899 –  }
900 –  else if( !strcmp( simnfo->mixingRule, "explicit") ){
901 –    the_ff->initForceField( EXPLICIT_MIXING_RULE );
902 –  }
903 –  else{
904 –    sprintf( painCave.errMsg,
905 –             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
906 –             simnfo->mixingRule );
907 –    painCave.isFatal = 1;
908 –    simError();
909 –  }
116  
117  
912 – #ifdef IS_MPI
913 –  strcpy( checkPointMsg,
914 –          "Successfully intialized the mixingRule for Fortran." );
915 –  MPIcheckPoint();
916 – #endif // is_mpi
118   }
119  
120  
121   void SimSetup::makeMolecules( void ){
122  
123    int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
124 <  molInit info;
124 >  molInit molInfo;
125    DirectionalAtom* dAtom;
126    LinkedAssign* extras;
127    LinkedAssign* current_extra;
# Line 945 | Line 146 | void SimSetup::makeMolecules( void ){
146    
147    atomOffset = 0;
148    excludeOffset = 0;
149 <  for(i=0; i<simnfo->n_mol; i++){
149 >  for(i=0; i<info->n_mol; i++){
150      
151      stampID = the_molecules[i].getStampID();
152  
153 <    info.nAtoms    = comp_stamps[stampID]->getNAtoms();
154 <    info.nBonds    = comp_stamps[stampID]->getNBonds();
155 <    info.nBends    = comp_stamps[stampID]->getNBends();
156 <    info.nTorsions = comp_stamps[stampID]->getNTorsions();
157 <    info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
153 >    molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
154 >    molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
155 >    molInfo.nBends    = comp_stamps[stampID]->getNBends();
156 >    molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
157 >    molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
158  
159 <    info.myAtoms = &the_atoms[atomOffset];
160 <    info.myExcludes = &the_excludes[excludeOffset];
161 <    info.myBonds = new Bond*[info.nBonds];
162 <    info.myBends = new Bend*[info.nBends];
163 <    info.myTorsions = new Torsion*[info.nTorsions];
159 >    molInfo.myAtoms = &the_atoms[atomOffset];
160 >    molInfo.myExcludes = &the_excludes[excludeOffset];
161 >    molInfo.myBonds = new Bond*[molInfo.nBonds];
162 >    molInfo.myBends = new Bend*[molInfo.nBends];
163 >    molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
164  
165 <    theBonds = new bond_pair[info.nBonds];
166 <    theBends = new bend_set[info.nBends];
167 <    theTorsions = new torsion_set[info.nTorsions];
165 >    theBonds = new bond_pair[molInfo.nBonds];
166 >    theBends = new bend_set[molInfo.nBends];
167 >    theTorsions = new torsion_set[molInfo.nTorsions];
168      
169      // make the Atoms
170      
171 <    for(j=0; j<info.nAtoms; j++){
171 >    for(j=0; j<molInfo.nAtoms; j++){
172        
173        currentAtom = comp_stamps[stampID]->getAtom( j );
174        if( currentAtom->haveOrientation() ){
175          
176          dAtom = new DirectionalAtom(j + atomOffset);
177 <        simnfo->n_oriented++;
178 <        info.myAtoms[j] = dAtom;
177 >        info->n_oriented++;
178 >        molInfo.myAtoms[j] = dAtom;
179          
180          ux = currentAtom->getOrntX();
181          uy = currentAtom->getOrntY();
# Line 992 | Line 193 | void SimSetup::makeMolecules( void ){
193          dAtom->setSUz( uz );
194        }
195        else{
196 <        info.myAtoms[j] = new GeneralAtom(j + atomOffset);
196 >        molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset);
197        }
198 <      info.myAtoms[j]->setType( currentAtom->getType() );
198 >      molInfo.myAtoms[j]->setType( currentAtom->getType() );
199      
200   #ifdef IS_MPI
201        
202 <      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
202 >      molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
203        
204   #endif // is_mpi
205      }
206      
207      // make the bonds
208 <    for(j=0; j<info.nBonds; j++){
208 >    for(j=0; j<molInfo.nBonds; j++){
209        
210        currentBond = comp_stamps[stampID]->getBond( j );
211        theBonds[j].a = currentBond->getA() + atomOffset;
# Line 1031 | Line 232 | void SimSetup::makeMolecules( void ){
232        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
233   #endif  //is_mpi
234      }
235 <    excludeOffset += info.nBonds;
235 >    excludeOffset += molInfo.nBonds;
236  
237      //make the bends
238 <    for(j=0; j<info.nBends; j++){
238 >    for(j=0; j<molInfo.nBends; j++){
239        
240        currentBend = comp_stamps[stampID]->getBend( j );
241        theBends[j].a = currentBend->getA() + atomOffset;
# Line 1117 | Line 318 | void SimSetup::makeMolecules( void ){
318        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
319   #endif  //is_mpi
320      }
321 <    excludeOffset += info.nBends;
321 >    excludeOffset += molInfo.nBends;
322  
323 <    for(j=0; j<info.nTorsions; j++){
323 >    for(j=0; j<molInfo.nTorsions; j++){
324        
325        currentTorsion = comp_stamps[stampID]->getTorsion( j );
326        theTorsions[j].a = currentTorsion->getA() + atomOffset;
# Line 1147 | Line 348 | void SimSetup::makeMolecules( void ){
348        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
349   #endif  //is_mpi
350      }
351 <    excludeOffset += info.nTorsions;
351 >    excludeOffset += molInfo.nTorsions;
352  
353      
354      // send the arrays off to the forceField for init.
355  
356 <    the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
357 <    the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
358 <    the_ff->initializeBends( info.nBends, info.myBends, theBends );
359 <    the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
356 >    the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
357 >    the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
358 >    the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
359 >    the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
360  
361  
362 <    the_molecules[i].initialize( info );
362 >    the_molecules[i].initialize( molInfo );
363  
364  
365 <    atomOffset += info.nAtoms;
365 >    atomOffset += molInfo.nAtoms;
366      delete[] theBonds;
367      delete[] theBends;
368      delete[] theTorsions;
# Line 1197 | Line 398 | void SimSetup::initFromBass( void ){
398      have_extra =1;
399  
400      n_cells = (int)temp3 - 1;
401 <    cellx = simnfo->boxLx / temp3;
402 <    celly = simnfo->boxLy / temp3;
403 <    cellz = simnfo->boxLz / temp3;
401 >    cellx = info->boxL[0] / temp3;
402 >    celly = info->boxL[1] / temp3;
403 >    cellz = info->boxL[2] / temp3;
404      n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
405      temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
406      n_per_extra = (int)ceil( temp1 );
# Line 1214 | Line 415 | void SimSetup::initFromBass( void ){
415    }
416    else{
417      n_cells = (int)temp3;
418 <    cellx = simnfo->boxLx / temp3;
419 <    celly = simnfo->boxLy / temp3;
420 <    cellz = simnfo->boxLz / temp3;
418 >    cellx = info->boxL[0] / temp3;
419 >    celly = info->boxL[1] / temp3;
420 >    cellz = info->boxL[2] / temp3;
421    }
422  
423    current_mol = 0;
# Line 1296 | Line 497 | void SimSetup::initFromBass( void ){
497    }
498  
499  
500 <  for( i=0; i<simnfo->n_atoms; i++ ){
501 <    simnfo->atoms[i]->set_vx( 0.0 );
502 <    simnfo->atoms[i]->set_vy( 0.0 );
503 <    simnfo->atoms[i]->set_vz( 0.0 );
500 >  for( i=0; i<info->n_atoms; i++ ){
501 >    info->atoms[i]->set_vx( 0.0 );
502 >    info->atoms[i]->set_vy( 0.0 );
503 >    info->atoms[i]->set_vz( 0.0 );
504    }
505   }
506  
# Line 1360 | Line 561 | void SimSetup::makeElement( double x, double y, double
561      current_comp++;
562    }
563   }
564 +
565 +
566 + void SimSetup::gatherInfo( void ){
567 +  int i,j,k;
568 +
569 +  ensembleCase = -1;
570 +  ffCase = -1;
571 +
572 +  // get the stamps and globals;
573 +  stamps = stamps;
574 +  globals = globals;
575 +
576 +  // set the easy ones first
577 +  info->target_temp = globals->getTargetTemp();
578 +  info->dt = globals->getDt();
579 +  info->run_time = globals->getRunTime();
580 +  n_components = globals->getNComponents();
581 +
582 +
583 +  // get the forceField
584 +
585 +  strcpy( force_field, globals->getForceField() );
586 +
587 +  if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
588 +  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
589 +  else{
590 +    sprintf( painCave.errMsg,
591 +             "SimSetup Error. Unrecognized force field -> %s\n",
592 +             force_field );
593 +    painCave.isFatal = 1;
594 +    simError();
595 +  }
596 +
597 +  // get the ensemble
598 +
599 +  strcpy( ensemble, globals->getEnsemble() );
600 +
601 +  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
602 +  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
603 +  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") )
604 +    ensembleCase = NPTi_ENS;
605 +  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
606 +  else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
607 +  else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
608 +  else{
609 +    sprintf( painCave.errMsg,
610 +             "SimSetup Warning. Unrecognized Ensemble -> %s, "
611 +             "reverting to NVE for this simulation.\n",
612 +             ensemble );
613 +    painCave.isFatal = 0;
614 +    simError();
615 +    strcpy( ensemble, "NVE" );
616 +    ensembleCase = NVE_ENS;
617 +  }  
618 +  strcpy( info->ensemble, ensemble );
619 +
620 +  // get the mixing rule
621 +
622 +  strcpy( info->mixingRule, globals->getMixingRule() );
623 +  info->usePBC = globals->getPBC();
624 +        
625 +  
626 +  // get the components and calculate the tot_nMol and indvidual n_mol
627 +
628 +  the_components = globals->getComponents();
629 +  components_nmol = new int[n_components];
630 +
631 +
632 +  if( !globals->haveNMol() ){
633 +    // we don't have the total number of molecules, so we assume it is
634 +    // given in each component
635 +
636 +    tot_nmol = 0;
637 +    for( i=0; i<n_components; i++ ){
638 +
639 +      if( !the_components[i]->haveNMol() ){
640 +        // we have a problem
641 +        sprintf( painCave.errMsg,
642 +                 "SimSetup Error. No global NMol or component NMol"
643 +                 " given. Cannot calculate the number of atoms.\n" );
644 +        painCave.isFatal = 1;
645 +        simError();
646 +      }
647 +
648 +      tot_nmol += the_components[i]->getNMol();
649 +      components_nmol[i] = the_components[i]->getNMol();
650 +    }
651 +  }
652 +  else{
653 +    sprintf( painCave.errMsg,
654 +             "SimSetup error.\n"
655 +             "\tSorry, the ability to specify total"
656 +             " nMols and then give molfractions in the components\n"
657 +             "\tis not currently supported."
658 +             " Please give nMol in the components.\n" );
659 +    painCave.isFatal = 1;
660 +    simError();
661 +  }
662 +
663 +  // set the status, sample, and thermal kick times
664 +  
665 +  if( globals->haveSampleTime() ){
666 +    info->sampleTime = globals->getSampleTime();
667 +    info->statusTime = info->sampleTime;
668 +    info->thermalTime = info->sampleTime;
669 +  }
670 +  else{
671 +    info->sampleTime = globals->getRunTime();
672 +    info->statusTime = info->sampleTime;
673 +    info->thermalTime = info->sampleTime;
674 +  }
675 +
676 +  if( globals->haveStatusTime() ){
677 +    info->statusTime = globals->getStatusTime();
678 +  }
679 +
680 +  if( globals->haveThermalTime() ){
681 +    info->thermalTime = globals->getThermalTime();
682 +  }
683 +
684 +  // check for the temperature set flag
685 +
686 +  if( globals->haveTempSet() ) info->setTemp = globals->getTempSet();
687 +
688 +  // get some of the tricky things that may still be in the globals
689 +
690 +  double boxVector[3];
691 +  if( globals->haveBox() ){
692 +    boxVector[0] = globals->getBox();
693 +    boxVector[1] = globals->getBox();
694 +    boxVector[2] = globals->getBox();
695 +    
696 +    info->setBox( boxVector );
697 +  }
698 +  else if( globals->haveDensity() ){
699 +
700 +    double vol;
701 +    vol = (double)tot_nmol / globals->getDensity();
702 +     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
703 +     boxVector[1] = boxVector[0];
704 +     boxVector[2] = boxVector[0];
705 +
706 +    info->setBox( boxVector );
707 +  }
708 +  else{
709 +    if( !globals->haveBoxX() ){
710 +      sprintf( painCave.errMsg,
711 +               "SimSetup error, no periodic BoxX size given.\n" );
712 +      painCave.isFatal = 1;
713 +      simError();
714 +    }
715 +    boxVector[0] = globals->getBoxX();
716 +
717 +    if( !globals->haveBoxY() ){
718 +      sprintf( painCave.errMsg,
719 +               "SimSetup error, no periodic BoxY size given.\n" );
720 +      painCave.isFatal = 1;
721 +      simError();
722 +    }
723 +    boxVector[1] = globals->getBoxY();
724 +
725 +    if( !globals->haveBoxZ() ){
726 +      sprintf( painCave.errMsg,
727 +               "SimSetup error, no periodic BoxZ size given.\n" );
728 +      painCave.isFatal = 1;
729 +      simError();
730 +    }
731 +    boxVector[2] = globals->getBoxZ();
732 +
733 +    info->setBox( boxVector );
734 +  }
735 +
736 +
737 +    
738 + #ifdef IS_MPI
739 +  strcpy( checkPointMsg, "Succesfully gathered all information from Bass\n" );
740 +  MPIcheckPoint();
741 + #endif // is_mpi
742 +
743 + }
744 +
745 +
746 + void SimSetup::finalInfoCheck( void ){
747 +  int index;
748 +  int usesDipoles;
749 +  
750 +
751 +  // check electrostatic parameters
752 +  
753 +  index = 0;
754 +  usesDipoles = 0;
755 +  while( (index < info->n_atoms) && !usesDipoles ){
756 +    usesDipoles = ((info->atoms)[index])->hasDipole();
757 +    index++;
758 +  }
759 +  
760 + #ifdef IS_MPI
761 +  int myUse = usesDipoles;
762 +  MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
763 + #endif //is_mpi
764 +
765 +
766 +  if (globals->getUseRF() ) {
767 +    info->useReactionField = 1;
768 +    
769 +    if( !globals->haveECR() ){
770 +      sprintf( painCave.errMsg,
771 +               "SimSetup Warning: using default value of 1/2 the smallest "
772 +               "box length for the electrostaticCutoffRadius.\n"
773 +               "I hope you have a very fast processor!\n");
774 +      painCave.isFatal = 0;
775 +      simError();
776 +      double smallest;
777 +      smallest = info->boxL[0];
778 +      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
779 +      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
780 +      info->ecr = 0.5 * smallest;
781 +    } else {
782 +      info->ecr        = globals->getECR();
783 +    }
784 +
785 +    if( !globals->haveEST() ){
786 +      sprintf( painCave.errMsg,
787 +               "SimSetup Warning: using default value of 0.05 * the "
788 +               "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
789 +               );
790 +      painCave.isFatal = 0;
791 +      simError();
792 +      info->est = 0.05 * info->ecr;
793 +    } else {
794 +      info->est        = globals->getEST();
795 +    }
796 +    
797 +    if(!globals->haveDielectric() ){
798 +      sprintf( painCave.errMsg,
799 +               "SimSetup Error: You are trying to use Reaction Field without"
800 +               "setting a dielectric constant!\n"
801 +               );
802 +      painCave.isFatal = 1;
803 +      simError();
804 +    }
805 +    info->dielectric = globals->getDielectric();  
806 +  }
807 +  else {
808 +    if (usesDipoles) {
809 +      
810 +      if( !globals->haveECR() ){
811 +        sprintf( painCave.errMsg,
812 +                 "SimSetup Warning: using default value of 1/2 the smallest "
813 +                 "box length for the electrostaticCutoffRadius.\n"
814 +                 "I hope you have a very fast processor!\n");
815 +        painCave.isFatal = 0;
816 +        simError();
817 +        double smallest;
818 +        smallest = info->boxL[0];
819 +        if (info->boxL[1] <= smallest) smallest = info->boxL[1];
820 +        if (info->boxL[2] <= smallest) smallest = info->boxL[2];
821 +        info->ecr = 0.5 * smallest;
822 +      } else {
823 +        info->ecr        = globals->getECR();
824 +      }
825 +      
826 +      if( !globals->haveEST() ){
827 +        sprintf( painCave.errMsg,
828 +                 "SimSetup Warning: using default value of 5%% of the "
829 +                 "electrostaticCutoffRadius for the "
830 +                 "electrostaticSkinThickness\n"
831 +                 );
832 +        painCave.isFatal = 0;
833 +        simError();
834 +        info->est = 0.05 * info->ecr;
835 +      } else {
836 +        info->est        = globals->getEST();
837 +      }
838 +    }
839 +  }  
840 +
841 + #ifdef IS_MPI
842 +  strcpy( checkPointMsg, "post processing checks out" );
843 +  MPIcheckPoint();
844 + #endif // is_mpi
845 +
846 + }
847 +
848 + void SimSetup::initSystemCoords( void ){
849 +
850 + if( globals->haveInitialConfig() ){
851 +
852 +     InitializeFromFile* fileInit;
853 + #ifdef IS_MPI // is_mpi
854 +     if( worldRank == 0 ){
855 + #endif //is_mpi
856 +   fileInit = new InitializeFromFile( globals->getInitialConfig() );
857 + #ifdef IS_MPI
858 +     }else fileInit = new InitializeFromFile( NULL );
859 + #endif
860 +   fileInit->read_xyz( info ); // default velocities on
861 +
862 +   delete fileInit;
863 + }
864 + else{
865 +
866 + #ifdef IS_MPI
867 +
868 +  // no init from bass
869 +  
870 +  sprintf( painCave.errMsg,
871 +           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
872 +  painCave.isFatal;
873 +  simError();
874 +  
875 + #else
876 +
877 +  initFromBass();
878 +
879 +
880 + #endif
881 + }
882 +
883 + #ifdef IS_MPI
884 +  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
885 +  MPIcheckPoint();
886 + #endif // is_mpi
887 +
888 + }
889 +
890 +
891 + void SimSetup::makeOutNames( void ){
892 +
893 + #ifdef IS_MPI
894 +  if( worldRank == 0 ){
895 + #endif // is_mpi
896 +    
897 +    if( globals->haveFinalConfig() ){
898 +      strcpy( info->finalName, globals->getFinalConfig() );
899 +    }
900 +    else{
901 +      strcpy( info->finalName, inFileName );
902 +      char* endTest;
903 +      int nameLength = strlen( info->finalName );
904 +      endTest = &(info->finalName[nameLength - 5]);
905 +      if( !strcmp( endTest, ".bass" ) ){
906 +        strcpy( endTest, ".eor" );
907 +      }
908 +      else if( !strcmp( endTest, ".BASS" ) ){
909 +        strcpy( endTest, ".eor" );
910 +      }
911 +      else{
912 +        endTest = &(info->finalName[nameLength - 4]);
913 +        if( !strcmp( endTest, ".bss" ) ){
914 +          strcpy( endTest, ".eor" );
915 +        }
916 +        else if( !strcmp( endTest, ".mdl" ) ){
917 +          strcpy( endTest, ".eor" );
918 +        }
919 +        else{
920 +          strcat( info->finalName, ".eor" );
921 +        }
922 +      }
923 +    }
924 +    
925 +    // make the sample and status out names
926 +    
927 +    strcpy( info->sampleName, inFileName );
928 +    char* endTest;
929 +    int nameLength = strlen( info->sampleName );
930 +    endTest = &(info->sampleName[nameLength - 5]);
931 +    if( !strcmp( endTest, ".bass" ) ){
932 +      strcpy( endTest, ".dump" );
933 +    }
934 +    else if( !strcmp( endTest, ".BASS" ) ){
935 +      strcpy( endTest, ".dump" );
936 +    }
937 +    else{
938 +      endTest = &(info->sampleName[nameLength - 4]);
939 +      if( !strcmp( endTest, ".bss" ) ){
940 +        strcpy( endTest, ".dump" );
941 +      }
942 +      else if( !strcmp( endTest, ".mdl" ) ){
943 +        strcpy( endTest, ".dump" );
944 +      }
945 +      else{
946 +        strcat( info->sampleName, ".dump" );
947 +      }
948 +    }
949 +    
950 +    strcpy( info->statusName, inFileName );
951 +    nameLength = strlen( info->statusName );
952 +    endTest = &(info->statusName[nameLength - 5]);
953 +    if( !strcmp( endTest, ".bass" ) ){
954 +      strcpy( endTest, ".stat" );
955 +    }
956 +    else if( !strcmp( endTest, ".BASS" ) ){
957 +      strcpy( endTest, ".stat" );
958 +    }
959 +    else{
960 +      endTest = &(info->statusName[nameLength - 4]);
961 +      if( !strcmp( endTest, ".bss" ) ){
962 +        strcpy( endTest, ".stat" );
963 +      }
964 +      else if( !strcmp( endTest, ".mdl" ) ){
965 +        strcpy( endTest, ".stat" );
966 +      }
967 +      else{
968 +        strcat( info->statusName, ".stat" );
969 +      }
970 +    }
971 +    
972 + #ifdef IS_MPI
973 +  }
974 + #endif // is_mpi
975 +
976 + }
977 +
978 +
979 + void SimSetup::sysObjectsCreation( void ){
980 +
981 +  int i;
982 +
983 +  // create the forceField
984 +
985 +  createFF();
986 +
987 +  // extract componentList
988 +
989 +  compList();
990 +
991 +  // calc the number of atoms, bond, bends, and torsions
992 +
993 +  calcSysValues();
994 +
995 + #ifdef IS_MPI
996 +  // divide the molecules among the processors
997 +  
998 +  mpiMolDivide();
999 + #endif //is_mpi
1000 +  
1001 +  // create the atom and SRI arrays. Also initialize Molecule Stamp ID's
1002 +  
1003 +  makeSysArrays();
1004 +
1005 +  // make and initialize the molecules (all but atomic coordinates)
1006 +  
1007 +  makeMolecules();
1008 +  info->identArray = new int[info->n_atoms];
1009 +  for(i=0; i<info->n_atoms; i++){
1010 +    info->identArray[i] = the_atoms[i]->getIdent();
1011 +  }
1012 +  
1013 +
1014 +
1015 + }
1016 +
1017 +
1018 + void SimSetup::createFF( void ){
1019 +
1020 +  switch( ffCase ){
1021 +
1022 +  case FF_DUFF:
1023 +    the_ff = new DUFF();
1024 +    break;
1025 +
1026 +  case FF_LJ:
1027 +    the_ff = new LJFF();
1028 +    break;
1029 +
1030 +  default:
1031 +    sprintf( painCave.errMsg,
1032 +             "SimSetup Error. Unrecognized force field in case statement.\n");
1033 +    painCave.isFatal = 1;
1034 +    simError();
1035 +  }
1036 +
1037 + #ifdef IS_MPI
1038 +  strcpy( checkPointMsg, "ForceField creation successful" );
1039 +  MPIcheckPoint();
1040 + #endif // is_mpi
1041 +
1042 + }
1043 +
1044 +
1045 + void SimSetup::compList( void ){
1046 +
1047 +  int i;
1048 +
1049 +  comp_stamps = new MoleculeStamp*[n_components];
1050 +
1051 +  // make an array of molecule stamps that match the components used.
1052 +  // also extract the used stamps out into a separate linked list
1053 +
1054 +  info->nComponents = n_components;
1055 +  info->componentsNmol = components_nmol;
1056 +  info->compStamps = comp_stamps;
1057 +  info->headStamp = new LinkedMolStamp();
1058 +  
1059 +  char* id;
1060 +  LinkedMolStamp* headStamp = info->headStamp;
1061 +  LinkedMolStamp* currentStamp = NULL;
1062 +  for( i=0; i<n_components; i++ ){
1063 +
1064 +    id = the_components[i]->getType();
1065 +    comp_stamps[i] = NULL;
1066 +    
1067 +    // check to make sure the component isn't already in the list
1068 +
1069 +    comp_stamps[i] = headStamp->match( id );
1070 +    if( comp_stamps[i] == NULL ){
1071 +      
1072 +      // extract the component from the list;
1073 +      
1074 +      currentStamp = stamps->extractMolStamp( id );
1075 +      if( currentStamp == NULL ){
1076 +        sprintf( painCave.errMsg,
1077 +                 "SimSetup error: Component \"%s\" was not found in the "
1078 +                 "list of declared molecules\n",
1079 +                 id );
1080 +        painCave.isFatal = 1;
1081 +        simError();
1082 +      }
1083 +      
1084 +      headStamp->add( currentStamp );
1085 +      comp_stamps[i] = headStamp->match( id );
1086 +    }
1087 +  }
1088 +
1089 + #ifdef IS_MPI
1090 +  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
1091 +  MPIcheckPoint();
1092 + #endif // is_mpi
1093 +
1094 +
1095 + }
1096 +
1097 + void SimSetup::calcSysValues( void ){
1098 +  int i, j, k;
1099 +
1100 +
1101 +  tot_atoms = 0;
1102 +  tot_bonds = 0;
1103 +  tot_bends = 0;
1104 +  tot_torsions = 0;
1105 +  for( i=0; i<n_components; i++ ){
1106 +    
1107 +    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
1108 +    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
1109 +    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
1110 +    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1111 +  }
1112 +
1113 +  tot_SRI = tot_bonds + tot_bends + tot_torsions;
1114 +
1115 +  info->n_atoms = tot_atoms;
1116 +  info->n_bonds = tot_bonds;
1117 +  info->n_bends = tot_bends;
1118 +  info->n_torsions = tot_torsions;
1119 +  info->n_SRI = tot_SRI;
1120 +  info->n_mol = tot_nmol;
1121 +  
1122 +  info->molMembershipArray = new int[tot_atoms];
1123 + }
1124 +
1125 +
1126 + #ifdef IS_MPI
1127 +
1128 + void SimSetup::mpiMolDivide( void ){
1129 +  
1130 +  int i, j, k;
1131 +  int localMol, allMol;
1132 +  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1133 +
1134 +  mpiSim = new mpiSimulation( info );
1135 +  
1136 +  globalIndex = mpiSim->divideLabor();
1137 +
1138 +  // set up the local variables
1139 +  
1140 +  mol2proc = mpiSim->getMolToProcMap();
1141 +  molCompType = mpiSim->getMolComponentType();
1142 +  
1143 +  allMol = 0;
1144 +  localMol = 0;
1145 +  local_atoms = 0;
1146 +  local_bonds = 0;
1147 +  local_bends = 0;
1148 +  local_torsions = 0;
1149 +  globalAtomIndex = 0;
1150 +
1151 +
1152 +  for( i=0; i<n_components; i++ ){
1153 +
1154 +    for( j=0; j<components_nmol[i]; j++ ){
1155 +      
1156 +      if( mol2proc[allMol] == worldRank ){
1157 +        
1158 +        local_atoms +=    comp_stamps[i]->getNAtoms();
1159 +        local_bonds +=    comp_stamps[i]->getNBonds();
1160 +        local_bends +=    comp_stamps[i]->getNBends();
1161 +        local_torsions += comp_stamps[i]->getNTorsions();
1162 +        localMol++;
1163 +      }      
1164 +      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
1165 +        info->molMembershipArray[globalAtomIndex] = allMol;
1166 +        globalAtomIndex++;
1167 +      }
1168 +
1169 +      allMol++;      
1170 +    }
1171 +  }
1172 +  local_SRI = local_bonds + local_bends + local_torsions;
1173 +  
1174 +  info->n_atoms = mpiSim->getMyNlocal();  
1175 +  
1176 +  if( local_atoms != info->n_atoms ){
1177 +    sprintf( painCave.errMsg,
1178 +             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1179 +             " localAtom (%d) are not equal.\n",
1180 +             info->n_atoms,
1181 +             local_atoms );
1182 +    painCave.isFatal = 1;
1183 +    simError();
1184 +  }
1185 +
1186 +  info->n_bonds = local_bonds;
1187 +  info->n_bends = local_bends;
1188 +  info->n_torsions = local_torsions;
1189 +  info->n_SRI = local_SRI;
1190 +  info->n_mol = localMol;
1191 +
1192 +  strcpy( checkPointMsg, "Passed nlocal consistency check." );
1193 +  MPIcheckPoint();
1194 + }
1195 +  
1196 + #endif // is_mpi
1197 +
1198 +
1199 + void SimSetup::makeSysArrays( void ){
1200 +  int i, j, k;
1201 +
1202 +
1203 +  // create the atom and short range interaction arrays
1204 +
1205 +  Atom::createArrays(info->n_atoms);
1206 +  the_atoms = new Atom*[info->n_atoms];
1207 +  the_molecules = new Molecule[info->n_mol];
1208 +  int molIndex;
1209 +
1210 +  // initialize the molecule's stampID's
1211 +
1212 + #ifdef IS_MPI
1213 +  
1214 +
1215 +  molIndex = 0;
1216 +  for(i=0; i<mpiSim->getTotNmol(); i++){
1217 +    
1218 +    if(mol2proc[i] == worldRank ){
1219 +      the_molecules[molIndex].setStampID( molCompType[i] );
1220 +      the_molecules[molIndex].setMyIndex( molIndex );
1221 +      the_molecules[molIndex].setGlobalIndex( i );
1222 +      molIndex++;
1223 +    }
1224 +  }
1225 +
1226 + #else // is_mpi
1227 +  
1228 +  molIndex = 0;
1229 +  globalAtomIndex = 0;
1230 +  for(i=0; i<n_components; i++){
1231 +    for(j=0; j<components_nmol[i]; j++ ){
1232 +      the_molecules[molIndex].setStampID( i );
1233 +      the_molecules[molIndex].setMyIndex( molIndex );
1234 +      the_molecules[molIndex].setGlobalIndex( molIndex );
1235 +      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
1236 +        info->molMembershipArray[globalAtomIndex] = molIndex;
1237 +        globalAtomIndex++;
1238 +      }
1239 +      molIndex++;
1240 +    }
1241 +  }
1242 +    
1243 +
1244 + #endif // is_mpi
1245 +
1246 +
1247 +  if( info->n_SRI ){
1248 +    
1249 +    Exclude::createArray(info->n_SRI);
1250 +    the_excludes = new Exclude*[info->n_SRI];
1251 +    for( int ex=0; ex<info->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
1252 +    info->globalExcludes = new int;
1253 +    info->n_exclude = info->n_SRI;
1254 +  }
1255 +  else{
1256 +    
1257 +    Exclude::createArray( 1 );
1258 +    the_excludes = new Exclude*;
1259 +    the_excludes[0] = new Exclude(0);
1260 +    the_excludes[0]->setPair( 0,0 );
1261 +    info->globalExcludes = new int;
1262 +    info->globalExcludes[0] = 0;
1263 +    info->n_exclude = 0;
1264 +  }
1265 +
1266 +  // set the arrays into the SimInfo object
1267 +
1268 +  info->atoms = the_atoms;
1269 +  info->molecules = the_molecules;
1270 +  info->nGlobalExcludes = 0;
1271 +  info->excludes = the_excludes;
1272 +
1273 +  the_ff->setSimInfo( info );
1274 +
1275 + }
1276 +
1277 + void SimSetup::makeIntegrator( void ){
1278 +
1279 +  NVT*  myNVT = NULL;
1280 +  NPTi* myNPTi = NULL;
1281 +  NPTf* myNPTf = NULL;
1282 +  NPTim* myNPTim = NULL;
1283 +  NPTfm* myNPTfm = NULL;
1284 +
1285 +  switch( ensembleCase ){
1286 +
1287 +  case NVE_ENS:
1288 +    new NVE( info, the_ff );
1289 +    break;
1290 +
1291 +  case NVT_ENS:
1292 +    myNVT = new NVT( info, the_ff );
1293 +    myNVT->setTargetTemp(globals->getTargetTemp());
1294 +
1295 +    if (globals->haveTauThermostat())
1296 +      myNVT->setTauThermostat(globals->getTauThermostat());
1297 +
1298 +    else {
1299 +      sprintf( painCave.errMsg,
1300 +               "SimSetup error: If you use the NVT\n"
1301 +               "    ensemble, you must set tauThermostat.\n");
1302 +      painCave.isFatal = 1;
1303 +      simError();
1304 +    }
1305 +    break;
1306 +
1307 +  case NPTi_ENS:
1308 +    myNPTi = new NPTi( info, the_ff );
1309 +    myNPTi->setTargetTemp( globals->getTargetTemp() );
1310 +
1311 +    if (globals->haveTargetPressure())
1312 +      myNPTi->setTargetPressure(globals->getTargetPressure());
1313 +    else {
1314 +      sprintf( painCave.errMsg,
1315 +               "SimSetup error: If you use a constant pressure\n"
1316 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1317 +      painCave.isFatal = 1;
1318 +      simError();
1319 +    }
1320 +    
1321 +    if( globals->haveTauThermostat() )
1322 +      myNPTi->setTauThermostat( globals->getTauThermostat() );
1323 +    else{
1324 +      sprintf( painCave.errMsg,
1325 +               "SimSetup error: If you use an NPT\n"
1326 +               "    ensemble, you must set tauThermostat.\n");
1327 +      painCave.isFatal = 1;
1328 +      simError();
1329 +    }
1330 +
1331 +    if( globals->haveTauBarostat() )
1332 +      myNPTi->setTauBarostat( globals->getTauBarostat() );
1333 +    else{
1334 +      sprintf( painCave.errMsg,
1335 +               "SimSetup error: If you use an NPT\n"
1336 +               "    ensemble, you must set tauBarostat.\n");
1337 +      painCave.isFatal = 1;
1338 +      simError();
1339 +    }
1340 +    break;
1341 +
1342 +  case NPTf_ENS:
1343 +    myNPTf = new NPTf( info, the_ff );
1344 +    myNPTf->setTargetTemp( globals->getTargetTemp());
1345 +
1346 +    if (globals->haveTargetPressure())
1347 +      myNPTf->setTargetPressure(globals->getTargetPressure());
1348 +    else {
1349 +      sprintf( painCave.errMsg,
1350 +               "SimSetup error: If you use a constant pressure\n"
1351 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1352 +      painCave.isFatal = 1;
1353 +      simError();
1354 +    }    
1355 +
1356 +    if( globals->haveTauThermostat() )
1357 +      myNPTf->setTauThermostat( globals->getTauThermostat() );
1358 +    else{
1359 +      sprintf( painCave.errMsg,
1360 +               "SimSetup error: If you use an NPT\n"
1361 +               "    ensemble, you must set tauThermostat.\n");
1362 +      painCave.isFatal = 1;
1363 +      simError();
1364 +    }
1365 +
1366 +    if( globals->haveTauBarostat() )
1367 +      myNPTf->setTauBarostat( globals->getTauBarostat() );
1368 +    else{
1369 +      sprintf( painCave.errMsg,
1370 +               "SimSetup error: If you use an NPT\n"
1371 +               "    ensemble, you must set tauBarostat.\n");
1372 +      painCave.isFatal = 1;
1373 +      simError();
1374 +    }
1375 +    break;
1376 +    
1377 +  case NPTim_ENS:
1378 +    myNPTim = new NPTim( info, the_ff );
1379 +    myNPTim->setTargetTemp( globals->getTargetTemp());
1380 +
1381 +    if (globals->haveTargetPressure())
1382 +      myNPTim->setTargetPressure(globals->getTargetPressure());
1383 +    else {
1384 +      sprintf( painCave.errMsg,
1385 +               "SimSetup error: If you use a constant pressure\n"
1386 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1387 +      painCave.isFatal = 1;
1388 +      simError();
1389 +    }
1390 +    
1391 +    if( globals->haveTauThermostat() )
1392 +      myNPTim->setTauThermostat( globals->getTauThermostat() );
1393 +    else{
1394 +      sprintf( painCave.errMsg,
1395 +               "SimSetup error: If you use an NPT\n"
1396 +               "    ensemble, you must set tauThermostat.\n");
1397 +      painCave.isFatal = 1;
1398 +      simError();
1399 +    }
1400 +
1401 +    if( globals->haveTauBarostat() )
1402 +      myNPTim->setTauBarostat( globals->getTauBarostat() );
1403 +    else{
1404 +      sprintf( painCave.errMsg,
1405 +               "SimSetup error: If you use an NPT\n"
1406 +               "    ensemble, you must set tauBarostat.\n");
1407 +      painCave.isFatal = 1;
1408 +      simError();
1409 +    }
1410 +    break;
1411 +
1412 +  case NPTfm_ENS:
1413 +    myNPTfm = new NPTfm( info, the_ff );
1414 +    myNPTfm->setTargetTemp( globals->getTargetTemp());
1415 +
1416 +    if (globals->haveTargetPressure())
1417 +      myNPTfm->setTargetPressure(globals->getTargetPressure());
1418 +    else {
1419 +      sprintf( painCave.errMsg,
1420 +               "SimSetup error: If you use a constant pressure\n"
1421 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1422 +      painCave.isFatal = 1;
1423 +      simError();
1424 +    }
1425 +    
1426 +    if( globals->haveTauThermostat() )
1427 +      myNPTfm->setTauThermostat( globals->getTauThermostat() );
1428 +    else{
1429 +      sprintf( painCave.errMsg,
1430 +               "SimSetup error: If you use an NPT\n"
1431 +               "    ensemble, you must set tauThermostat.\n");
1432 +      painCave.isFatal = 1;
1433 +      simError();
1434 +    }
1435 +
1436 +    if( globals->haveTauBarostat() )
1437 +      myNPTfm->setTauBarostat( globals->getTauBarostat() );
1438 +    else{
1439 +      sprintf( painCave.errMsg,
1440 +               "SimSetup error: If you use an NPT\n"
1441 +               "    ensemble, you must set tauBarostat.\n");
1442 +      painCave.isFatal = 1;
1443 +      simError();
1444 +    }
1445 +    break;
1446 +
1447 +  default:
1448 +    sprintf( painCave.errMsg,
1449 +             "SimSetup Error. Unrecognized ensemble in case statement.\n");
1450 +    painCave.isFatal = 1;
1451 +    simError();
1452 +  }
1453 +
1454 + }
1455 +
1456 + void SimSetup::initFortran( void ){
1457 +
1458 +  info->refreshSim();
1459 +  
1460 +  if( !strcmp( info->mixingRule, "standard") ){
1461 +    the_ff->initForceField( LB_MIXING_RULE );
1462 +  }
1463 +  else if( !strcmp( info->mixingRule, "explicit") ){
1464 +    the_ff->initForceField( EXPLICIT_MIXING_RULE );
1465 +  }
1466 +  else{
1467 +    sprintf( painCave.errMsg,
1468 +             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1469 +             info->mixingRule );
1470 +    painCave.isFatal = 1;
1471 +    simError();
1472 +  }
1473 +
1474 +
1475 + #ifdef IS_MPI
1476 +  strcpy( checkPointMsg,
1477 +          "Successfully intialized the mixingRule for Fortran." );
1478 +  MPIcheckPoint();
1479 + #endif // is_mpi
1480 +
1481 + }

Diff Legend

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