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 424 by mmeineke, Thu Mar 27 20:36:16 2003 UTC vs.
Revision 616 by mmeineke, Tue Jul 15 18:52:16 2003 UTC

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

Diff Legend

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