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