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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 582 by mmeineke, Wed Jul 9 15:33:46 2003 UTC vs.
Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC

-+
+#include <algorithm>
+#include <cstdlib>
+#include <iostream>
+#include <cmath>
-+
+#include <string>
+#include "SimSetup.hpp"
+#include "parse_me.h"
+// some defines for ensemble and Forcefield  cases
-<
+#define NVE_ENS  0
-<
+#define NVT_ENS  1
-<
+#define NPTi_ENS 2
-<
+#define NPTf_ENS 3
->
+#define NVE_ENS        0
->
+#define NVT_ENS        1
->
+#define NPTi_ENS       2
->
+#define NPTf_ENS       3
->
+#define NPTim_ENS      4
->
+#define NPTfm_ENS      5
-–
+#define FF_DUFF 0
+#define FF_LJ   1
-+
+#define FF_EAM  2
-+
+using namespace std;
+SimSetup::SimSetup(){
-+
-+
+  isInfoArray = 0;
-+
+  nInfo = 1;
-+
+  stamps = new MakeStamps();
+  globals = new Globals();
-+
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "SimSetup creation successful" );
+  MPIcheckPoint();
+  delete globals;
-+
+void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
-+
+    info = the_info;
-+
+    nInfo = theNinfo;
-+
+    isInfoArray = 1;
-+
+}
-+
-+
+void SimSetup::parseFile( char* fileName ){
+#ifdef IS_MPI
+#endif // is_mpi
-<
+void SimSetup::createSim( void ){
->
+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
-+
-+
+  if( !isInfoArray ) 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 k,l;
-+
+  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;
-+
-+
+  for(k=0; k<nInfo; k++){
-+
-+
+    the_ff->setSimInfo( &(info[k]) );
-+
-+
+    atomOffset = 0;
-+
+    excludeOffset = 0;
-+
+    for(i=0; i<info[k].n_mol; i++){
-+
-+
+      stampID = info[k].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 = &(info[k].atoms[atomOffset]);
-+
+      molInfo.myExcludes = &(info[k].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[k].getConfiguration() );
-+
+          info[k].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),
-+
+                                                info[k].getConfiguration() );
-+
+        }
-+
+        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 = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-+
+        tempEx = exJ;
-+
+        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-+
-+
+        info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
-+
+#else  // isn't MPI
-+
-+
+        info[k].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 = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-+
+        tempEx = exJ;
-+
+        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-+
-+
+        info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
-+
+#else  // isn't MPI
-+
+        info[k].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 = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-+
+        tempEx = exJ;
-+
+        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-+
-+
+        info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
-+
+#else  // isn't MPI
-+
+        info[k].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 );
-+
-+
-+
+      info[k].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;
-+
-+
+  double vel[3];
-+
+  vel[0] = 0.0;
-+
+  vel[1] = 0.0;
-+
+  vel[2] = 0.0;
-+
-+
+  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[0].boxL[0] / temp3;
-+
+    celly = info[0].boxL[1] / temp3;
-+
+    cellz = info[0].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[0].boxL[0] / temp3;
-+
+    celly = info[0].boxL[1] / temp3;
-+
+    cellz = info[0].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[0].n_atoms; i++ ){
-+
+    info[0].atoms[i]->setVel( vel );
-+
+  }
-+
+}
-+
-+
+void SimSetup::makeElement( double x, double y, double z ){
-+
-+
+  int k;
-+
+  AtomStamp* current_atom;
-+
+  DirectionalAtom* dAtom;
-+
+  double rotMat[3][3];
-+
+  double pos[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();
-+
+    }
-+
-+
+    pos[0] = x + current_atom->getPosX();
-+
+    pos[1] = y + current_atom->getPosY();
-+
+    pos[2] = z + current_atom->getPosZ();
-+
-+
+    info[0].atoms[current_atom_ndx]->setPos( pos );
-+
-+
+    if( info[0].atoms[current_atom_ndx]->isDirectional() ){
-+
-+
+      dAtom = (DirectionalAtom *)info[0].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;
-–
+  // set the easy ones first
-–
+  simnfo->target_temp = the_globals->getTargetTemp();
-–
+  simnfo->dt = the_globals->getDt();
-–
+  simnfo->run_time = the_globals->getRunTime();
-<
+  // get the ones we know are there, yet still may need some work.
-<
+  n_components = the_globals->getNComponents();
-<
+  strcpy( force_field, the_globals->getForceField() );
->
+  for( i=0; i<nInfo; i++){
->
+    info[i].target_temp = globals->getTargetTemp();
->
+    info[i].dt = globals->getDt();
->
+    info[i].run_time = globals->getRunTime();
->
+  }
->
+  n_components = globals->getNComponents();
-+
-+
+  // 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 if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
+  else{
+    sprintf( painCave.errMsg,
+             "SimSetup Error. Unrecognized force field -> %s\n",
+    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") )
+    ensembleCase = NPTi_ENS;
+  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
-+
+  else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
-+
+  else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
+  else{
+    sprintf( painCave.errMsg,
+             "SimSetup Warning. Unrecognized Ensemble -> %s, "
+    strcpy( ensemble, "NVE" );
+    ensembleCase = NVE_ENS;
-<
+  strcpy( simnfo->ensemble, ensemble );
->
->
+  for(i=0; i<nInfo; i++){
->
->
+    strcpy( info[i].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();
->
+    strcpy( info[i].mixingRule, globals->getMixingRule() );
->
+    info[i].usePBC = globals->getPBC();
-<
-<
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "ForceField creation successful" );
-<
+  MPIcheckPoint();
-<
+#endif // is_mpi
-<
->
+  // 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++ ){
->
+  for(i=0; i<nInfo; i++){
-<
+    id = the_components[i]->getType();
-<
+    comp_stamps[i] = NULL;
->
+    if( globals->haveSampleTime() ){
->
+      info[i].sampleTime = globals->getSampleTime();
->
+      info[i].statusTime = info[i].sampleTime;
->
+      info[i].thermalTime = info[i].sampleTime;
->
+    }
->
+    else{
->
+      info[i].sampleTime = globals->getRunTime();
->
+      info[i].statusTime = info[i].sampleTime;
->
+      info[i].thermalTime = info[i].sampleTime;
->
+    }
-<
+    // check to make sure the component isn't already in the list
->
+    if( globals->haveStatusTime() ){
->
+      info[i].statusTime = globals->getStatusTime();
->
+    }
->
->
+    if( globals->haveThermalTime() ){
->
+      info[i].thermalTime = globals->getThermalTime();
->
+    }
-<
+    comp_stamps[i] = headStamp->match( id );
-<
+    if( comp_stamps[i] == NULL ){
->
+    // check for the temperature set flag
->
->
+    if( globals->haveTempSet() ) info[i].setTemp = globals->getTempSet();
->
->
+    // get some of the tricky things that may still be in the globals
->
->
+    double boxVector[3];
->
+    if( globals->haveBox() ){
->
+      boxVector[0] = globals->getBox();
->
+      boxVector[1] = globals->getBox();
->
+      boxVector[2] = globals->getBox();
-<
+      // extract the component from the list;
->
+      info[i].setBox( boxVector );
->
+    }
->
+    else if( globals->haveDensity() ){
-<
+      currentStamp = the_stamps->extractMolStamp( id );
-<
+      if( currentStamp == NULL ){
->
+      double vol;
->
+      vol = (double)tot_nmol / globals->getDensity();
->
+      boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
->
+      boxVector[1] = boxVector[0];
->
+      boxVector[2] = boxVector[0];
->
->
+      info[i].setBox( boxVector );
->
+  }
->
+    else{
->
+      if( !globals->haveBoxX() ){
+        sprintf( painCave.errMsg,
-<
+                 "SimSetup error: Component \"%s\" was not found in the "
-<
+                 "list of declared molecules\n",
-<
+                 id );
->
+                 "SimSetup error, no periodic BoxX size given.\n" );
+        painCave.isFatal = 1;
+        simError();
-+
+      boxVector[0] = globals->getBoxX();
-<
+      headStamp->add( currentStamp );
-<
+      comp_stamps[i] = headStamp->match( id );
->
+      if( !globals->haveBoxY() ){
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error, no periodic BoxY size given.\n" );
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
+      boxVector[1] = globals->getBoxY();
->
->
+      if( !globals->haveBoxZ() ){
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error, no periodic BoxZ size given.\n" );
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
+      boxVector[2] = globals->getBoxZ();
->
->
+      info[i].setBox( boxVector );
-–
+  }
-+
+  }
-+
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
->
+  strcpy( checkPointMsg, "Succesfully gathered all information from Bass\n" );
+  MPIcheckPoint();
+#endif // is_mpi
-–
-+
+}
-<
+  // caclulate the number of atoms, bonds, bends and torsions
->
+void SimSetup::finalInfoCheck( void ){
->
+  int index;
->
+  int usesDipoles;
->
+  int i;
-<
+  tot_atoms = 0;
-<
+  tot_bonds = 0;
-<
+  tot_bends = 0;
-<
+  tot_torsions = 0;
-<
+  for( i=0; i<n_components; i++ ){
->
+  for(i=0; i<nInfo; i++){
->
+    // check electrostatic parameters
-<
+    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();
-<
+  }
-<
-<
+  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;
->
+    index = 0;
->
+    usesDipoles = 0;
->
+    while( (index < info[i].n_atoms) && !usesDipoles ){
->
+      usesDipoles = (info[i].atoms[index])->hasDipole();
->
+      index++;
->
+    }
->
->
+#ifdef IS_MPI
->
+    int myUse = usesDipoles;
->
+    MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
->
+#endif //is_mpi
->
->
+    double theEcr, theEst;
-<
+  simnfo->molMembershipArray = new int[tot_atoms];
->
+    if (globals->getUseRF() ) {
->
+      info[i].useReactionField = 1;
->
->
+      if( !globals->haveECR() ){
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup Warning: using default value of 1/2 the smallest "
->
+                 "box length for the electrostaticCutoffRadius.\n"
->
+                 "I hope you have a very fast processor!\n");
->
+        painCave.isFatal = 0;
->
+        simError();
->
+        double smallest;
->
+        smallest = info[i].boxL[0];
->
+        if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
->
+        if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
->
+        theEcr = 0.5 * smallest;
->
+      } else {
->
+        theEcr = globals->getECR();
->
+      }
->
->
+      if( !globals->haveEST() ){
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup Warning: using default value of 0.05 * the "
->
+                 "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
->
+                 );
->
+        painCave.isFatal = 0;
->
+        simError();
->
+        theEst = 0.05 * theEcr;
->
+      } else {
->
+        theEst= globals->getEST();
->
+      }
->
->
+      info[i].setEcr( theEcr, theEst );
->
->
+      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();
->
+      }
->
+      info[i].dielectric = globals->getDielectric();
->
+    }
->
+    else {
->
+      if (usesDipoles) {
->
->
+        if( !globals->haveECR() ){
->
+          sprintf( painCave.errMsg,
->
+                   "SimSetup Warning: using default value of 1/2 the smallest "
->
+                   "box length for the electrostaticCutoffRadius.\n"
->
+                   "I hope you have a very fast processor!\n");
->
+          painCave.isFatal = 0;
->
+          simError();
->
+          double smallest;
->
+          smallest = info[i].boxL[0];
->
+          if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
->
+          if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
->
+          theEcr = 0.5 * smallest;
->
+        } else {
->
+          theEcr = globals->getECR();
->
+        }
->
->
+        if( !globals->haveEST() ){
->
+          sprintf( painCave.errMsg,
->
+                   "SimSetup Warning: using default value of 0.05 * the "
->
+                   "electrostaticCutoffRadius for the "
->
+                   "electrostaticSkinThickness\n"
->
+                   );
->
+          painCave.isFatal = 0;
->
+          simError();
->
+          theEst = 0.05 * theEcr;
->
+        } else {
->
+          theEst= globals->getEST();
->
+        }
->
->
+        info[i].setEcr( theEcr, theEst );
->
+      }
->
+    }
->
+  }
+#ifdef IS_MPI
-+
+  strcpy( checkPointMsg, "post processing checks out" );
-+
+  MPIcheckPoint();
-+
+#endif // is_mpi
-<
+  // divide the molecules among processors here.
-<
-<
+  mpiSim = new mpiSimulation( simnfo );
-<
-<
+  globalIndex = mpiSim->divideLabor();
->
+}
-<
+  // set up the local variables
->
+void SimSetup::initSystemCoords( void ){
->
+  int i;
-<
+  int localMol, allMol;
-<
+  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
->
+  std::cerr << "Setting atom Coords\n";
-<
+  int* mol2proc = mpiSim->getMolToProcMap();
-<
+  int* molCompType = mpiSim->getMolComponentType();
->
+  (info[0].getConfiguration())->createArrays( info[0].n_atoms );
-<
+  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++) {
-<
+        simnfo->molMembershipArray[globalAtomIndex] = allMol;
-<
+        globalAtomIndex++;
-<
+      }
-<
-<
+      allMol++;
-<
+    }
-<
+  }
-<
+  local_SRI = local_bonds + local_bends + local_torsions;
->
+  for(i=0; i<info[0].n_atoms; i++) info[0].atoms[i]->setCoords();
-<
+  simnfo->n_atoms = mpiSim->getMyNlocal();
-<
-<
+  if( local_atoms != simnfo->n_atoms ){
->
+  if( globals->haveInitialConfig() ){
->
->
+    InitializeFromFile* fileInit;
->
+#ifdef IS_MPI // is_mpi
->
+    if( worldRank == 0 ){
->
+#endif //is_mpi
->
+      fileInit = new InitializeFromFile( globals->getInitialConfig() );
->
+#ifdef IS_MPI
->
+    }else fileInit = new InitializeFromFile( NULL );
->
+#endif
->
+    fileInit->readInit( info ); // default velocities on
->
->
+    delete fileInit;
->
+  }
->
+  else{
->
->
+#ifdef IS_MPI
->
->
+    // no init from bass
->
+    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;
->
+             "Cannot intialize a parallel simulation without an initial configuration file.\n" );
->
+    painCave.isFatal;
+    simError();
-+
-+
+#else
-+
-+
+    initFromBass();
-+
-+
-+
+#endif
-–
-–
+  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();
-<
->
+#ifdef IS_MPI
->
+  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
->
+  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
->
+void SimSetup::makeOutNames( void ){
-+
+  int k;
-–
+  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++;
-<
+    }
-<
+  }
-<
->
+  for(k=0; k<nInfo; k++){
-+
+#ifdef IS_MPI
-+
+    if( worldRank == 0 ){
+#endif // is_mpi
-+
-+
+      if( globals->haveFinalConfig() ){
-+
+        strcpy( info[k].finalName, globals->getFinalConfig() );
-+
+      }
-+
+      else{
-+
+        strcpy( info[k].finalName, inFileName );
-+
+        char* endTest;
-+
+        int nameLength = strlen( info[k].finalName );
-+
+        endTest = &(info[k].finalName[nameLength - 5]);
-+
+        if( !strcmp( endTest, ".bass" ) ){
-+
+          strcpy( endTest, ".eor" );
-+
+        }
-+
+        else if( !strcmp( endTest, ".BASS" ) ){
-+
+          strcpy( endTest, ".eor" );
-+
+        }
-+
+        else{
-+
+          endTest = &(info[k].finalName[nameLength - 4]);
-+
+          if( !strcmp( endTest, ".bss" ) ){
-+
+            strcpy( endTest, ".eor" );
-+
+          }
-+
+          else if( !strcmp( endTest, ".mdl" ) ){
-+
+            strcpy( endTest, ".eor" );
-+
+          }
-+
+          else{
-+
+            strcat( info[k].finalName, ".eor" );
-+
+          }
-+
+        }
-+
+      }
-+
-+
+      // make the sample and status out names
-+
-+
+      strcpy( info[k].sampleName, inFileName );
-+
+      char* endTest;
-+
+      int nameLength = strlen( info[k].sampleName );
-+
+      endTest = &(info[k].sampleName[nameLength - 5]);
-+
+      if( !strcmp( endTest, ".bass" ) ){
-+
+        strcpy( endTest, ".dump" );
-+
+      }
-+
+      else if( !strcmp( endTest, ".BASS" ) ){
-+
+        strcpy( endTest, ".dump" );
-+
+      }
-+
+      else{
-+
+        endTest = &(info[k].sampleName[nameLength - 4]);
-+
+        if( !strcmp( endTest, ".bss" ) ){
-+
+          strcpy( endTest, ".dump" );
-+
+        }
-+
+        else if( !strcmp( endTest, ".mdl" ) ){
-+
+          strcpy( endTest, ".dump" );
-+
+        }
-+
+        else{
-+
+          strcat( info[k].sampleName, ".dump" );
-+
+        }
-+
+      }
-+
-+
+      strcpy( info[k].statusName, inFileName );
-+
+      nameLength = strlen( info[k].statusName );
-+
+      endTest = &(info[k].statusName[nameLength - 5]);
-+
+      if( !strcmp( endTest, ".bass" ) ){
-+
+        strcpy( endTest, ".stat" );
-+
+      }
-+
+      else if( !strcmp( endTest, ".BASS" ) ){
-+
+        strcpy( endTest, ".stat" );
-+
+      }
-+
+      else{
-+
+        endTest = &(info[k].statusName[nameLength - 4]);
-+
+        if( !strcmp( endTest, ".bss" ) ){
-+
+          strcpy( endTest, ".stat" );
-+
+        }
-+
+        else if( !strcmp( endTest, ".mdl" ) ){
-+
+          strcpy( endTest, ".stat" );
-+
+        }
-+
+        else{
-+
+          strcat( info[k].statusName, ".stat" );
-+
+        }
-+
+      }
-+
-+
+#ifdef IS_MPI
-+
+    }
-+
+#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;
-<
+  }
->
+void SimSetup::sysObjectsCreation( void ){
->
->
+  int i,k;
->
->
+  // create the forceField
->
->
+  createFF();
-<
+  // set the arrays into the SimInfo object
->
+  // extract componentList
-<
+  simnfo->atoms = the_atoms;
-<
+  simnfo->molecules = the_molecules;
-<
+  simnfo->nGlobalExcludes = 0;
-<
+  simnfo->excludes = the_excludes;
->
+  compList();
-+
+  // calc the number of atoms, bond, bends, and torsions
-<
+  // get some of the tricky things that may still be in the globals
->
+  calcSysValues();
-<
+  double boxVector[3];
-<
+  if( the_globals->haveBox() ){
-<
+    boxVector[0] = the_globals->getBox();
-<
+    boxVector[1] = the_globals->getBox();
-<
+    boxVector[2] = the_globals->getBox();
-<
-<
+    simnfo->setBox( boxVector );
->
+#ifdef IS_MPI
->
+  // divide the molecules among the processors
->
->
+  mpiMolDivide();
->
+#endif //is_mpi
->
->
+  // create the atom and SRI arrays. Also initialize Molecule Stamp ID's
->
->
+  makeSysArrays();
->
->
+  // make and initialize the molecules (all but atomic coordinates)
->
->
+  makeMolecules();
->
->
+  for(k=0; k<nInfo; k++){
->
+    info[k].identArray = new int[info[k].n_atoms];
->
+    for(i=0; i<info[k].n_atoms; i++){
->
+      info[k].identArray[i] = info[k].atoms[i]->getIdent();
->
+    }
-<
+  else if( the_globals->haveDensity() ){
->
+}
-–
+    double vol;
-–
+    vol = (double)tot_nmol / the_globals->getDensity();
-–
+     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
-–
+     boxVector[1] = boxVector[0];
-–
+     boxVector[2] = boxVector[0];
-<
+    simnfo->setBox( boxVector );
-<
+  }
-<
+  else{
-<
+    if( !the_globals->haveBoxX() ){
-<
+      sprintf( painCave.errMsg,
-<
+               "SimSetup error, no periodic BoxX size given.\n" );
-<
+      painCave.isFatal = 1;
-<
+      simError();
-<
+    }
-<
+    boxVector[0] = the_globals->getBoxX();
->
+void SimSetup::createFF( void ){
-<
+    if( !the_globals->haveBoxY() ){
-<
+      sprintf( painCave.errMsg,
-<
+               "SimSetup error, no periodic BoxY size given.\n" );
-<
+      painCave.isFatal = 1;
-<
+      simError();
-<
+    }
-<
+    boxVector[1] = the_globals->getBoxY();
->
+  switch( ffCase ){
-<
+    if( !the_globals->haveBoxZ() ){
-<
+      sprintf( painCave.errMsg,
-<
+               "SimSetup error, no periodic BoxZ size given.\n" );
-<
+      painCave.isFatal = 1;
-<
+      simError();
-<
+    }
-<
+    boxVector[2] = the_globals->getBoxZ();
->
+  case FF_DUFF:
->
+    the_ff = new DUFF();
->
+    break;
-<
+    simnfo->setBox( boxVector );
->
+  case FF_LJ:
->
+    the_ff = new LJFF();
->
+    break;
->
->
+  case FF_EAM:
->
+    the_ff = new EAM_FF();
->
+    break;
->
->
+  default:
->
+    sprintf( painCave.errMsg,
->
+             "SimSetup Error. Unrecognized force field in case statement.\n");
->
+    painCave.isFatal = 1;
->
+    simError();
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "Box size set up" );
->
+  strcpy( checkPointMsg, "ForceField creation successful" );
+  MPIcheckPoint();
+#endif // is_mpi
-+
+}
-–
+  // initialize the arrays
-<
+  the_ff->setSimInfo( simnfo );
->
+void SimSetup::compList( void ){
-<
+  makeMolecules();
-<
+  simnfo->identArray = new int[simnfo->n_atoms];
-<
+  for(i=0; i<simnfo->n_atoms; i++){
-<
+    simnfo->identArray[i] = the_atoms[i]->getIdent();
-<
+  }
->
+  int i;
->
+  char* id;
->
+  LinkedMolStamp* headStamp = new LinkedMolStamp();
->
+  LinkedMolStamp* currentStamp = NULL;
->
+  comp_stamps = new MoleculeStamp*[n_components];
-<
+  if (the_globals->getUseRF() ) {
-<
+    simnfo->useReactionField = 1;
->
+  // make an array of molecule stamps that match the components used.
->
+  // also extract the used stamps out into a separate linked list
-<
+    if( !the_globals->haveECR() ){
-<
+      sprintf( painCave.errMsg,
-<
+               "SimSetup Warning: using default value of 1/2 the smallest "
-<
+               "box length for the electrostaticCutoffRadius.\n"
-<
+               "I hope you have a very fast processor!\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;
-<
+    } else {
-<
+      simnfo->ecr        = the_globals->getECR();
-<
+    }
->
+  for(i=0; i<nInfo; i++){
->
+    info[i].nComponents = n_components;
->
+    info[i].componentsNmol = components_nmol;
->
+    info[i].compStamps = comp_stamps;
->
+    info[i].headStamp = headStamp;
->
+  }
->
-<
+    if( !the_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;
-<
+    } else {
-<
+      simnfo->est        = the_globals->getEST();
-<
+    }
->
+  for( i=0; i<n_components; i++ ){
->
->
+    id = the_components[i]->getType();
->
+    comp_stamps[i] = NULL;
-<
+    if(!the_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 {
-<
+    if (usesDipoles) {
->
+    // check to make sure the component isn't already in the list
->
->
+    comp_stamps[i] = headStamp->match( id );
->
+    if( comp_stamps[i] == NULL ){
-<
+      if( !the_globals->haveECR() ){
-<
+        sprintf( painCave.errMsg,
-<
+                 "SimSetup Warning: using default value of 1/2 the smallest "
-<
+                 "box length for the electrostaticCutoffRadius.\n"
-<
+                 "I hope you have a very fast processor!\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;
-<
+      } else {
-<
+        simnfo->ecr        = the_globals->getECR();
-<
+      }
->
+      // extract the component from the list;
-<
+      if( !the_globals->haveEST() ){
-<
+        sprintf( painCave.errMsg,
-<
+                 "SimSetup Warning: using default value of 5%% of the "
-<
+                 "electrostaticCutoffRadius for the "
-<
+                 "electrostaticSkinThickness\n"
-<
+                 );
-<
+        painCave.isFatal = 0;
-<
+        simError();
-<
+        simnfo->est = 0.05 * simnfo->ecr;
-<
+      } else {
-<
+        simnfo->est        = the_globals->getEST();
->
+      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 );
-<
+  }
->
+  }
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "electrostatic parameters check out" );
->
+  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
+  MPIcheckPoint();
+#endif // is_mpi
-–
+ if( the_globals->haveInitialConfig() ){
-–
-–
+     InitializeFromFile* fileInit;
-–
+#ifdef IS_MPI // is_mpi
-–
+     if( worldRank == 0 ){
-–
+#endif //is_mpi
-–
+   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
-–
+#ifdef IS_MPI
-–
+     }else fileInit = new InitializeFromFile( NULL );
-–
+#endif
-–
+   fileInit->read_xyz( simnfo ); // default velocities on
-<
+   delete fileInit;
-<
+ }
-<
+ else{
->
+}
-<
+#ifdef IS_MPI
-<
-<
+  // no init from bass
->
+void SimSetup::calcSysValues( void ){
->
+  int i, j, k;
-<
+  sprintf( painCave.errMsg,
-<
+           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
-<
+  painCave.isFatal;
-<
+  simError();
->
+  int *molMembershipArray;
-<
+#else
-<
-<
+  initFromBass();
-<
-<
-<
+#endif
-<
+ }
->
+  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();
->
+  }
->
->
+  tot_SRI = tot_bonds + tot_bends + tot_torsions;
->
+  molMembershipArray = new int[tot_atoms];
->
->
+  for(i=0; i<nInfo; i++){
->
+    info[i].n_atoms = tot_atoms;
->
+    info[i].n_bonds = tot_bonds;
->
+    info[i].n_bends = tot_bends;
->
+    info[i].n_torsions = tot_torsions;
->
+    info[i].n_SRI = tot_SRI;
->
+    info[i].n_mol = tot_nmol;
->
->
+    info[i].molMembershipArray = molMembershipArray;
->
+  }
->
+}
+#ifdef IS_MPI
-<
+  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
->
->
+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[0].molMembershipArray[globalAtomIndex] = allMol;
->
+        globalAtomIndex++;
->
+      }
->
->
+      allMol++;
->
+    }
->
+  }
->
+  local_SRI = local_bonds + local_bends + local_torsions;
->
->
+  info[0].n_atoms = mpiSim->getMyNlocal();
->
->
+  if( local_atoms != info[0].n_atoms ){
->
+    sprintf( painCave.errMsg,
->
+             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
->
+             " localAtom (%d) are not equal.\n",
->
+             info[0].n_atoms,
->
+             local_atoms );
->
+    painCave.isFatal = 1;
->
+    simError();
->
+  }
->
->
+  info[0].n_bonds = local_bonds;
->
+  info[0].n_bends = local_bends;
->
+  info[0].n_torsions = local_torsions;
->
+  info[0].n_SRI = local_SRI;
->
+  info[0].n_mol = localMol;
->
->
+  strcpy( checkPointMsg, "Passed nlocal consistency check." );
+  MPIcheckPoint();
-+
+}
-+
+#endif // is_mpi
-<
-<
-<
->
+void SimSetup::makeSysArrays( void ){
->
+  int i, j, k, l;
-+
+  Atom** the_atoms;
-+
+  Molecule* the_molecules;
-+
+  Exclude** the_excludes;
-+
-+
+  for(l=0; l<nInfo; l++){
-+
-+
+    // create the atom and short range interaction arrays
-+
-+
+    the_atoms = new Atom*[info[l].n_atoms];
-+
+    the_molecules = new Molecule[info[l].n_mol];
-+
+    int molIndex;
-+
-+
+    // initialize the molecule's stampID's
-+
+#ifdef IS_MPI
-–
+  if( worldRank == 0 ){
-–
+#endif // is_mpi
-<
+    if( the_globals->haveFinalConfig() ){
-<
+      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
-<
+    }
-<
+    else{
-<
+      strcpy( simnfo->finalName, inFileName );
-<
+      char* endTest;
-<
+      int nameLength = strlen( simnfo->finalName );
-<
+      endTest = &(simnfo->finalName[nameLength - 5]);
-<
+      if( !strcmp( endTest, ".bass" ) ){
-<
+        strcpy( endTest, ".eor" );
->
->
+    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 if( !strcmp( endTest, ".BASS" ) ){
-–
+        strcpy( endTest, ".eor" );
-–
+      }
-–
+      else{
-–
+        endTest = &(simnfo->finalName[nameLength - 4]);
-–
+        if( !strcmp( endTest, ".bss" ) ){
-–
+          strcpy( endTest, ".eor" );
-–
+        }
-–
+        else if( !strcmp( endTest, ".mdl" ) ){
-–
+          strcpy( endTest, ".eor" );
-–
+        }
-–
+        else{
-–
+          strcat( simnfo->finalName, ".eor" );
-–
+        }
-–
+      }
-<
+    // make the sample and status out names
->
+#else // is_mpi
-<
+    strcpy( simnfo->sampleName, inFileName );
-<
+    char* endTest;
-<
+    int nameLength = strlen( simnfo->sampleName );
-<
+    endTest = &(simnfo->sampleName[nameLength - 5]);
-<
+    if( !strcmp( endTest, ".bass" ) ){
-<
+      strcpy( endTest, ".dump" );
-<
+    }
-<
+    else if( !strcmp( endTest, ".BASS" ) ){
-<
+      strcpy( endTest, ".dump" );
-<
+    }
-<
+    else{
-<
+      endTest = &(simnfo->sampleName[nameLength - 4]);
-<
+      if( !strcmp( endTest, ".bss" ) ){
-<
+        strcpy( endTest, ".dump" );
->
+    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[l].molMembershipArray[globalAtomIndex] = molIndex;
->
+          globalAtomIndex++;
->
+        }
->
+        molIndex++;
-–
+      else if( !strcmp( endTest, ".mdl" ) ){
-–
+        strcpy( endTest, ".dump" );
-–
+      }
-–
+      else{
-–
+        strcat( simnfo->sampleName, ".dump" );
-–
+      }
-<
+    strcpy( simnfo->statusName, inFileName );
-<
+    nameLength = strlen( simnfo->statusName );
-<
+    endTest = &(simnfo->statusName[nameLength - 5]);
-<
+    if( !strcmp( endTest, ".bass" ) ){
-<
+      strcpy( endTest, ".stat" );
-<
+    }
-<
+    else if( !strcmp( endTest, ".BASS" ) ){
-<
+      strcpy( endTest, ".stat" );
-<
+    }
-<
+    else{
-<
+      endTest = &(simnfo->statusName[nameLength - 4]);
-<
+      if( !strcmp( endTest, ".bss" ) ){
-<
+        strcpy( endTest, ".stat" );
->
->
+#endif // is_mpi
->
->
->
+    if( info[l].n_SRI ){
->
->
+      Exclude::createArray(info[l].n_SRI);
->
+      the_excludes = new Exclude*[info[l].n_SRI];
->
+      for( int ex=0; ex<info[l].n_SRI; ex++){
->
+        the_excludes[ex] = new Exclude(ex);
-<
+      else if( !strcmp( endTest, ".mdl" ) ){
-<
+        strcpy( endTest, ".stat" );
-<
+      }
-<
+      else{
-<
+        strcat( simnfo->statusName, ".stat" );
-<
+      }
->
+      info[l].globalExcludes = new int;
->
+      info[l].n_exclude = info[l].n_SRI;
-+
+    else{
-<
+#ifdef IS_MPI
-<
+  }
-<
+#endif // is_mpi
-<
-<
+  // set the status, sample, and themal kick times
-<
-<
+  if( the_globals->haveSampleTime() ){
-<
+    simnfo->sampleTime = the_globals->getSampleTime();
-<
+    simnfo->statusTime = simnfo->sampleTime;
-<
+    simnfo->thermalTime = simnfo->sampleTime;
-<
+  }
-<
+  else{
-<
+    simnfo->sampleTime = the_globals->getRunTime();
-<
+    simnfo->statusTime = simnfo->sampleTime;
-<
+    simnfo->thermalTime = simnfo->sampleTime;
-<
+  }
->
+      Exclude::createArray( 1 );
->
+      the_excludes = new Exclude*;
->
+      the_excludes[0] = new Exclude(0);
->
+      the_excludes[0]->setPair( 0,0 );
->
+      info[l].globalExcludes = new int;
->
+      info[l].globalExcludes[0] = 0;
->
+      info[l].n_exclude = 0;
->
+    }
-<
+  if( the_globals->haveStatusTime() ){
-<
+    simnfo->statusTime = the_globals->getStatusTime();
-<
+  }
->
+    // set the arrays into the SimInfo object
-<
+  if( the_globals->haveThermalTime() ){
-<
+    simnfo->thermalTime = the_globals->getThermalTime();
->
+    info[l].atoms = the_atoms;
->
+    info[l].molecules = the_molecules;
->
+    info[l].nGlobalExcludes = 0;
->
+    info[l].excludes = the_excludes;
->
->
+    the_ff->setSimInfo( info );
->
-+
+}
-<
+  // check for the temperature set flag
->
+void SimSetup::makeIntegrator( void ){
-<
+  if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
->
+  int k;
-+
+  NVT<RealIntegrator>*  myNVT = NULL;
-+
+  NPTi<RealIntegrator>* myNPTi = NULL;
-+
+  NPTf<RealIntegrator>* myNPTf = NULL;
-+
+  NPTim<RealIntegrator>* myNPTim = NULL;
-+
+  NPTfm<RealIntegrator>* myNPTfm = NULL;
-+
-+
+  for(k=0; k<nInfo; k++){
-+
-+
+    switch( ensembleCase ){
-+
-+
+    case NVE_ENS:
-+
+        if (haveZConstraint){
-+
+         setupZConstraint();
-+
+           new ZConstraint<NVE<RealIntegrator> >( &(info[k]), the_ff );
-+
+        }
-<
+  // make the integrator
-<
-<
-<
+  NVT* myNVT = NULL;
-<
+  switch( ensembleCase ){
->
+        else
->
+        new NVE<RealIntegrator>( &(info[k]), the_ff );
->
+      break;
->
->
+    case NVT_ENS:
->
+        if (haveZConstraint){
->
+         setupZConstraint();
->
+           myNVT = new ZConstraint<NVT<RealIntegrator> >( &(info[k]), the_ff );
->
+        }
->
+        else
->
+        myNVT = new NVT<RealIntegrator>( &(info[k]), the_ff );
-<
+  case NVE_ENS:
-<
+    new NVE( simnfo, the_ff );
-<
+    break;
->
+      myNVT->setTargetTemp(globals->getTargetTemp());
->
->
+      if (globals->haveTauThermostat())
->
+        myNVT->setTauThermostat(globals->getTauThermostat());
->
->
+      else {
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use the NVT\n"
->
+                 "    ensemble, you must set tauThermostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
+      break;
->
->
+    case NPTi_ENS:
->
+        if (haveZConstraint){
->
+         setupZConstraint();
->
+           myNPTi = new ZConstraint<NPTi<RealIntegrator> >( &(info[k]), the_ff );
->
+        }
->
+        else
->
+        myNPTi = new NPTi<RealIntegrator>( &(info[k]), the_ff );
-<
+  case NVT_ENS:
-<
+    myNVT = new NVT( simnfo, the_ff );
-<
+    myNVT->setTargetTemp(the_globals->getTargetTemp());
-<
-<
+    if (the_globals->haveTauThermostat())
-<
+      myNVT->setTauThermostat(the_globals->getTauThermostat());
-<
+//     else if (the_globals->haveQmass())
-<
+//       myNVT->setQmass(the_globals->getQmass());
-<
+    else {
->
+        myNPTi->setTargetTemp( globals->getTargetTemp() );
->
->
+      if (globals->haveTargetPressure())
->
+        myNPTi->setTargetPressure(globals->getTargetPressure());
->
+      else {
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use a constant pressure\n"
->
+                 "    ensemble, you must set targetPressure in the BASS file.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauThermostat() )
->
+        myNPTi->setTauThermostat( globals->getTauThermostat() );
->
+      else{
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use an NPT\n"
->
+                 "    ensemble, you must set tauThermostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauBarostat() )
->
+        myNPTi->setTauBarostat( globals->getTauBarostat() );
->
+      else{
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use an NPT\n"
->
+                 "    ensemble, you must set tauBarostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
+      break;
->
->
+    case NPTf_ENS:
->
+        if (haveZConstraint){
->
+         setupZConstraint();
->
+           myNPTf = new ZConstraint<NPTf<RealIntegrator> >( &(info[k]), the_ff );
->
+        }
->
+        else
->
+        myNPTf = new NPTf<RealIntegrator>( &(info[k]), the_ff );
->
->
+      myNPTf->setTargetTemp( globals->getTargetTemp());
->
->
+      if (globals->haveTargetPressure())
->
+        myNPTf->setTargetPressure(globals->getTargetPressure());
->
+      else {
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use a constant pressure\n"
->
+                 "    ensemble, you must set targetPressure in the BASS file.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauThermostat() )
->
+        myNPTf->setTauThermostat( globals->getTauThermostat() );
->
+      else{
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use an NPT\n"
->
+               "    ensemble, you must set tauThermostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauBarostat() )
->
+        myNPTf->setTauBarostat( globals->getTauBarostat() );
->
+      else{
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use an NPT\n"
->
+                 "    ensemble, you must set tauBarostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
+      break;
->
->
+    case NPTim_ENS:
->
+        if (haveZConstraint){
->
+         setupZConstraint();
->
+           myNPTim = new ZConstraint<NPTim<RealIntegrator> >( &(info[k]), the_ff );
->
+        }
->
+        else
->
+        myNPTim = new NPTim<RealIntegrator>( &(info[k]), the_ff );
->
->
+        myNPTim->setTargetTemp( globals->getTargetTemp());
->
->
+      if (globals->haveTargetPressure())
->
+        myNPTim->setTargetPressure(globals->getTargetPressure());
->
+      else {
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use a constant pressure\n"
->
+                 "    ensemble, you must set targetPressure in the BASS file.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauThermostat() )
->
+        myNPTim->setTauThermostat( globals->getTauThermostat() );
->
+      else{
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use an NPT\n"
->
+                 "    ensemble, you must set tauThermostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauBarostat() )
->
+        myNPTim->setTauBarostat( globals->getTauBarostat() );
->
+      else{
+      sprintf( painCave.errMsg,
-<
+               "SimSetup error: If you use the NVT\n"
-<
+               "    ensemble, you must set either tauThermostat or qMass.\n"
-<
+               "    Neither of these was found in the BASS file.\n");
->
+               "SimSetup error: If you use an NPT\n"
->
+               "    ensemble, you must set tauBarostat.\n");
+      painCave.isFatal = 1;
+      simError();
-<
+    }
-<
+    break;
->
+      }
->
+      break;
->
->
+    case NPTfm_ENS:
->
+        if (haveZConstraint){
->
+         setupZConstraint();
->
+           myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >( &(info[k]), the_ff );
->
+        }
->
+        else
->
+        myNPTfm = new NPTfm<RealIntegrator>( &(info[k]), the_ff );
-<
+  default:
-<
+    sprintf( painCave.errMsg,
-<
+             "SimSetup Error. Unrecognized ensemble in case statement.\n");
-<
+    painCave.isFatal = 1;
-<
+    simError();
->
+        myNPTfm->setTargetTemp( globals->getTargetTemp());
->
->
+      if (globals->haveTargetPressure())
->
+        myNPTfm->setTargetPressure(globals->getTargetPressure());
->
+      else {
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use a constant pressure\n"
->
+                 "    ensemble, you must set targetPressure in the BASS file.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauThermostat() )
->
+        myNPTfm->setTauThermostat( globals->getTauThermostat() );
->
+      else{
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use an NPT\n"
->
+                 "    ensemble, you must set tauThermostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
->
+      if( globals->haveTauBarostat() )
->
+        myNPTfm->setTauBarostat( globals->getTauBarostat() );
->
+      else{
->
+        sprintf( painCave.errMsg,
->
+                 "SimSetup error: If you use an NPT\n"
->
+                 "    ensemble, you must set tauBarostat.\n");
->
+        painCave.isFatal = 1;
->
+        simError();
->
+      }
->
+      break;
->
->
+    default:
->
+      sprintf( painCave.errMsg,
->
+               "SimSetup Error. Unrecognized ensemble in case statement.\n");
->
+      painCave.isFatal = 1;
->
+      simError();
->
+    }
-+
+}
-+
+void SimSetup::initFortran( void ){
-<
+#ifdef IS_MPI
-<
+  mpiSim->mpiRefresh();
-<
+#endif
-<
-<
+  // initialize the Fortran
-<
-<
-<
+  simnfo->refreshSim();
->
+  info[0].refreshSim();
-<
+  if( !strcmp( simnfo->mixingRule, "standard") ){
->
+  if( !strcmp( info[0].mixingRule, "standard") ){
+    the_ff->initForceField( LB_MIXING_RULE );
-<
+  else if( !strcmp( simnfo->mixingRule, "explicit") ){
->
+  else if( !strcmp( info[0].mixingRule, "explicit") ){
+    the_ff->initForceField( EXPLICIT_MIXING_RULE );
+  else{
+    sprintf( painCave.errMsg,
+             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
-<
+             simnfo->mixingRule );
->
+             info[0].mixingRule );
+    painCave.isFatal = 1;
+    simError();
+          "Successfully intialized the mixingRule for Fortran." );
+  MPIcheckPoint();
+#endif // is_mpi
-+
-+
+void SimSetup::setupZConstraint()
-+
+{
-+
+  int k;
-<
+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++){
->
+  for(k=0; k<nInfo; k++){
-<
+    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++){
->
+    if(globals->haveZConsTime()){
-<
+      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() );
->
+      //add sample time of z-constraint  into SimInfo's property list
->
+      DoubleData* zconsTimeProp = new DoubleData();
->
+      zconsTimeProp->setID("zconstime");
->
+      zconsTimeProp->setData(globals->getZConsTime());
->
+      info[k].addProperty(zconsTimeProp);
->
+    }
->
+    else{
->
+      sprintf( painCave.errMsg,
->
+               "ZConstraint error: If you use an ZConstraint\n"
->
+               " , you must set sample time.\n");
->
+      painCave.isFatal = 1;
->
+      simError();
->
+    }
-<
+#ifdef IS_MPI
->
+    if(globals->haveIndexOfAllZConsMols()){
->
->
+      //add index of z-constraint molecules into SimInfo's property list
->
+      vector<int> tempIndex = globals->getIndexOfAllZConsMols();
-<
+      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
->
+      //sort the index
->
+      sort(tempIndex.begin(), tempIndex.end());
-<
+#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
->
+      IndexData* zconsIndex = new IndexData();
->
+      zconsIndex->setID("zconsindex");
->
+      zconsIndex->setIndexData(tempIndex);
->
+      info[k].addProperty(zconsIndex);
-<
+    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;
-<
+      }
->
+    else{
->
+      sprintf( painCave.errMsg,
->
+               "SimSetup error: If you use an ZConstraint\n"
->
+               " , you must set index of z-constraint molecules.\n");
->
+      painCave.isFatal = 1;
->
+      simError();
-–
+      // 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;
->
+    //Determine the name of ouput file and add it into SimInfo's property list
->
+    //Be careful, do not use inFileName, since it is a pointer which
->
+    //point to a string at master node, and slave nodes do not contain that string
->
->
+    string zconsOutput(info[k].finalName);
->
->
+    zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
->
->
+    StringData* zconsFilename = new StringData();
->
+    zconsFilename->setID("zconsfilename");
->
+    zconsFilename->setData(zconsOutput);
->
->
+    info[k].addProperty(zconsFilename);
-–
-–
+#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 582 by mmeineke, Wed Jul 9 15:33:46 2003 UTC vs. Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 582 by mmeineke, Wed Jul 9 15:33:46 2003 UTC vs.
Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC