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 481 by gezelter, Tue Apr 8 21:35:49 2003 UTC vs.
Revision 661 by tim, Fri Aug 1 16:18:13 2003 UTC

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

Diff Legend

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