[ViewVC] Diff of: group/trunk/OOPSE/libmdtools/SimSetup.cpp

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 605 by gezelter, Tue Jul 15 03:27:24 2003 UTC vs.
Revision 621 by gezelter, Wed Jul 16 02:11:02 2003 UTC

+void SimSetup::createSim( void ){
-–
+  MakeStamps *the_stamps;
-–
+  Globals* the_globals;
+  int i, j, k, globalAtomIndex;
-<
+  int ensembleCase;
-<
+  int ffCase;
->
+  // gather all of the information from the Bass file
-+
+  gatherInfo();
-+
-+
+  // creation of complex system objects
-+
-+
+  sysObjectsCreation();
-+
-+
+  // check on the post processing info
-+
-+
+  finalInfoCheck();
-+
-+
+  // initialize the system coordinates
-+
-+
+  initSystemCoords();
-+
-+
-+
+  // make the output filenames
-+
-+
+  makeOutNames();
-+
-+
+  // make the integrator
-+
-+
+  makeIntegrator();
-+
-+
+#ifdef IS_MPI
-+
+  mpiSim->mpiRefresh();
-+
+#endif
-+
-+
+  // initialize the Fortran
-+
-+
+  initFortran();
-+
-+
-+
-+
+}
-+
-+
-+
+void SimSetup::makeMolecules( void ){
-+
-+
+  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
-+
+  molInit molInfo;
-+
+  DirectionalAtom* dAtom;
-+
+  LinkedAssign* extras;
-+
+  LinkedAssign* current_extra;
-+
+  AtomStamp* currentAtom;
-+
+  BondStamp* currentBond;
-+
+  BendStamp* currentBend;
-+
+  TorsionStamp* currentTorsion;
-+
-+
+  bond_pair* theBonds;
-+
+  bend_set* theBends;
-+
+  torsion_set* theTorsions;
-+
-+
-+
+  //init the forceField paramters
-+
-+
+  the_ff->readParams();
-+
-+
-+
+  // init the atoms
-+
-+
+  double ux, uy, uz, u, uSqr;
-+
-+
+  atomOffset = 0;
-+
+  excludeOffset = 0;
-+
+  for(i=0; i<info->n_mol; i++){
-+
-+
+    stampID = the_molecules[i].getStampID();
-+
-+
+    molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
-+
+    molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
-+
+    molInfo.nBends    = comp_stamps[stampID]->getNBends();
-+
+    molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
-+
+    molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
-+
-+
+    molInfo.myAtoms = &the_atoms[atomOffset];
-+
+    molInfo.myExcludes = &the_excludes[excludeOffset];
-+
+    molInfo.myBonds = new Bond*[molInfo.nBonds];
-+
+    molInfo.myBends = new Bend*[molInfo.nBends];
-+
+    molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
-+
-+
+    theBonds = new bond_pair[molInfo.nBonds];
-+
+    theBends = new bend_set[molInfo.nBends];
-+
+    theTorsions = new torsion_set[molInfo.nTorsions];
-+
-+
+    // make the Atoms
-+
-+
+    for(j=0; j<molInfo.nAtoms; j++){
-+
-+
+      currentAtom = comp_stamps[stampID]->getAtom( j );
-+
+      if( currentAtom->haveOrientation() ){
-+
-+
+        dAtom = new DirectionalAtom(j + atomOffset);
-+
+        info->n_oriented++;
-+
+        molInfo.myAtoms[j] = dAtom;
-+
-+
+        ux = currentAtom->getOrntX();
-+
+        uy = currentAtom->getOrntY();
-+
+        uz = currentAtom->getOrntZ();
-+
-+
+        uSqr = (ux * ux) + (uy * uy) + (uz * uz);
-+
-+
+        u = sqrt( uSqr );
-+
+        ux = ux / u;
-+
+        uy = uy / u;
-+
+        uz = uz / u;
-+
-+
+        dAtom->setSUx( ux );
-+
+        dAtom->setSUy( uy );
-+
+        dAtom->setSUz( uz );
-+
+      }
-+
+      else{
-+
+        molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset);
-+
+      }
-+
+      molInfo.myAtoms[j]->setType( currentAtom->getType() );
-+
-+
+#ifdef IS_MPI
-+
-+
+      molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
-+
-+
+#endif // is_mpi
-+
+    }
-+
-+
+    // make the bonds
-+
+    for(j=0; j<molInfo.nBonds; j++){
-+
-+
+      currentBond = comp_stamps[stampID]->getBond( j );
-+
+      theBonds[j].a = currentBond->getA() + atomOffset;
-+
+      theBonds[j].b = currentBond->getB() + atomOffset;
-+
-+
+      exI = theBonds[j].a;
-+
+      exJ = theBonds[j].b;
-+
-+
+      // exclude_I must always be the smaller of the pair
-+
+      if( exI > exJ ){
-+
+        tempEx = exI;
-+
+        exI = exJ;
-+
+        exJ = tempEx;
-+
+      }
-+
+#ifdef IS_MPI
-+
+      tempEx = exI;
-+
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
-+
+      tempEx = exJ;
-+
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
-+
-+
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
-+
+#else  // isn't MPI
-+
-+
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
-+
+#endif  //is_mpi
-+
+    }
-+
+    excludeOffset += molInfo.nBonds;
-+
-+
+    //make the bends
-+
+    for(j=0; j<molInfo.nBends; j++){
-+
-+
+      currentBend = comp_stamps[stampID]->getBend( j );
-+
+      theBends[j].a = currentBend->getA() + atomOffset;
-+
+      theBends[j].b = currentBend->getB() + atomOffset;
-+
+      theBends[j].c = currentBend->getC() + atomOffset;
-+
-+
+      if( currentBend->haveExtras() ){
-+
-+
+        extras = currentBend->getExtras();
-+
+        current_extra = extras;
-+
-+
+        while( current_extra != NULL ){
-+
+          if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
-+
-+
+            switch( current_extra->getType() ){
-+
-+
+            case 0:
-+
+              theBends[j].ghost =
-+
+                current_extra->getInt() + atomOffset;
-+
+              theBends[j].isGhost = 1;
-+
+              break;
-+
-+
+            case 1:
-+
+              theBends[j].ghost =
-+
+                (int)current_extra->getDouble() + atomOffset;
-+
+              theBends[j].isGhost = 1;
-+
+              break;
-+
-+
+            default:
-+
+              sprintf( painCave.errMsg,
-+
+                       "SimSetup Error: ghostVectorSource was neither a "
-+
+                       "double nor an int.\n"
-+
+                       "-->Bend[%d] in %s\n",
-+
+                       j, comp_stamps[stampID]->getID() );
-+
+              painCave.isFatal = 1;
-+
+              simError();
-+
+            }
-+
+          }
-+
-+
+          else{
-+
-+
+            sprintf( painCave.errMsg,
-+
+                     "SimSetup Error: unhandled bend assignment:\n"
-+
+                     "    -->%s in Bend[%d] in %s\n",
-+
+                     current_extra->getlhs(),
-+
+                     j, comp_stamps[stampID]->getID() );
-+
+            painCave.isFatal = 1;
-+
+            simError();
-+
+          }
-+
-+
+          current_extra = current_extra->getNext();
-+
+        }
-+
+      }
-+
-+
+      if( !theBends[j].isGhost ){
-+
-+
+        exI = theBends[j].a;
-+
+        exJ = theBends[j].c;
-+
+      }
-+
+      else{
-+
-+
+        exI = theBends[j].a;
-+
+        exJ = theBends[j].b;
-+
+      }
-+
-+
+      // exclude_I must always be the smaller of the pair
-+
+      if( exI > exJ ){
-+
+        tempEx = exI;
-+
+        exI = exJ;
-+
+        exJ = tempEx;
-+
+      }
-+
+#ifdef IS_MPI
-+
+      tempEx = exI;
-+
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
-+
+      tempEx = exJ;
-+
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
-+
-+
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
-+
+#else  // isn't MPI
-+
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
-+
+#endif  //is_mpi
-+
+    }
-+
+    excludeOffset += molInfo.nBends;
-+
-+
+    for(j=0; j<molInfo.nTorsions; j++){
-+
-+
+      currentTorsion = comp_stamps[stampID]->getTorsion( j );
-+
+      theTorsions[j].a = currentTorsion->getA() + atomOffset;
-+
+      theTorsions[j].b = currentTorsion->getB() + atomOffset;
-+
+      theTorsions[j].c = currentTorsion->getC() + atomOffset;
-+
+      theTorsions[j].d = currentTorsion->getD() + atomOffset;
-+
-+
+      exI = theTorsions[j].a;
-+
+      exJ = theTorsions[j].d;
-+
-+
+      // exclude_I must always be the smaller of the pair
-+
+      if( exI > exJ ){
-+
+        tempEx = exI;
-+
+        exI = exJ;
-+
+        exJ = tempEx;
-+
+      }
-+
+#ifdef IS_MPI
-+
+      tempEx = exI;
-+
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
-+
+      tempEx = exJ;
-+
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
-+
-+
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
-+
+#else  // isn't MPI
-+
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
-+
+#endif  //is_mpi
-+
+    }
-+
+    excludeOffset += molInfo.nTorsions;
-+
-+
-+
+    // send the arrays off to the forceField for init.
-+
-+
+    the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
-+
+    the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
-+
+    the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
-+
+    the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
-+
-+
-+
+    the_molecules[i].initialize( molInfo );
-+
-+
-+
+    atomOffset += molInfo.nAtoms;
-+
+    delete[] theBonds;
-+
+    delete[] theBends;
-+
+    delete[] theTorsions;
-+
+  }
-+
-+
+#ifdef IS_MPI
-+
+  sprintf( checkPointMsg, "all molecules initialized succesfully" );
-+
+  MPIcheckPoint();
-+
+#endif // is_mpi
-+
-+
+  // clean up the forcefield
-+
+  the_ff->calcRcut();
-+
+  the_ff->cleanMe();
-+
-+
+}
-+
-+
+void SimSetup::initFromBass( void ){
-+
-+
+  int i, j, k;
-+
+  int n_cells;
-+
+  double cellx, celly, cellz;
-+
+  double temp1, temp2, temp3;
-+
+  int n_per_extra;
-+
+  int n_extra;
-+
+  int have_extra, done;
-+
-+
+  temp1 = (double)tot_nmol / 4.0;
-+
+  temp2 = pow( temp1, ( 1.0 / 3.0 ) );
-+
+  temp3 = ceil( temp2 );
-+
-+
+  have_extra =0;
-+
+  if( temp2 < temp3 ){ // we have a non-complete lattice
-+
+    have_extra =1;
-+
-+
+    n_cells = (int)temp3 - 1;
-+
+    cellx = info->boxL[0] / temp3;
-+
+    celly = info->boxL[1] / temp3;
-+
+    cellz = info->boxL[2] / temp3;
-+
+    n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
-+
+    temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
-+
+    n_per_extra = (int)ceil( temp1 );
-+
-+
+    if( n_per_extra > 4){
-+
+      sprintf( painCave.errMsg,
-+
+               "SimSetup error. There has been an error in constructing"
-+
+               " the non-complete lattice.\n" );
-+
+      painCave.isFatal = 1;
-+
+      simError();
-+
+    }
-+
+  }
-+
+  else{
-+
+    n_cells = (int)temp3;
-+
+    cellx = info->boxL[0] / temp3;
-+
+    celly = info->boxL[1] / temp3;
-+
+    cellz = info->boxL[2] / temp3;
-+
+  }
-+
-+
+  current_mol = 0;
-+
+  current_comp_mol = 0;
-+
+  current_comp = 0;
-+
+  current_atom_ndx = 0;
-+
-+
+  for( i=0; i < n_cells ; i++ ){
-+
+    for( j=0; j < n_cells; j++ ){
-+
+      for( k=0; k < n_cells; k++ ){
-+
-+
+        makeElement( i * cellx,
-+
+                     j * celly,
-+
+                     k * cellz );
-+
-+
+        makeElement( i * cellx + 0.5 * cellx,
-+
+                     j * celly + 0.5 * celly,
-+
+                     k * cellz );
-+
-+
+        makeElement( i * cellx,
-+
+                     j * celly + 0.5 * celly,
-+
+                     k * cellz + 0.5 * cellz );
-+
-+
+        makeElement( i * cellx + 0.5 * cellx,
-+
+                     j * celly,
-+
+                     k * cellz + 0.5 * cellz );
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  if( have_extra ){
-+
+    done = 0;
-+
-+
+    int start_ndx;
-+
+    for( i=0; i < (n_cells+1) && !done; i++ ){
-+
+      for( j=0; j < (n_cells+1) && !done; j++ ){
-+
-+
+        if( i < n_cells ){
-+
-+
+          if( j < n_cells ){
-+
+            start_ndx = n_cells;
-+
+          }
-+
+          else start_ndx = 0;
-+
+        }
-+
+        else start_ndx = 0;
-+
-+
+        for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
-+
-+
+          makeElement( i * cellx,
-+
+                       j * celly,
-+
+                       k * cellz );
-+
+          done = ( current_mol >= tot_nmol );
-+
-+
+          if( !done && n_per_extra > 1 ){
-+
+            makeElement( i * cellx + 0.5 * cellx,
-+
+                         j * celly + 0.5 * celly,
-+
+                         k * cellz );
-+
+            done = ( current_mol >= tot_nmol );
-+
+          }
-+
-+
+          if( !done && n_per_extra > 2){
-+
+            makeElement( i * cellx,
-+
+                         j * celly + 0.5 * celly,
-+
+                         k * cellz + 0.5 * cellz );
-+
+            done = ( current_mol >= tot_nmol );
-+
+          }
-+
-+
+          if( !done && n_per_extra > 3){
-+
+            makeElement( i * cellx + 0.5 * cellx,
-+
+                         j * celly,
-+
+                         k * cellz + 0.5 * cellz );
-+
+            done = ( current_mol >= tot_nmol );
-+
+          }
-+
+        }
-+
+      }
-+
+    }
-+
+  }
-+
-+
-+
+  for( i=0; i<info->n_atoms; i++ ){
-+
+    info->atoms[i]->set_vx( 0.0 );
-+
+    info->atoms[i]->set_vy( 0.0 );
-+
+    info->atoms[i]->set_vz( 0.0 );
-+
+  }
-+
+}
-+
-+
+void SimSetup::makeElement( double x, double y, double z ){
-+
-+
+  int k;
-+
+  AtomStamp* current_atom;
-+
+  DirectionalAtom* dAtom;
-+
+  double rotMat[3][3];
-+
-+
+  for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
-+
-+
+    current_atom = comp_stamps[current_comp]->getAtom( k );
-+
+    if( !current_atom->havePosition() ){
-+
+      sprintf( painCave.errMsg,
-+
+               "SimSetup:initFromBass error.\n"
-+
+               "\tComponent %s, atom %s does not have a position specified.\n"
-+
+               "\tThe initialization routine is unable to give a start"
-+
+               " position.\n",
-+
+               comp_stamps[current_comp]->getID(),
-+
+               current_atom->getType() );
-+
+      painCave.isFatal = 1;
-+
+      simError();
-+
+    }
-+
-+
+    the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
-+
+    the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
-+
+    the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
-+
-+
+    if( the_atoms[current_atom_ndx]->isDirectional() ){
-+
-+
+      dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
-+
-+
+      rotMat[0][0] = 1.0;
-+
+      rotMat[0][1] = 0.0;
-+
+      rotMat[0][2] = 0.0;
-+
-+
+      rotMat[1][0] = 0.0;
-+
+      rotMat[1][1] = 1.0;
-+
+      rotMat[1][2] = 0.0;
-+
-+
+      rotMat[2][0] = 0.0;
-+
+      rotMat[2][1] = 0.0;
-+
+      rotMat[2][2] = 1.0;
-+
-+
+      dAtom->setA( rotMat );
-+
+    }
-+
-+
+    current_atom_ndx++;
-+
+  }
-+
-+
+  current_mol++;
-+
+  current_comp_mol++;
-+
-+
+  if( current_comp_mol >= components_nmol[current_comp] ){
-+
-+
+    current_comp_mol = 0;
-+
+    current_comp++;
-+
+  }
-+
+}
-+
-+
-+
+void SimSetup::gatherInfo( void ){
-+
+  int i,j,k;
-+
+  ensembleCase = -1;
+  ffCase = -1;
+  // get the stamps and globals;
-<
+  the_stamps = stamps;
-<
+  the_globals = globals;
->
+  stamps = stamps;
->
+  globals = globals;
+  // set the easy ones first
-<
+  simnfo->target_temp = the_globals->getTargetTemp();
-<
+  simnfo->dt = the_globals->getDt();
-<
+  simnfo->run_time = the_globals->getRunTime();
->
+  info->target_temp = globals->getTargetTemp();
->
+  info->dt = globals->getDt();
->
+  info->run_time = globals->getRunTime();
->
+  n_components = globals->getNComponents();
-–
+  // get the ones we know are there, yet still may need some work.
-–
+  n_components = the_globals->getNComponents();
-–
+  strcpy( force_field, the_globals->getForceField() );
-+
+  // get the forceField
-+
-+
+  strcpy( force_field, globals->getForceField() );
-+
+  if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
+  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
+  else{
+    simError();
-<
+  // get the ensemble:
-<
+  strcpy( ensemble, the_globals->getEnsemble() );
->
+  // get the ensemble
-+
+  strcpy( ensemble, globals->getEnsemble() );
-+
+  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
+  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
+  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") )
+    strcpy( ensemble, "NVE" );
+    ensembleCase = NVE_ENS;
-<
+  strcpy( simnfo->ensemble, ensemble );
->
+  strcpy( info->ensemble, ensemble );
-+
+  // get the mixing rule
-<
+//   if( !strcasecmp( ensemble, "NPT" ) ) {
-<
+//     the_extendedsystem = new ExtendedSystem( simnfo );
-<
+//     the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
-<
+//     if (the_globals->haveTargetPressure())
-<
+//       the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
-<
+//     else {
-<
+//       sprintf( painCave.errMsg,
-<
+//                "SimSetup error: If you use the constant pressure\n"
-<
+//                "    ensemble, you must set targetPressure.\n"
-<
+//                "    This was found in the BASS file.\n");
-<
+//       painCave.isFatal = 1;
-<
+//       simError();
-<
+//     }
-<
-<
+//     if (the_globals->haveTauThermostat())
-<
+//       the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
-<
+//     else if (the_globals->haveQmass())
-<
+//       the_extendedsystem->setQmass(the_globals->getQmass());
-<
+//     else {
-<
+//       sprintf( painCave.errMsg,
-<
+//                "SimSetup error: If you use one of the constant temperature\n"
-<
+//                "    ensembles, you must set either tauThermostat or qMass.\n"
-<
+//                "    Neither of these was found in the BASS file.\n");
-<
+//       painCave.isFatal = 1;
-<
+//       simError();
-<
+//     }
-<
-<
+//     if (the_globals->haveTauBarostat())
-<
+//       the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
-<
+//     else {
-<
+//       sprintf( painCave.errMsg,
-<
+//                "SimSetup error: If you use the constant pressure\n"
-<
+//                "    ensemble, you must set tauBarostat.\n"
-<
+//                "    This was found in the BASS file.\n");
-<
+//       painCave.isFatal = 1;
-<
+//       simError();
-<
+//     }
-<
-<
+//   } else if ( !strcasecmp( ensemble, "NVT") ) {
-<
+//     the_extendedsystem = new ExtendedSystem( simnfo );
-<
+//     the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
-<
-<
+//     if (the_globals->haveTauThermostat())
-<
+//       the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
-<
+//     else if (the_globals->haveQmass())
-<
+//       the_extendedsystem->setQmass(the_globals->getQmass());
-<
+//     else {
-<
+//       sprintf( painCave.errMsg,
-<
+//                "SimSetup error: If you use one of the constant temperature\n"
-<
+//                "    ensembles, you must set either tauThermostat or qMass.\n"
-<
+//                "    Neither of these was found in the BASS file.\n");
-<
+//       painCave.isFatal = 1;
-<
+//       simError();
-<
+//     }
-<
-<
+  strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
-<
+  simnfo->usePBC = the_globals->getPBC();
-<
-<
+  int usesDipoles = 0;
-<
+  switch( ffCase ){
-<
-<
+  case FF_DUFF:
-<
+    the_ff = new DUFF();
-<
+    usesDipoles = 1;
-<
+    break;
-<
-<
+  case FF_LJ:
-<
+    the_ff = new LJFF();
-<
+    break;
-<
-<
+  default:
-<
+    sprintf( painCave.errMsg,
-<
+             "SimSetup Error. Unrecognized force field in case statement.\n");
-<
+    painCave.isFatal = 1;
-<
+    simError();
-<
+  }
-<
-<
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "ForceField creation successful" );
-<
+  MPIcheckPoint();
-<
+#endif // is_mpi
-<
->
+  strcpy( info->mixingRule, globals->getMixingRule() );
->
+  info->usePBC = globals->getPBC();
->
->
+  // get the components and calculate the tot_nMol and indvidual n_mol
-<
+  the_components = the_globals->getComponents();
->
->
+  the_components = globals->getComponents();
+  components_nmol = new int[n_components];
-–
+  comp_stamps = new MoleculeStamp*[n_components];
-<
+  if( !the_globals->haveNMol() ){
->
->
+  if( !globals->haveNMol() ){
+    // we don't have the total number of molecules, so we assume it is
+    // given in each component
+             " Please give nMol in the components.\n" );
+    painCave.isFatal = 1;
+    simError();
-–
-–
-–
+    //     tot_nmol = the_globals->getNMol();
-–
-–
+    //   //we have the total number of molecules, now we check for molfractions
-–
+    //     for( i=0; i<n_components; i++ ){
-–
-–
+    //       if( !the_components[i]->haveMolFraction() ){
-–
-–
+    //  if( !the_components[i]->haveNMol() ){
-–
+    //    //we have a problem
-–
+    //    std::cerr << "SimSetup error. Neither molFraction nor "
-–
+    //              << " nMol was given in component
-–
-<
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "Have the number of components" );
-<
+  MPIcheckPoint();
-<
+#endif // is_mpi
-<
-<
+  // make an array of molecule stamps that match the components used.
-<
+  // also extract the used stamps out into a separate linked list
-<
-<
+  simnfo->nComponents = n_components;
-<
+  simnfo->componentsNmol = components_nmol;
-<
+  simnfo->compStamps = comp_stamps;
-<
+  simnfo->headStamp = new LinkedMolStamp();
->
+  // set the status, sample, and thermal kick times
-<
+  char* id;
-<
+  LinkedMolStamp* headStamp = simnfo->headStamp;
-<
+  LinkedMolStamp* currentStamp = NULL;
-<
+  for( i=0; i<n_components; i++ ){
-<
-<
+    id = the_components[i]->getType();
-<
+    comp_stamps[i] = NULL;
-<
-<
+    // check to make sure the component isn't already in the list
-<
-<
+    comp_stamps[i] = headStamp->match( id );
-<
+    if( comp_stamps[i] == NULL ){
-<
-<
+      // extract the component from the list;
-<
-<
+      currentStamp = the_stamps->extractMolStamp( id );
-<
+      if( currentStamp == NULL ){
-<
+        sprintf( painCave.errMsg,
-<
+                 "SimSetup error: Component \"%s\" was not found in the "
-<
+                 "list of declared molecules\n",
-<
+                 id );
-<
+        painCave.isFatal = 1;
-<
+        simError();
-<
+      }
-<
-<
+      headStamp->add( currentStamp );
-<
+      comp_stamps[i] = headStamp->match( id );
-<
+    }
->
+  if( globals->haveSampleTime() ){
->
+    info->sampleTime = globals->getSampleTime();
->
+    info->statusTime = info->sampleTime;
->
+    info->thermalTime = info->sampleTime;
-+
+  else{
-+
+    info->sampleTime = globals->getRunTime();
-+
+    info->statusTime = info->sampleTime;
-+
+    info->thermalTime = info->sampleTime;
-+
+  }
-<
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
-<
+  MPIcheckPoint();
-<
+#endif // is_mpi
-<
-<
-<
-<
-<
+  // caclulate the number of atoms, bonds, bends and torsions
-<
-<
+  tot_atoms = 0;
-<
+  tot_bonds = 0;
-<
+  tot_bends = 0;
-<
+  tot_torsions = 0;
-<
+  for( i=0; i<n_components; i++ ){
-<
-<
+    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
-<
+    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
-<
+    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
-<
+    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
->
+  if( globals->haveStatusTime() ){
->
+    info->statusTime = globals->getStatusTime();
-<
+  tot_SRI = tot_bonds + tot_bends + tot_torsions;
-<
-<
+  simnfo->n_atoms = tot_atoms;
-<
+  simnfo->n_bonds = tot_bonds;
-<
+  simnfo->n_bends = tot_bends;
-<
+  simnfo->n_torsions = tot_torsions;
-<
+  simnfo->n_SRI = tot_SRI;
-<
+  simnfo->n_mol = tot_nmol;
-<
-<
+  simnfo->molMembershipArray = new int[tot_atoms];
-<
-<
+#ifdef IS_MPI
-<
-<
+  // divide the molecules among processors here.
-<
-<
+  mpiSim = new mpiSimulation( simnfo );
-<
-<
+  globalIndex = mpiSim->divideLabor();
-<
-<
+  // set up the local variables
-<
-<
+  int localMol, allMol;
-<
+  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
-<
-<
+  int* mol2proc = mpiSim->getMolToProcMap();
-<
+  int* molCompType = mpiSim->getMolComponentType();
-<
-<
+  allMol = 0;
-<
+  localMol = 0;
-<
+  local_atoms = 0;
-<
+  local_bonds = 0;
-<
+  local_bends = 0;
-<
+  local_torsions = 0;
-<
+  globalAtomIndex = 0;
->
+  if( globals->haveThermalTime() ){
->
+    info->thermalTime = globals->getThermalTime();
->
+  }
-+
+  // check for the temperature set flag
-<
+  for( i=0; i<n_components; i++ ){
->
+  if( globals->haveTempSet() ) info->setTemp = globals->getTempSet();
-–
+    for( j=0; j<components_nmol[i]; j++ ){
-–
-–
+      if( mol2proc[allMol] == worldRank ){
-–
-–
+        local_atoms +=    comp_stamps[i]->getNAtoms();
-–
+        local_bonds +=    comp_stamps[i]->getNBonds();
-–
+        local_bends +=    comp_stamps[i]->getNBends();
-–
+        local_torsions += comp_stamps[i]->getNTorsions();
-–
+        localMol++;
-–
+      }
-–
+      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
-–
+        simnfo->molMembershipArray[globalAtomIndex] = allMol;
-–
+        globalAtomIndex++;
-–
+      }
-–
-–
+      allMol++;
-–
+    }
-–
+  }
-–
+  local_SRI = local_bonds + local_bends + local_torsions;
-–
-–
+  simnfo->n_atoms = mpiSim->getMyNlocal();
-–
-–
+  if( local_atoms != simnfo->n_atoms ){
-–
+    sprintf( painCave.errMsg,
-–
+             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
-–
+             " localAtom (%d) are not equal.\n",
-–
+             simnfo->n_atoms,
-–
+             local_atoms );
-–
+    painCave.isFatal = 1;
-–
+    simError();
-–
+  }
-–
-–
+  simnfo->n_bonds = local_bonds;
-–
+  simnfo->n_bends = local_bends;
-–
+  simnfo->n_torsions = local_torsions;
-–
+  simnfo->n_SRI = local_SRI;
-–
+  simnfo->n_mol = localMol;
-–
-–
+  strcpy( checkPointMsg, "Passed nlocal consistency check." );
-–
+  MPIcheckPoint();
-–
-–
-–
+#endif // is_mpi
-–
-–
-–
+  // create the atom and short range interaction arrays
-–
-–
+  Atom::createArrays(simnfo->n_atoms);
-–
+  the_atoms = new Atom*[simnfo->n_atoms];
-–
+  the_molecules = new Molecule[simnfo->n_mol];
-–
+  int molIndex;
-–
-–
+  // initialize the molecule's stampID's
-–
-–
+#ifdef IS_MPI
-–
-–
-–
+  molIndex = 0;
-–
+  for(i=0; i<mpiSim->getTotNmol(); i++){
-–
-–
+    if(mol2proc[i] == worldRank ){
-–
+      the_molecules[molIndex].setStampID( molCompType[i] );
-–
+      the_molecules[molIndex].setMyIndex( molIndex );
-–
+      the_molecules[molIndex].setGlobalIndex( i );
-–
+      molIndex++;
-–
+    }
-–
+  }
-–
-–
+#else // is_mpi
-–
-–
+  molIndex = 0;
-–
+  globalAtomIndex = 0;
-–
+  for(i=0; i<n_components; i++){
-–
+    for(j=0; j<components_nmol[i]; j++ ){
-–
+      the_molecules[molIndex].setStampID( i );
-–
+      the_molecules[molIndex].setMyIndex( molIndex );
-–
+      the_molecules[molIndex].setGlobalIndex( molIndex );
-–
+      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
-–
+        simnfo->molMembershipArray[globalAtomIndex] = molIndex;
-–
+        globalAtomIndex++;
-–
+      }
-–
+      molIndex++;
-–
+    }
-–
+  }
-–
-–
-–
+#endif // is_mpi
-–
-–
-–
+  if( simnfo->n_SRI ){
-–
-–
+    Exclude::createArray(simnfo->n_SRI);
-–
+    the_excludes = new Exclude*[simnfo->n_SRI];
-–
+    for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
-–
+    simnfo->globalExcludes = new int;
-–
+    simnfo->n_exclude = simnfo->n_SRI;
-–
+  }
-–
+  else{
-–
-–
+    Exclude::createArray( 1 );
-–
+    the_excludes = new Exclude*;
-–
+    the_excludes[0] = new Exclude(0);
-–
+    the_excludes[0]->setPair( 0,0 );
-–
+    simnfo->globalExcludes = new int;
-–
+    simnfo->globalExcludes[0] = 0;
-–
+    simnfo->n_exclude = 0;
-–
+  }
-–
-–
+  // set the arrays into the SimInfo object
-–
-–
+  simnfo->atoms = the_atoms;
-–
+  simnfo->molecules = the_molecules;
-–
+  simnfo->nGlobalExcludes = 0;
-–
+  simnfo->excludes = the_excludes;
-–
-–
+  // get some of the tricky things that may still be in the globals
+  double boxVector[3];
-<
+  if( the_globals->haveBox() ){
-<
+    boxVector[0] = the_globals->getBox();
-<
+    boxVector[1] = the_globals->getBox();
-<
+    boxVector[2] = the_globals->getBox();
->
+  if( globals->haveBox() ){
->
+    boxVector[0] = globals->getBox();
->
+    boxVector[1] = globals->getBox();
->
+    boxVector[2] = globals->getBox();
-<
+    simnfo->setBox( boxVector );
->
+    info->setBox( boxVector );
-<
+  else if( the_globals->haveDensity() ){
->
+  else if( globals->haveDensity() ){
+    double vol;
-<
+    vol = (double)tot_nmol / the_globals->getDensity();
->
+    vol = (double)tot_nmol / globals->getDensity();
+     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
+     boxVector[1] = boxVector[0];
+     boxVector[2] = boxVector[0];
-<
+    simnfo->setBox( boxVector );
->
+    info->setBox( boxVector );
+  else{
-<
+    if( !the_globals->haveBoxX() ){
->
+    if( !globals->haveBoxX() ){
+      sprintf( painCave.errMsg,
+               "SimSetup error, no periodic BoxX size given.\n" );
+      painCave.isFatal = 1;
+      simError();
-<
+    boxVector[0] = the_globals->getBoxX();
->
+    boxVector[0] = globals->getBoxX();
-<
+    if( !the_globals->haveBoxY() ){
->
+    if( !globals->haveBoxY() ){
+      sprintf( painCave.errMsg,
+               "SimSetup error, no periodic BoxY size given.\n" );
+      painCave.isFatal = 1;
+      simError();
-<
+    boxVector[1] = the_globals->getBoxY();
->
+    boxVector[1] = globals->getBoxY();
-<
+    if( !the_globals->haveBoxZ() ){
->
+    if( !globals->haveBoxZ() ){
+      sprintf( painCave.errMsg,
+               "SimSetup error, no periodic BoxZ size given.\n" );
+      painCave.isFatal = 1;
+      simError();
-<
+    boxVector[2] = the_globals->getBoxZ();
->
+    boxVector[2] = globals->getBoxZ();
-<
+    simnfo->setBox( boxVector );
->
+    info->setBox( boxVector );
-+
-+
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "Box size set up" );
->
+  strcpy( checkPointMsg, "Succesfully gathered all information from Bass\n" );
+  MPIcheckPoint();
+#endif // is_mpi
-+
+}
-–
+  // initialize the arrays
-<
+  the_ff->setSimInfo( simnfo );
->
+void SimSetup::finalInfoCheck( void ){
->
+  int index;
->
+  int usesDipoles;
->
-<
+  makeMolecules();
-<
+  simnfo->identArray = new int[simnfo->n_atoms];
-<
+  for(i=0; i<simnfo->n_atoms; i++){
-<
+    simnfo->identArray[i] = the_atoms[i]->getIdent();
->
+  // check electrostatic parameters
->
->
+  index = 0;
->
+  usesDipoles = 0;
->
+  while( (index < info->n_atoms) && !usesDipoles ){
->
+    usesDipoles = ((info->atoms)[index])->hasDipole();
->
+    index++;
-<
+  if (the_globals->getUseRF() ) {
-<
+    simnfo->useReactionField = 1;
-<
-<
+    if( !the_globals->haveECR() ){
->
+#ifdef IS_MPI
->
+  int myUse = usesDipoles;
->
+  MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
->
+#endif //is_mpi
->
->
->
+  if (globals->getUseRF() ) {
->
+    info->useReactionField = 1;
->
->
+    if( !globals->haveECR() ){
+      sprintf( painCave.errMsg,
+               "SimSetup Warning: using default value of 1/2 the smallest "
+               "box length for the electrostaticCutoffRadius.\n"
+      painCave.isFatal = 0;
+      simError();
+      double smallest;
-<
+      smallest = simnfo->boxLx;
-<
+      if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
-<
+      if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
-<
+      simnfo->ecr = 0.5 * smallest;
->
+      smallest = info->boxL[0];
->
+      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
->
+      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
->
+      info->ecr = 0.5 * smallest;
+    } else {
-<
+      simnfo->ecr        = the_globals->getECR();
->
+      info->ecr        = globals->getECR();
-<
+    if( !the_globals->haveEST() ){
->
+    if( !globals->haveEST() ){
+      sprintf( painCave.errMsg,
+               "SimSetup Warning: using default value of 0.05 * the "
+               "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
+               );
+      painCave.isFatal = 0;
+      simError();
-<
+      simnfo->est = 0.05 * simnfo->ecr;
->
+      info->est = 0.05 * info->ecr;
+    } else {
-<
+      simnfo->est        = the_globals->getEST();
->
+      info->est        = globals->getEST();
-<
+    if(!the_globals->haveDielectric() ){
->
+    if(!globals->haveDielectric() ){
+      sprintf( painCave.errMsg,
+               "SimSetup Error: You are trying to use Reaction Field without"
+               "setting a dielectric constant!\n"
+      painCave.isFatal = 1;
+      simError();
-<
+    simnfo->dielectric = the_globals->getDielectric();
-<
+  } else {
->
+    info->dielectric = globals->getDielectric();
->
+  }
->
+  else {
+    if (usesDipoles) {
-<
+      if( !the_globals->haveECR() ){
->
+      if( !globals->haveECR() ){
+        sprintf( painCave.errMsg,
+                 "SimSetup Warning: using default value of 1/2 the smallest "
+                 "box length for the electrostaticCutoffRadius.\n"
+        painCave.isFatal = 0;
+        simError();
+        double smallest;
-<
+        smallest = simnfo->boxLx;
-<
+        if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
-<
+        if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
-<
+        simnfo->ecr = 0.5 * smallest;
->
+        smallest = info->boxL[0];
->
+        if (info->boxL[1] <= smallest) smallest = info->boxL[1];
->
+        if (info->boxL[2] <= smallest) smallest = info->boxL[2];
->
+        info->ecr = 0.5 * smallest;
+      } else {
-<
+        simnfo->ecr        = the_globals->getECR();
->
+        info->ecr        = globals->getECR();
-<
+      if( !the_globals->haveEST() ){
->
+      if( !globals->haveEST() ){
+        sprintf( painCave.errMsg,
+                 "SimSetup Warning: using default value of 5%% of the "
+                 "electrostaticCutoffRadius for the "
+                 );
+        painCave.isFatal = 0;
+        simError();
-<
+        simnfo->est = 0.05 * simnfo->ecr;
->
+        info->est = 0.05 * info->ecr;
+      } else {
-<
+        simnfo->est        = the_globals->getEST();
->
+        info->est        = globals->getEST();
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "electrostatic parameters check out" );
->
+  strcpy( checkPointMsg, "post processing checks out" );
+  MPIcheckPoint();
+#endif // is_mpi
-<
+ if( the_globals->haveInitialConfig() ){
->
+}
->
->
+void SimSetup::initSystemCoords( void ){
->
->
+ if( globals->haveInitialConfig() ){
+     InitializeFromFile* fileInit;
+#ifdef IS_MPI // is_mpi
+     if( worldRank == 0 ){
+#endif //is_mpi
-<
+   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
->
+   fileInit = new InitializeFromFile( globals->getInitialConfig() );
+#ifdef IS_MPI
+     }else fileInit = new InitializeFromFile( NULL );
+#endif
-<
+   fileInit->read_xyz( simnfo ); // default velocities on
->
+   fileInit->read_xyz( info ); // default velocities on
+   delete fileInit;
+  MPIcheckPoint();
+#endif // is_mpi
-+
+}
-–
-–
-–
-<
->
+void SimSetup::makeOutNames( void ){
->
+#ifdef IS_MPI
+  if( worldRank == 0 ){
+#endif // is_mpi
-<
+    if( the_globals->haveFinalConfig() ){
-<
+      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
->
+    if( globals->haveFinalConfig() ){
->
+      strcpy( info->finalName, globals->getFinalConfig() );
+    else{
-<
+      strcpy( simnfo->finalName, inFileName );
->
+      strcpy( info->finalName, inFileName );
+      char* endTest;
-<
+      int nameLength = strlen( simnfo->finalName );
-<
+      endTest = &(simnfo->finalName[nameLength - 5]);
->
+      int nameLength = strlen( info->finalName );
->
+      endTest = &(info->finalName[nameLength - 5]);
+      if( !strcmp( endTest, ".bass" ) ){
+        strcpy( endTest, ".eor" );
+        strcpy( endTest, ".eor" );
+      else{
-<
+        endTest = &(simnfo->finalName[nameLength - 4]);
->
+        endTest = &(info->finalName[nameLength - 4]);
+        if( !strcmp( endTest, ".bss" ) ){
+          strcpy( endTest, ".eor" );
+          strcpy( endTest, ".eor" );
+        else{
-<
+          strcat( simnfo->finalName, ".eor" );
->
+          strcat( info->finalName, ".eor" );
+    // make the sample and status out names
-<
+    strcpy( simnfo->sampleName, inFileName );
->
+    strcpy( info->sampleName, inFileName );
+    char* endTest;
-<
+    int nameLength = strlen( simnfo->sampleName );
-<
+    endTest = &(simnfo->sampleName[nameLength - 5]);
->
+    int nameLength = strlen( info->sampleName );
->
+    endTest = &(info->sampleName[nameLength - 5]);
+    if( !strcmp( endTest, ".bass" ) ){
+      strcpy( endTest, ".dump" );
+      strcpy( endTest, ".dump" );
+    else{
-<
+      endTest = &(simnfo->sampleName[nameLength - 4]);
->
+      endTest = &(info->sampleName[nameLength - 4]);
+      if( !strcmp( endTest, ".bss" ) ){
+        strcpy( endTest, ".dump" );
+        strcpy( endTest, ".dump" );
+      else{
-<
+        strcat( simnfo->sampleName, ".dump" );
->
+        strcat( info->sampleName, ".dump" );
-<
+    strcpy( simnfo->statusName, inFileName );
-<
+    nameLength = strlen( simnfo->statusName );
-<
+    endTest = &(simnfo->statusName[nameLength - 5]);
->
+    strcpy( info->statusName, inFileName );
->
+    nameLength = strlen( info->statusName );
->
+    endTest = &(info->statusName[nameLength - 5]);
+    if( !strcmp( endTest, ".bass" ) ){
+      strcpy( endTest, ".stat" );
+      strcpy( endTest, ".stat" );
+    else{
-<
+      endTest = &(simnfo->statusName[nameLength - 4]);
->
+      endTest = &(info->statusName[nameLength - 4]);
+      if( !strcmp( endTest, ".bss" ) ){
+        strcpy( endTest, ".stat" );
+        strcpy( endTest, ".stat" );
+      else{
-<
+        strcat( simnfo->statusName, ".stat" );
->
+        strcat( info->statusName, ".stat" );
+#ifdef IS_MPI
+#endif // is_mpi
-+
-+
+}
-+
-+
-+
+void SimSetup::sysObjectsCreation( void ){
-+
-+
+  int i;
-+
-+
+  // create the forceField
-+
-+
+  createFF();
-+
-+
+  // extract componentList
-+
-+
+  compList();
-+
-+
+  // calc the number of atoms, bond, bends, and torsions
-+
-+
+  calcSysValues();
-+
-+
+#ifdef IS_MPI
-+
+  // divide the molecules among the processors
-<
+  // set the status, sample, and themal kick times
->
+  mpiMolDivide();
->
+#endif //is_mpi
-<
+  if( the_globals->haveSampleTime() ){
-<
+    simnfo->sampleTime = the_globals->getSampleTime();
-<
+    simnfo->statusTime = simnfo->sampleTime;
-<
+    simnfo->thermalTime = simnfo->sampleTime;
->
+  // create the atom and SRI arrays. Also initialize Molecule Stamp ID's
->
->
+  makeSysArrays();
->
->
+  // make and initialize the molecules (all but atomic coordinates)
->
->
+  makeMolecules();
->
+  info->identArray = new int[info->n_atoms];
->
+  for(i=0; i<info->n_atoms; i++){
->
+    info->identArray[i] = the_atoms[i]->getIdent();
-<
+  else{
-<
+    simnfo->sampleTime = the_globals->getRunTime();
-<
+    simnfo->statusTime = simnfo->sampleTime;
-<
+    simnfo->thermalTime = simnfo->sampleTime;
->
->
->
->
+}
->
->
->
+void SimSetup::createFF( void ){
->
->
+  switch( ffCase ){
->
->
+  case FF_DUFF:
->
+    the_ff = new DUFF();
->
+    break;
->
->
+  case FF_LJ:
->
+    the_ff = new LJFF();
->
+    break;
->
->
+  default:
->
+    sprintf( painCave.errMsg,
->
+             "SimSetup Error. Unrecognized force field in case statement.\n");
->
+    painCave.isFatal = 1;
->
+    simError();
-<
+  if( the_globals->haveStatusTime() ){
-<
+    simnfo->statusTime = the_globals->getStatusTime();
->
+#ifdef IS_MPI
->
+  strcpy( checkPointMsg, "ForceField creation successful" );
->
+  MPIcheckPoint();
->
+#endif // is_mpi
->
->
+}
->
->
->
+void SimSetup::compList( void ){
->
->
+  int i;
->
->
+  comp_stamps = new MoleculeStamp*[n_components];
->
->
+  // make an array of molecule stamps that match the components used.
->
+  // also extract the used stamps out into a separate linked list
->
->
+  info->nComponents = n_components;
->
+  info->componentsNmol = components_nmol;
->
+  info->compStamps = comp_stamps;
->
+  info->headStamp = new LinkedMolStamp();
->
->
+  char* id;
->
+  LinkedMolStamp* headStamp = info->headStamp;
->
+  LinkedMolStamp* currentStamp = NULL;
->
+  for( i=0; i<n_components; i++ ){
->
->
+    id = the_components[i]->getType();
->
+    comp_stamps[i] = NULL;
->
->
+    // check to make sure the component isn't already in the list
->
->
+    comp_stamps[i] = headStamp->match( id );
->
+    if( comp_stamps[i] == NULL ){
->
->
+      // extract the component from the list;
->
->
+      currentStamp = stamps->extractMolStamp( id );
->
+      if( currentStamp == NULL ){
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: Component \"%s\" was not found in the "
->
+                 "list of declared molecules\n",
->
+                 id );
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      headStamp->add( currentStamp );
->
+      comp_stamps[i] = headStamp->match( id );
->
+    }
-<
+  if( the_globals->haveThermalTime() ){
-<
+    simnfo->thermalTime = the_globals->getThermalTime();
->
+#ifdef IS_MPI
->
+  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
->
+  MPIcheckPoint();
->
+#endif // is_mpi
->
->
->
+}
->
->
+void SimSetup::calcSysValues( void ){
->
+  int i, j, k;
->
->
->
+  tot_atoms = 0;
->
+  tot_bonds = 0;
->
+  tot_bends = 0;
->
+  tot_torsions = 0;
->
+  for( i=0; i<n_components; i++ ){
->
->
+    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
->
+    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
->
+    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
->
+    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
-<
+  // check for the temperature set flag
->
+  tot_SRI = tot_bonds + tot_bends + tot_torsions;
-<
+  if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
->
+  info->n_atoms = tot_atoms;
->
+  info->n_bonds = tot_bonds;
->
+  info->n_bends = tot_bends;
->
+  info->n_torsions = tot_torsions;
->
+  info->n_SRI = tot_SRI;
->
+  info->n_mol = tot_nmol;
->
->
+  info->molMembershipArray = new int[tot_atoms];
->
+}
-<
+  // make the integrator
->
+#ifdef IS_MPI
->
->
+void SimSetup::mpiMolDivide( void ){
-+
+  int i, j, k;
-+
+  int localMol, allMol;
-+
+  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
-+
-+
+  mpiSim = new mpiSimulation( info );
-+
+  globalIndex = mpiSim->divideLabor();
-+
-+
+  // set up the local variables
-+
-+
+  mol2proc = mpiSim->getMolToProcMap();
-+
+  molCompType = mpiSim->getMolComponentType();
-+
-+
+  allMol = 0;
-+
+  localMol = 0;
-+
+  local_atoms = 0;
-+
+  local_bonds = 0;
-+
+  local_bends = 0;
-+
+  local_torsions = 0;
-+
+  globalAtomIndex = 0;
-+
-+
-+
+  for( i=0; i<n_components; i++ ){
-+
-+
+    for( j=0; j<components_nmol[i]; j++ ){
-+
-+
+      if( mol2proc[allMol] == worldRank ){
-+
-+
+        local_atoms +=    comp_stamps[i]->getNAtoms();
-+
+        local_bonds +=    comp_stamps[i]->getNBonds();
-+
+        local_bends +=    comp_stamps[i]->getNBends();
-+
+        local_torsions += comp_stamps[i]->getNTorsions();
-+
+        localMol++;
-+
+      }
-+
+      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
-+
+        info->molMembershipArray[globalAtomIndex] = allMol;
-+
+        globalAtomIndex++;
-+
+      }
-+
-+
+      allMol++;
-+
+    }
-+
+  }
-+
+  local_SRI = local_bonds + local_bends + local_torsions;
-+
-+
+  info->n_atoms = mpiSim->getMyNlocal();
-+
-+
+  if( local_atoms != info->n_atoms ){
-+
+    sprintf( painCave.errMsg,
-+
+             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
-+
+             " localAtom (%d) are not equal.\n",
-+
+             info->n_atoms,
-+
+             local_atoms );
-+
+    painCave.isFatal = 1;
-+
+    simError();
-+
+  }
-+
-+
+  info->n_bonds = local_bonds;
-+
+  info->n_bends = local_bends;
-+
+  info->n_torsions = local_torsions;
-+
+  info->n_SRI = local_SRI;
-+
+  info->n_mol = localMol;
-+
-+
+  strcpy( checkPointMsg, "Passed nlocal consistency check." );
-+
+  MPIcheckPoint();
-+
+}
-+
-+
+#endif // is_mpi
-+
-+
-+
+void SimSetup::makeSysArrays( void ){
-+
+  int i, j, k;
-+
-+
-+
+  // create the atom and short range interaction arrays
-+
-+
+  Atom::createArrays(info->n_atoms);
-+
+  the_atoms = new Atom*[info->n_atoms];
-+
+  the_molecules = new Molecule[info->n_mol];
-+
+  int molIndex;
-+
-+
+  // initialize the molecule's stampID's
-+
-+
+#ifdef IS_MPI
-+
-+
-+
+  molIndex = 0;
-+
+  for(i=0; i<mpiSim->getTotNmol(); i++){
-+
-+
+    if(mol2proc[i] == worldRank ){
-+
+      the_molecules[molIndex].setStampID( molCompType[i] );
-+
+      the_molecules[molIndex].setMyIndex( molIndex );
-+
+      the_molecules[molIndex].setGlobalIndex( i );
-+
+      molIndex++;
-+
+    }
-+
+  }
-+
-+
+#else // is_mpi
-+
-+
+  molIndex = 0;
-+
+  globalAtomIndex = 0;
-+
+  for(i=0; i<n_components; i++){
-+
+    for(j=0; j<components_nmol[i]; j++ ){
-+
+      the_molecules[molIndex].setStampID( i );
-+
+      the_molecules[molIndex].setMyIndex( molIndex );
-+
+      the_molecules[molIndex].setGlobalIndex( molIndex );
-+
+      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
-+
+        info->molMembershipArray[globalAtomIndex] = molIndex;
-+
+        globalAtomIndex++;
-+
+      }
-+
+      molIndex++;
-+
+    }
-+
+  }
-+
-+
-+
+#endif // is_mpi
-+
-+
-+
+  if( info->n_SRI ){
-+
-+
+    Exclude::createArray(info->n_SRI);
-+
+    the_excludes = new Exclude*[info->n_SRI];
-+
+    for( int ex=0; ex<info->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
-+
+    info->globalExcludes = new int;
-+
+    info->n_exclude = info->n_SRI;
-+
+  }
-+
+  else{
-+
-+
+    Exclude::createArray( 1 );
-+
+    the_excludes = new Exclude*;
-+
+    the_excludes[0] = new Exclude(0);
-+
+    the_excludes[0]->setPair( 0,0 );
-+
+    info->globalExcludes = new int;
-+
+    info->globalExcludes[0] = 0;
-+
+    info->n_exclude = 0;
-+
+  }
-+
-+
+  // set the arrays into the SimInfo object
-+
-+
+  info->atoms = the_atoms;
-+
+  info->molecules = the_molecules;
-+
+  info->nGlobalExcludes = 0;
-+
+  info->excludes = the_excludes;
-+
-+
+  the_ff->setSimInfo( info );
-+
-+
+}
-+
-+
+void SimSetup::makeIntegrator( void ){
-+
+  NVT*  myNVT = NULL;
+  NPTi* myNPTi = NULL;
+  NPTf* myNPTf = NULL;
+  switch( ensembleCase ){
+  case NVE_ENS:
-<
+    new NVE( simnfo, the_ff );
->
+    new NVE( info, the_ff );
+    break;
+  case NVT_ENS:
-<
+    myNVT = new NVT( simnfo, the_ff );
-<
+    myNVT->setTargetTemp(the_globals->getTargetTemp());
->
+    myNVT = new NVT( info, the_ff );
->
+    myNVT->setTargetTemp(globals->getTargetTemp());
-<
+    if (the_globals->haveTauThermostat())
-<
+      myNVT->setTauThermostat(the_globals->getTauThermostat());
->
+    if (globals->haveTauThermostat())
->
+      myNVT->setTauThermostat(globals->getTauThermostat());
+    else {
+      sprintf( painCave.errMsg,
+    break;
+  case NPTi_ENS:
-<
+    myNPTi = new NPTi( simnfo, the_ff );
-<
+    myNPTi->setTargetTemp( the_globals->getTargetTemp());
->
+    myNPTi = new NPTi( info, the_ff );
->
+    myNPTi->setTargetTemp( globals->getTargetTemp() );
-<
+    if (the_globals->haveTargetPressure())
-<
+      myNPTi->setTargetPressure(the_globals->getTargetPressure());
->
+    if (globals->haveTargetPressure())
->
+      myNPTi->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+      simError();
-<
+    if( the_globals->haveTauThermostat() )
-<
+      myNPTi->setTauThermostat( the_globals->getTauThermostat() );
->
+    if( globals->haveTauThermostat() )
->
+      myNPTi->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+      simError();
-<
+    if( the_globals->haveTauBarostat() )
-<
+      myNPTi->setTauBarostat( the_globals->getTauBarostat() );
->
+    if( globals->haveTauBarostat() )
->
+      myNPTi->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+    break;
+  case NPTf_ENS:
-<
+    myNPTf = new NPTf( simnfo, the_ff );
-<
+    myNPTf->setTargetTemp( the_globals->getTargetTemp());
->
+    myNPTf = new NPTf( info, the_ff );
->
+    myNPTf->setTargetTemp( globals->getTargetTemp());
-<
+    if (the_globals->haveTargetPressure())
-<
+      myNPTf->setTargetPressure(the_globals->getTargetPressure());
->
+    if (globals->haveTargetPressure())
->
+      myNPTf->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+      simError();
-<
+    if( the_globals->haveTauThermostat() )
-<
+      myNPTf->setTauThermostat( the_globals->getTauThermostat() );
->
+    if( globals->haveTauThermostat() )
->
+      myNPTf->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+      simError();
-<
+    if( the_globals->haveTauBarostat() )
-<
+      myNPTf->setTauBarostat( the_globals->getTauBarostat() );
->
+    if( globals->haveTauBarostat() )
->
+      myNPTf->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+    break;
+  case NPTim_ENS:
-<
+    myNPTim = new NPTim( simnfo, the_ff );
-<
+    myNPTim->setTargetTemp( the_globals->getTargetTemp());
->
+    myNPTim = new NPTim( info, the_ff );
->
+    myNPTim->setTargetTemp( globals->getTargetTemp());
-<
+    if (the_globals->haveTargetPressure())
-<
+      myNPTim->setTargetPressure(the_globals->getTargetPressure());
->
+    if (globals->haveTargetPressure())
->
+      myNPTim->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+      simError();
-<
+    if( the_globals->haveTauThermostat() )
-<
+      myNPTim->setTauThermostat( the_globals->getTauThermostat() );
->
+    if( globals->haveTauThermostat() )
->
+      myNPTim->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+      simError();
-<
+    if( the_globals->haveTauBarostat() )
-<
+      myNPTim->setTauBarostat( the_globals->getTauBarostat() );
->
+    if( globals->haveTauBarostat() )
->
+      myNPTim->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+    break;
+  case NPTfm_ENS:
-<
+    myNPTfm = new NPTfm( simnfo, the_ff );
-<
+    myNPTfm->setTargetTemp( the_globals->getTargetTemp());
->
+    myNPTfm = new NPTfm( info, the_ff );
->
+    myNPTfm->setTargetTemp( globals->getTargetTemp());
-<
+    if (the_globals->haveTargetPressure())
-<
+      myNPTfm->setTargetPressure(the_globals->getTargetPressure());
->
+    if (globals->haveTargetPressure())
->
+      myNPTfm->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+      simError();
-<
+    if( the_globals->haveTauThermostat() )
-<
+      myNPTfm->setTauThermostat( the_globals->getTauThermostat() );
->
+    if( globals->haveTauThermostat() )
->
+      myNPTfm->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+      simError();
-<
+    if( the_globals->haveTauBarostat() )
-<
+      myNPTfm->setTauBarostat( the_globals->getTauBarostat() );
->
+    if( globals->haveTauBarostat() )
->
+      myNPTfm->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+    simError();
-+
+}
-<
+#ifdef IS_MPI
-<
+  mpiSim->mpiRefresh();
-<
+#endif
->
+void SimSetup::initFortran( void ){
-<
+  // initialize the Fortran
-<
-<
-<
+  simnfo->refreshSim();
->
+  info->refreshSim();
-<
+  if( !strcmp( simnfo->mixingRule, "standard") ){
->
+  if( !strcmp( info->mixingRule, "standard") ){
+    the_ff->initForceField( LB_MIXING_RULE );
-<
+  else if( !strcmp( simnfo->mixingRule, "explicit") ){
->
+  else if( !strcmp( info->mixingRule, "explicit") ){
+    the_ff->initForceField( EXPLICIT_MIXING_RULE );
+  else{
+    sprintf( painCave.errMsg,
+             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
-<
+             simnfo->mixingRule );
->
+             info->mixingRule );
+    painCave.isFatal = 1;
+    simError();
+          "Successfully intialized the mixingRule for Fortran." );
+  MPIcheckPoint();
+#endif // is_mpi
-–
+}
-–
-–
+void SimSetup::makeMolecules( void ){
-–
-–
+  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
-–
+  molInit info;
-–
+  DirectionalAtom* dAtom;
-–
+  LinkedAssign* extras;
-–
+  LinkedAssign* current_extra;
-–
+  AtomStamp* currentAtom;
-–
+  BondStamp* currentBond;
-–
+  BendStamp* currentBend;
-–
+  TorsionStamp* currentTorsion;
-–
-–
+  bond_pair* theBonds;
-–
+  bend_set* theBends;
-–
+  torsion_set* theTorsions;
-–
-–
-–
+  //init the forceField paramters
-–
-–
+  the_ff->readParams();
-–
-–
-–
+  // init the atoms
-–
-–
+  double ux, uy, uz, u, uSqr;
-–
-–
+  atomOffset = 0;
-–
+  excludeOffset = 0;
-–
+  for(i=0; i<simnfo->n_mol; i++){
-–
-–
+    stampID = the_molecules[i].getStampID();
-–
-–
+    info.nAtoms    = comp_stamps[stampID]->getNAtoms();
-–
+    info.nBonds    = comp_stamps[stampID]->getNBonds();
-–
+    info.nBends    = comp_stamps[stampID]->getNBends();
-–
+    info.nTorsions = comp_stamps[stampID]->getNTorsions();
-–
+    info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
-–
-–
+    info.myAtoms = &the_atoms[atomOffset];
-–
+    info.myExcludes = &the_excludes[excludeOffset];
-–
+    info.myBonds = new Bond*[info.nBonds];
-–
+    info.myBends = new Bend*[info.nBends];
-–
+    info.myTorsions = new Torsion*[info.nTorsions];
-–
-–
+    theBonds = new bond_pair[info.nBonds];
-–
+    theBends = new bend_set[info.nBends];
-–
+    theTorsions = new torsion_set[info.nTorsions];
-–
-–
+    // make the Atoms
-–
-–
+    for(j=0; j<info.nAtoms; j++){
-–
-–
+      currentAtom = comp_stamps[stampID]->getAtom( j );
-–
+      if( currentAtom->haveOrientation() ){
-–
-–
+        dAtom = new DirectionalAtom(j + atomOffset);
-–
+        simnfo->n_oriented++;
-–
+        info.myAtoms[j] = dAtom;
-–
-–
+        ux = currentAtom->getOrntX();
-–
+        uy = currentAtom->getOrntY();
-–
+        uz = currentAtom->getOrntZ();
-–
-–
+        uSqr = (ux * ux) + (uy * uy) + (uz * uz);
-–
-–
+        u = sqrt( uSqr );
-–
+        ux = ux / u;
-–
+        uy = uy / u;
-–
+        uz = uz / u;
-–
-–
+        dAtom->setSUx( ux );
-–
+        dAtom->setSUy( uy );
-–
+        dAtom->setSUz( uz );
-–
+      }
-–
+      else{
-–
+        info.myAtoms[j] = new GeneralAtom(j + atomOffset);
-–
+      }
-–
+      info.myAtoms[j]->setType( currentAtom->getType() );
-–
-–
+#ifdef IS_MPI
-–
-–
+      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
-–
-–
+#endif // is_mpi
-–
+    }
-–
-–
+    // make the bonds
-–
+    for(j=0; j<info.nBonds; j++){
-–
-–
+      currentBond = comp_stamps[stampID]->getBond( j );
-–
+      theBonds[j].a = currentBond->getA() + atomOffset;
-–
+      theBonds[j].b = currentBond->getB() + atomOffset;
-–
-–
+      exI = theBonds[j].a;
-–
+      exJ = theBonds[j].b;
-–
-–
+      // exclude_I must always be the smaller of the pair
-–
+      if( exI > exJ ){
-–
+        tempEx = exI;
-–
+        exI = exJ;
-–
+        exJ = tempEx;
-–
+      }
-–
+#ifdef IS_MPI
-–
+      tempEx = exI;
-–
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
-–
+      tempEx = exJ;
-–
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
-–
-–
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
-–
+#else  // isn't MPI
-–
-–
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
-–
+#endif  //is_mpi
-–
+    }
-–
+    excludeOffset += info.nBonds;
-–
-–
+    //make the bends
-–
+    for(j=0; j<info.nBends; j++){
-–
-–
+      currentBend = comp_stamps[stampID]->getBend( j );
-–
+      theBends[j].a = currentBend->getA() + atomOffset;
-–
+      theBends[j].b = currentBend->getB() + atomOffset;
-–
+      theBends[j].c = currentBend->getC() + atomOffset;
-–
-–
+      if( currentBend->haveExtras() ){
-–
-–
+        extras = currentBend->getExtras();
-–
+        current_extra = extras;
-–
-–
+        while( current_extra != NULL ){
-–
+          if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
-–
-–
+            switch( current_extra->getType() ){
-–
-–
+            case 0:
-–
+              theBends[j].ghost =
-–
+                current_extra->getInt() + atomOffset;
-–
+              theBends[j].isGhost = 1;
-–
+              break;
-–
-–
+            case 1:
-–
+              theBends[j].ghost =
-–
+                (int)current_extra->getDouble() + atomOffset;
-–
+              theBends[j].isGhost = 1;
-–
+              break;
-–
-–
+            default:
-–
+              sprintf( painCave.errMsg,
-–
+                       "SimSetup Error: ghostVectorSource was neither a "
-–
+                       "double nor an int.\n"
-–
+                       "-->Bend[%d] in %s\n",
-–
+                       j, comp_stamps[stampID]->getID() );
-–
+              painCave.isFatal = 1;
-–
+              simError();
-–
+            }
-–
+          }
-–
-–
+          else{
-–
-–
+            sprintf( painCave.errMsg,
-–
+                     "SimSetup Error: unhandled bend assignment:\n"
-–
+                     "    -->%s in Bend[%d] in %s\n",
-–
+                     current_extra->getlhs(),
-–
+                     j, comp_stamps[stampID]->getID() );
-–
+            painCave.isFatal = 1;
-–
+            simError();
-–
+          }
-–
-–
+          current_extra = current_extra->getNext();
-–
+        }
-–
+      }
-–
-–
+      if( !theBends[j].isGhost ){
-–
-–
+        exI = theBends[j].a;
-–
+        exJ = theBends[j].c;
-–
+      }
-–
+      else{
-–
-–
+        exI = theBends[j].a;
-–
+        exJ = theBends[j].b;
-–
+      }
-–
-–
+      // exclude_I must always be the smaller of the pair
-–
+      if( exI > exJ ){
-–
+        tempEx = exI;
-–
+        exI = exJ;
-–
+        exJ = tempEx;
-–
+      }
-–
+#ifdef IS_MPI
-–
+      tempEx = exI;
-–
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
-–
+      tempEx = exJ;
-–
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
-–
-–
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
-–
+#else  // isn't MPI
-–
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
-–
+#endif  //is_mpi
-–
+    }
-–
+    excludeOffset += info.nBends;
-–
-–
+    for(j=0; j<info.nTorsions; j++){
-–
-–
+      currentTorsion = comp_stamps[stampID]->getTorsion( j );
-–
+      theTorsions[j].a = currentTorsion->getA() + atomOffset;
-–
+      theTorsions[j].b = currentTorsion->getB() + atomOffset;
-–
+      theTorsions[j].c = currentTorsion->getC() + atomOffset;
-–
+      theTorsions[j].d = currentTorsion->getD() + atomOffset;
-–
-–
+      exI = theTorsions[j].a;
-–
+      exJ = theTorsions[j].d;
-–
-–
+      // exclude_I must always be the smaller of the pair
-–
+      if( exI > exJ ){
-–
+        tempEx = exI;
-–
+        exI = exJ;
-–
+        exJ = tempEx;
-–
+      }
-–
+#ifdef IS_MPI
-–
+      tempEx = exI;
-–
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
-–
+      tempEx = exJ;
-–
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
-–
-–
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
-–
+#else  // isn't MPI
-–
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
-–
+#endif  //is_mpi
-–
+    }
-–
+    excludeOffset += info.nTorsions;
-–
-–
-–
+    // send the arrays off to the forceField for init.
-–
-–
+    the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
-–
+    the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
-–
+    the_ff->initializeBends( info.nBends, info.myBends, theBends );
-–
+    the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
-–
-–
-–
+    the_molecules[i].initialize( info );
-–
-–
-–
+    atomOffset += info.nAtoms;
-–
+    delete[] theBonds;
-–
+    delete[] theBends;
-–
+    delete[] theTorsions;
-–
+  }
-–
-–
+#ifdef IS_MPI
-–
+  sprintf( checkPointMsg, "all molecules initialized succesfully" );
-–
+  MPIcheckPoint();
-–
+#endif // is_mpi
-–
-–
+  // clean up the forcefield
-–
+  the_ff->calcRcut();
-–
+  the_ff->cleanMe();
-–
-–
-–
+void SimSetup::initFromBass( void ){
-–
-–
+  int i, j, k;
-–
+  int n_cells;
-–
+  double cellx, celly, cellz;
-–
+  double temp1, temp2, temp3;
-–
+  int n_per_extra;
-–
+  int n_extra;
-–
+  int have_extra, done;
-–
-–
+  temp1 = (double)tot_nmol / 4.0;
-–
+  temp2 = pow( temp1, ( 1.0 / 3.0 ) );
-–
+  temp3 = ceil( temp2 );
-–
-–
+  have_extra =0;
-–
+  if( temp2 < temp3 ){ // we have a non-complete lattice
-–
+    have_extra =1;
-–
-–
+    n_cells = (int)temp3 - 1;
-–
+    cellx = simnfo->boxLx / temp3;
-–
+    celly = simnfo->boxLy / temp3;
-–
+    cellz = simnfo->boxLz / temp3;
-–
+    n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
-–
+    temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
-–
+    n_per_extra = (int)ceil( temp1 );
-–
-–
+    if( n_per_extra > 4){
-–
+      sprintf( painCave.errMsg,
-–
+               "SimSetup error. There has been an error in constructing"
-–
+               " the non-complete lattice.\n" );
-–
+      painCave.isFatal = 1;
-–
+      simError();
-–
+    }
-–
+  }
-–
+  else{
-–
+    n_cells = (int)temp3;
-–
+    cellx = simnfo->boxLx / temp3;
-–
+    celly = simnfo->boxLy / temp3;
-–
+    cellz = simnfo->boxLz / temp3;
-–
+  }
-–
-–
+  current_mol = 0;
-–
+  current_comp_mol = 0;
-–
+  current_comp = 0;
-–
+  current_atom_ndx = 0;
-–
-–
+  for( i=0; i < n_cells ; i++ ){
-–
+    for( j=0; j < n_cells; j++ ){
-–
+      for( k=0; k < n_cells; k++ ){
-–
-–
+        makeElement( i * cellx,
-–
+                     j * celly,
-–
+                     k * cellz );
-–
-–
+        makeElement( i * cellx + 0.5 * cellx,
-–
+                     j * celly + 0.5 * celly,
-–
+                     k * cellz );
-–
-–
+        makeElement( i * cellx,
-–
+                     j * celly + 0.5 * celly,
-–
+                     k * cellz + 0.5 * cellz );
-–
-–
+        makeElement( i * cellx + 0.5 * cellx,
-–
+                     j * celly,
-–
+                     k * cellz + 0.5 * cellz );
-–
+      }
-–
+    }
-–
+  }
-–
-–
+  if( have_extra ){
-–
+    done = 0;
-–
-–
+    int start_ndx;
-–
+    for( i=0; i < (n_cells+1) && !done; i++ ){
-–
+      for( j=0; j < (n_cells+1) && !done; j++ ){
-–
-–
+        if( i < n_cells ){
-–
-–
+          if( j < n_cells ){
-–
+            start_ndx = n_cells;
-–
+          }
-–
+          else start_ndx = 0;
-–
+        }
-–
+        else start_ndx = 0;
-–
-–
+        for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
-–
-–
+          makeElement( i * cellx,
-–
+                       j * celly,
-–
+                       k * cellz );
-–
+          done = ( current_mol >= tot_nmol );
-–
-–
+          if( !done && n_per_extra > 1 ){
-–
+            makeElement( i * cellx + 0.5 * cellx,
-–
+                         j * celly + 0.5 * celly,
-–
+                         k * cellz );
-–
+            done = ( current_mol >= tot_nmol );
-–
+          }
-–
-–
+          if( !done && n_per_extra > 2){
-–
+            makeElement( i * cellx,
-–
+                         j * celly + 0.5 * celly,
-–
+                         k * cellz + 0.5 * cellz );
-–
+            done = ( current_mol >= tot_nmol );
-–
+          }
-–
-–
+          if( !done && n_per_extra > 3){
-–
+            makeElement( i * cellx + 0.5 * cellx,
-–
+                         j * celly,
-–
+                         k * cellz + 0.5 * cellz );
-–
+            done = ( current_mol >= tot_nmol );
-–
+          }
-–
+        }
-–
+      }
-–
+    }
-–
+  }
-–
-–
-–
+  for( i=0; i<simnfo->n_atoms; i++ ){
-–
+    simnfo->atoms[i]->set_vx( 0.0 );
-–
+    simnfo->atoms[i]->set_vy( 0.0 );
-–
+    simnfo->atoms[i]->set_vz( 0.0 );
-–
+  }
-–
+}
-–
-–
+void SimSetup::makeElement( double x, double y, double z ){
-–
-–
+  int k;
-–
+  AtomStamp* current_atom;
-–
+  DirectionalAtom* dAtom;
-–
+  double rotMat[3][3];
-–
-–
+  for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
-–
-–
+    current_atom = comp_stamps[current_comp]->getAtom( k );
-–
+    if( !current_atom->havePosition() ){
-–
+      sprintf( painCave.errMsg,
-–
+               "SimSetup:initFromBass error.\n"
-–
+               "\tComponent %s, atom %s does not have a position specified.\n"
-–
+               "\tThe initialization routine is unable to give a start"
-–
+               " position.\n",
-–
+               comp_stamps[current_comp]->getID(),
-–
+               current_atom->getType() );
-–
+      painCave.isFatal = 1;
-–
+      simError();
-–
+    }
-–
-–
+    the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
-–
+    the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
-–
+    the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
-–
-–
+    if( the_atoms[current_atom_ndx]->isDirectional() ){
-–
-–
+      dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
-–
-–
+      rotMat[0][0] = 1.0;
-–
+      rotMat[0][1] = 0.0;
-–
+      rotMat[0][2] = 0.0;
-–
-–
+      rotMat[1][0] = 0.0;
-–
+      rotMat[1][1] = 1.0;
-–
+      rotMat[1][2] = 0.0;
-–
-–
+      rotMat[2][0] = 0.0;
-–
+      rotMat[2][1] = 0.0;
-–
+      rotMat[2][2] = 1.0;
-–
-–
+      dAtom->setA( rotMat );
-–
+    }
-–
-–
+    current_atom_ndx++;
-–
+  }
-–
-–
+  current_mol++;
-–
+  current_comp_mol++;
-–
-–
+  if( current_comp_mol >= components_nmol[current_comp] ){
-–
-–
+    current_comp_mol = 0;
-–
+    current_comp++;
-–
+  }
-–
+}

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 605 by gezelter, Tue Jul 15 03:27:24 2003 UTC vs. Revision 621 by gezelter, Wed Jul 16 02:11:02 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 605 by gezelter, Tue Jul 15 03:27:24 2003 UTC vs.
Revision 621 by gezelter, Wed Jul 16 02:11:02 2003 UTC