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 604 by gezelter, Tue Jul 15 03:08:00 2003 UTC vs.
Revision 660 by tim, Thu Jul 31 19:59:34 2003 UTC

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

Diff Legend

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