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 414 by mmeineke, Wed Mar 26 22:02:36 2003 UTC vs.
Revision 614 by mmeineke, Tue Jul 15 17:57:04 2003 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines