mdtools/interface_implementation/SimSetup.cpp

#include <cstdlib>
#include <iostream>
#include <cmath>

#include "SimSetup.hpp"
#include "parse_me.h"
#include "LRI.hpp"
#include "Integrator.hpp"
#ifdef IS_MPI
#include "mpiInterface.h"
#include "bassDiag.hpp"
#endif

SimSetup::SimSetup(){
  stamps = new MakeStamps();
  globals = new Globals();
}

SimSetup::~SimSetup(){
  delete stamps;
  delete globals;
}

void SimSetup::parseFile( char* fileName ){

  inFileName = fileName;
  set_interface_stamps( stamps, globals );
#ifdef IS_MPI
  mpiEventInit();
#endif
  yacc_BASS( fileName );
#ifdef IS_MPI
  throwMPIEvent(NULL);
#endif

}

#ifdef IS_MPI
void SimSetup::receiveParse(void){

    set_interface_stamps( stamps, globals );
    mpiEventInit();
    mpiEventLoop();

}


void SimSetup::testMe(void){
  bassDiag* dumpMe = new bassDiag(globals,stamps);
  dumpMe->dumpStamps();
  delete dumpMe;
}
#endif
void SimSetup::createSim( void ){

  MakeStamps *the_stamps;
  Globals* the_globals;
  int i;

  // 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() );
  strcpy( ensemble, the_globals->getEnsemble() );

  if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
  else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
  else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
  else{
    std::cerr<< "SimSetup Error. Unrecognized force field -> "
             << force_field << "\n";
    exit(8);
  }

  // get the components and calculate the tot_nMol and indvidual n_mol
  the_components = the_globals->getComponents();
  components_nmol = new int[n_components];
  comp_stamps = new MoleculeStamp*[n_components];

  if( !the_globals->haveNMol() ){
    // we don't have the total number of molecules, so we assume it is
    // given in each component

    tot_nmol = 0;
    for( i=0; i<n_components; i++ ){

      if( !the_components[i]->haveNMol() ){
        // we have a problem
        std::cerr << "SimSetup Error. No global NMol or component NMol"
                  << " given. Cannot calculate the number of atoms.\n";
        exit( 8 );
      }

      tot_nmol += the_components[i]->getNMol();
      components_nmol[i] = the_components[i]->getNMol();
    }
  }
  else{
    std::cerr << "NOT A SUPPORTED FEATURE\n";

//     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

  }

  // make an array of molecule stamps that match the components used.

  for( i=0; i<n_components; i++ ){

    comp_stamps[i] =
      the_stamps->getMolecule( the_components[i]->getType() );
  }


  // 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();
  }

  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;

  // create the atom and short range interaction arrays

  the_atoms = new Atom*[tot_atoms];
  the_molecules = new Molecule[tot_nmol];


  if( tot_SRI ){
    the_sris = new SRI*[tot_SRI];
    the_excludes = new ex_pair[tot_SRI];
  }

  // set the arrays into the SimInfo object

  simnfo->atoms = the_atoms;
  simnfo->sr_interactions = the_sris;
  simnfo->n_exclude = tot_SRI;
  simnfo->excludes = the_excludes;


  // initialize the arrays

  the_ff->setSimInfo( simnfo );

  makeAtoms();

  if( tot_bonds ){
    makeBonds();
  }

  if( tot_bends ){
    makeBends();
  }

  if( tot_torsions ){
    makeTorsions();
  }

  //  makeMolecules();

  // get some of the tricky things that may still be in the globals

  if( simnfo->n_dipoles ){

    if( !the_globals->haveRRF() ){
      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
      exit(8);
    }
    if( !the_globals->haveDielectric() ){
      std::cerr << "SimSetup Error, system has dipoles, but no"
                << " dielectric was set.\n";
      exit(8);
    }

    simnfo->rRF        = the_globals->getRRF();
    simnfo->dielectric = the_globals->getDielectric();
  }

  if( the_globals->haveBox() ){
    simnfo->box_x = the_globals->getBox();
    simnfo->box_y = the_globals->getBox();
    simnfo->box_z = the_globals->getBox();
  }
  else if( the_globals->haveDensity() ){

    double vol;
    vol = (double)tot_nmol / the_globals->getDensity();
    simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
    simnfo->box_y = simnfo->box_x;
    simnfo->box_z = simnfo->box_x;
  }
  else{
    if( !the_globals->haveBoxX() ){
      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
      exit(8);
    }
    simnfo->box_x = the_globals->getBoxX();

    if( !the_globals->haveBoxY() ){
      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
      exit(8);
    }
    simnfo->box_y = the_globals->getBoxY();

    if( !the_globals->haveBoxZ() ){
      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
      exit(8);
    }
    simnfo->box_z = the_globals->getBoxZ();
  }


//   if( the_globals->haveInitialConfig() ){
//        InitializeFromFile* fileInit;
//     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );

//     fileInit->read_xyz( simnfo ); // default velocities on

//     delete fileInit;
//   }
//   else{

    initFromBass();


//   }

//   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" );
//     }
//     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

//   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" );
//     }
//     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" );
//     }
//     else if( !strcmp( endTest, ".mdl" ) ){
//       strcpy( endTest, ".stat" );
//     }
//     else{
//       strcat( simnfo->statusName, ".stat" );
//     }
//   }


  // 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;
  }

  if( the_globals->haveStatusTime() ){
    simnfo->statusTime = the_globals->getStatusTime();
  }

  if( the_globals->haveThermalTime() ){
    simnfo->thermalTime = the_globals->getThermalTime();
  }

  // check for the temperature set flag

  if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();


  // make the longe range forces and the integrator

  new AllLong( simnfo );

  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
  if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
  if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
}

void SimSetup::makeAtoms( void ){

  int i, j, k, index;
  double ux, uy, uz, uSqr, u;
  AtomStamp* current_atom;
  DirectionalAtom* dAtom;
  int molIndex, molStart, molEnd, nMemb;


  molIndex = 0;
  index = 0;
  for( i=0; i<n_components; i++ ){

    for( j=0; j<components_nmol[i]; j++ ){

      molStart = index;
      nMemb = comp_stamps[i]->getNAtoms();
      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){

        current_atom = comp_stamps[i]->getAtom( k );
        if( current_atom->haveOrientation() ){

          dAtom = new DirectionalAtom;
          simnfo->n_oriented++;
          the_atoms[index] = dAtom;

          ux = current_atom->getOrntX();
          uy = current_atom->getOrntY();
          uz = current_atom->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{
          the_atoms[index] = new GeneralAtom;
        }
        the_atoms[index]->setType( current_atom->getType() );
        the_atoms[index]->setIndex( index );

        // increment the index and repeat;
        index++;
      }

      molEnd = index -1;
      the_molecules[molIndex].setNMembers( nMemb );
      the_molecules[molIndex].setStartAtom( molStart );
      the_molecules[molIndex].setEndAtom( molEnd );
      molIndex++;

    }
  }

  the_ff->initializeAtoms();
}

void SimSetup::makeBonds( void ){

  int i, j, k, index, offset;
  bond_pair* the_bonds;
  BondStamp* current_bond;

  the_bonds = new bond_pair[tot_bonds];
  index = 0;
  offset = 0;
  for( i=0; i<n_components; i++ ){

    for( j=0; j<components_nmol[i]; j++ ){

      for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){

        current_bond = comp_stamps[i]->getBond( k );
        the_bonds[index].a = current_bond->getA() + offset;
        the_bonds[index].b = current_bond->getB() + offset;

        the_excludes[index].i = the_bonds[index].a;
        the_excludes[index].j = the_bonds[index].b;

        // increment the index and repeat;
        index++;
      }
      offset += comp_stamps[i]->getNAtoms();
    }
  }

  the_ff->initializeBonds( the_bonds );
}

void SimSetup::makeBends( void ){

  int i, j, k, index, offset;
  bend_set* the_bends;
  BendStamp* current_bend;

  the_bends = new bend_set[tot_bends];
  index = 0;
  offset = 0;
  for( i=0; i<n_components; i++ ){

    for( j=0; j<components_nmol[i]; j++ ){

      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){

        current_bend = comp_stamps[i]->getBend( k );
        the_bends[index].a = current_bend->getA() + offset;
        the_bends[index].b = current_bend->getB() + offset;
        the_bends[index].c = current_bend->getC() + offset;

        the_excludes[index + tot_bonds].i = the_bends[index].a;
        the_excludes[index + tot_bonds].j = the_bends[index].c;

        // increment the index and repeat;
        index++;
      }
      offset += comp_stamps[i]->getNAtoms();
    }
  }

  the_ff->initializeBends( the_bends );
}

void SimSetup::makeTorsions( void ){

  int i, j, k, index, offset;
  torsion_set* the_torsions;
  TorsionStamp* current_torsion;

  the_torsions = new torsion_set[tot_torsions];
  index = 0;
  offset = 0;
  for( i=0; i<n_components; i++ ){

    for( j=0; j<components_nmol[i]; j++ ){

      for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){

        current_torsion = comp_stamps[i]->getTorsion( k );
        the_torsions[index].a = current_torsion->getA() + offset;
        the_torsions[index].b = current_torsion->getB() + offset;
        the_torsions[index].c = current_torsion->getC() + offset;
        the_torsions[index].d = current_torsion->getD() + offset;

        the_excludes[index + tot_bonds + tot_bends].i = the_torsions[index].a;
        the_excludes[index + tot_bonds + tot_bends].j = the_torsions[index].d;

        // increment the index and repeat;
        index++;
      }
      offset += comp_stamps[i]->getNAtoms();
    }
  }

  the_ff->initializeTorsions( the_torsions );
}

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->box_x / temp3;
    celly = simnfo->box_y / temp3;
    cellz = simnfo->box_z / 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){
      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
      exit(8);
    }
  }
  else{
    n_cells = (int)temp3;
    cellx = simnfo->box_x / temp3;
    celly = simnfo->box_y / temp3;
    cellz = simnfo->box_z / 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() ){
      std::cerr << "Component " << comp_stamps[current_comp]->getID()
                << ", atom " << current_atom->getType()
                << " does not have a position specified.\n"
                << "The initialization routine is unable to give a start"
                << " position.\n";
      exit(8);
    }

    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++;
  }
}
Revision:	131
Committed:	Wed Oct 9 22:29:40 2002 UTC (21 years, 9 months ago) by chuckv
File size:	17128 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	mmeineke	10	#include <cstdlib>
2			#include <iostream>
3			#include <cmath>
4
5			#include "SimSetup.hpp"
6			#include "parse_me.h"
7			#include "LRI.hpp"
8			#include "Integrator.hpp"
9	chuckv	131	#ifdef IS_MPI
10	chuckv	124	#include "mpiInterface.h"
11	chuckv	131	#include "bassDiag.hpp"
12			#endif
13	mmeineke	10
14			SimSetup::SimSetup(){
15			stamps = new MakeStamps();
16			globals = new Globals();
17			}
18
19			SimSetup::~SimSetup(){
20			delete stamps;
21			delete globals;
22			}
23
24			void SimSetup::parseFile( char* fileName ){
25
26			inFileName = fileName;
27			set_interface_stamps( stamps, globals );
28	chuckv	131	#ifdef IS_MPI
29	chuckv	124	mpiEventInit();
30			#endif
31	mmeineke	10	yacc_BASS( fileName );
32	chuckv	131	#ifdef IS_MPI
33	chuckv	124	throwMPIEvent(NULL);
34			#endif
35
36	mmeineke	10	}
37
38	chuckv	131	#ifdef IS_MPI
39	chuckv	124	void SimSetup::receiveParse(void){
40
41			set_interface_stamps( stamps, globals );
42			mpiEventInit();
43			mpiEventLoop();
44
45			}
46
47	chuckv	131
48	chuckv	128	void SimSetup::testMe(void){
49			bassDiag* dumpMe = new bassDiag(globals,stamps);
50			dumpMe->dumpStamps();
51			delete dumpMe;
52			}
53	chuckv	131	#endif
54	mmeineke	10	void SimSetup::createSim( void ){
55
56			MakeStamps *the_stamps;
57			Globals* the_globals;
58			int i;
59
60			// get the stamps and globals;
61			the_stamps = stamps;
62			the_globals = globals;
63
64			// set the easy ones first
65			simnfo->target_temp = the_globals->getTargetTemp();
66			simnfo->dt = the_globals->getDt();
67			simnfo->run_time = the_globals->getRunTime();
68
69			// get the ones we know are there, yet still may need some work.
70			n_components = the_globals->getNComponents();
71			strcpy( force_field, the_globals->getForceField() );
72			strcpy( ensemble, the_globals->getEnsemble() );
73	chuckv	124
74	mmeineke	10	if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
75			else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
76			else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
77			else{
78	chuckv	124	std::cerr<< "SimSetup Error. Unrecognized force field -> "
79	mmeineke	10	<< force_field << "\n";
80			exit(8);
81			}
82
83			// get the components and calculate the tot_nMol and indvidual n_mol
84			the_components = the_globals->getComponents();
85			components_nmol = new int[n_components];
86			comp_stamps = new MoleculeStamp*[n_components];
87	chuckv	124
88	mmeineke	10	if( !the_globals->haveNMol() ){
89	chuckv	124	// we don't have the total number of molecules, so we assume it is
90	mmeineke	10	// given in each component
91
92			tot_nmol = 0;
93			for( i=0; i<n_components; i++ ){
94	chuckv	124
95	mmeineke	10	if( !the_components[i]->haveNMol() ){
96			// we have a problem
97			std::cerr << "SimSetup Error. No global NMol or component NMol"
98			<< " given. Cannot calculate the number of atoms.\n";
99			exit( 8 );
100			}
101
102			tot_nmol += the_components[i]->getNMol();
103			components_nmol[i] = the_components[i]->getNMol();
104			}
105			}
106			else{
107			std::cerr << "NOT A SUPPORTED FEATURE\n";
108	chuckv	124
109	mmeineke	10	// tot_nmol = the_globals->getNMol();
110	chuckv	124
111	mmeineke	10	// //we have the total number of molecules, now we check for molfractions
112			// for( i=0; i<n_components; i++ ){
113	chuckv	124
114	mmeineke	10	// if( !the_components[i]->haveMolFraction() ){
115	chuckv	124
116	mmeineke	10	// if( !the_components[i]->haveNMol() ){
117			// //we have a problem
118			// std::cerr << "SimSetup error. Neither molFraction nor "
119			// << " nMol was given in component
120
121			}
122
123			// make an array of molecule stamps that match the components used.
124
125			for( i=0; i<n_components; i++ ){
126
127	chuckv	124	comp_stamps[i] =
128	mmeineke	10	the_stamps->getMolecule( the_components[i]->getType() );
129			}
130
131
132	chuckv	124
133	mmeineke	10	// caclulate the number of atoms, bonds, bends and torsions
134
135			tot_atoms = 0;
136			tot_bonds = 0;
137			tot_bends = 0;
138			tot_torsions = 0;
139			for( i=0; i<n_components; i++ ){
140	chuckv	124
141	mmeineke	10	tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
142			tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
143			tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
144			tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
145			}
146	chuckv	124
147	mmeineke	10	tot_SRI = tot_bonds + tot_bends + tot_torsions;
148	chuckv	124
149	mmeineke	10	simnfo->n_atoms = tot_atoms;
150			simnfo->n_bonds = tot_bonds;
151			simnfo->n_bends = tot_bends;
152			simnfo->n_torsions = tot_torsions;
153			simnfo->n_SRI = tot_SRI;
154
155			// create the atom and short range interaction arrays
156	chuckv	124
157	mmeineke	10	the_atoms = new Atom*[tot_atoms];
158	mmeineke	113	the_molecules = new Molecule[tot_nmol];
159	chuckv	124
160
161	mmeineke	10	if( tot_SRI ){
162			the_sris = new SRI*[tot_SRI];
163			the_excludes = new ex_pair[tot_SRI];
164			}
165
166			// set the arrays into the SimInfo object
167
168			simnfo->atoms = the_atoms;
169			simnfo->sr_interactions = the_sris;
170			simnfo->n_exclude = tot_SRI;
171			simnfo->excludes = the_excludes;
172
173	chuckv	124
174	mmeineke	10	// initialize the arrays
175	chuckv	124
176	mmeineke	10	the_ff->setSimInfo( simnfo );
177	chuckv	124
178	mmeineke	10	makeAtoms();
179
180			if( tot_bonds ){
181			makeBonds();
182			}
183
184			if( tot_bends ){
185			makeBends();
186			}
187
188			if( tot_torsions ){
189			makeTorsions();
190			}
191
192			// makeMolecules();
193
194			// get some of the tricky things that may still be in the globals
195
196			if( simnfo->n_dipoles ){
197
198			if( !the_globals->haveRRF() ){
199			std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
200			exit(8);
201			}
202			if( !the_globals->haveDielectric() ){
203			std::cerr << "SimSetup Error, system has dipoles, but no"
204			<< " dielectric was set.\n";
205			exit(8);
206			}
207
208			simnfo->rRF = the_globals->getRRF();
209			simnfo->dielectric = the_globals->getDielectric();
210			}
211
212			if( the_globals->haveBox() ){
213			simnfo->box_x = the_globals->getBox();
214			simnfo->box_y = the_globals->getBox();
215			simnfo->box_z = the_globals->getBox();
216			}
217			else if( the_globals->haveDensity() ){
218	chuckv	124
219	mmeineke	10	double vol;
220			vol = (double)tot_nmol / the_globals->getDensity();
221			simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
222			simnfo->box_y = simnfo->box_x;
223			simnfo->box_z = simnfo->box_x;
224			}
225			else{
226			if( !the_globals->haveBoxX() ){
227			std::cerr << "SimSetup error, no periodic BoxX size given.\n";
228			exit(8);
229			}
230			simnfo->box_x = the_globals->getBoxX();
231
232			if( !the_globals->haveBoxY() ){
233			std::cerr << "SimSetup error, no periodic BoxY size given.\n";
234			exit(8);
235			}
236			simnfo->box_y = the_globals->getBoxY();
237
238			if( !the_globals->haveBoxZ() ){
239			std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
240			exit(8);
241			}
242			simnfo->box_z = the_globals->getBoxZ();
243			}
244
245	chuckv	124
246			// if( the_globals->haveInitialConfig() ){
247			// InitializeFromFile* fileInit;
248			// fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
249
250			// fileInit->read_xyz( simnfo ); // default velocities on
251
252			// delete fileInit;
253			// }
254			// else{
255
256	mmeineke	10	initFromBass();
257
258
259	chuckv	124	// }
260
261			// if( the_globals->haveFinalConfig() ){
262			// strcpy( simnfo->finalName, the_globals->getFinalConfig() );
263			// }
264			// else{
265			// strcpy( simnfo->finalName, inFileName );
266			// char* endTest;
267			// int nameLength = strlen( simnfo->finalName );
268			// endTest = &(simnfo->finalName[nameLength - 5]);
269			// if( !strcmp( endTest, ".bass" ) ){
270			// strcpy( endTest, ".eor" );
271			// }
272			// else if( !strcmp( endTest, ".BASS" ) ){
273			// strcpy( endTest, ".eor" );
274			// }
275			// else{
276			// endTest = &(simnfo->finalName[nameLength - 4]);
277			// if( !strcmp( endTest, ".bss" ) ){
278			// strcpy( endTest, ".eor" );
279			// }
280			// else if( !strcmp( endTest, ".mdl" ) ){
281			// strcpy( endTest, ".eor" );
282			// }
283			// else{
284			// strcat( simnfo->finalName, ".eor" );
285			// }
286			// }
287			// }
288
289			// // make the sample and status out names
290
291			// strcpy( simnfo->sampleName, inFileName );
292			// char* endTest;
293			// int nameLength = strlen( simnfo->sampleName );
294			// endTest = &(simnfo->sampleName[nameLength - 5]);
295			// if( !strcmp( endTest, ".bass" ) ){
296			// strcpy( endTest, ".dump" );
297			// }
298			// else if( !strcmp( endTest, ".BASS" ) ){
299			// strcpy( endTest, ".dump" );
300			// }
301			// else{
302			// endTest = &(simnfo->sampleName[nameLength - 4]);
303			// if( !strcmp( endTest, ".bss" ) ){
304			// strcpy( endTest, ".dump" );
305			// }
306			// else if( !strcmp( endTest, ".mdl" ) ){
307			// strcpy( endTest, ".dump" );
308			// }
309			// else{
310			// strcat( simnfo->sampleName, ".dump" );
311			// }
312			// }
313
314			// strcpy( simnfo->statusName, inFileName );
315			// nameLength = strlen( simnfo->statusName );
316			// endTest = &(simnfo->statusName[nameLength - 5]);
317			// if( !strcmp( endTest, ".bass" ) ){
318			// strcpy( endTest, ".stat" );
319			// }
320			// else if( !strcmp( endTest, ".BASS" ) ){
321			// strcpy( endTest, ".stat" );
322			// }
323			// else{
324			// endTest = &(simnfo->statusName[nameLength - 4]);
325			// if( !strcmp( endTest, ".bss" ) ){
326			// strcpy( endTest, ".stat" );
327			// }
328			// else if( !strcmp( endTest, ".mdl" ) ){
329			// strcpy( endTest, ".stat" );
330			// }
331			// else{
332			// strcat( simnfo->statusName, ".stat" );
333			// }
334			// }
335
336
337	mmeineke	10	// set the status, sample, and themal kick times
338
339			if( the_globals->haveSampleTime() ){
340	chuckv	124	simnfo->sampleTime = the_globals->getSampleTime();
341	mmeineke	10	simnfo->statusTime = simnfo->sampleTime;
342			simnfo->thermalTime = simnfo->sampleTime;
343			}
344			else{
345	chuckv	124	simnfo->sampleTime = the_globals->getRunTime();
346	mmeineke	10	simnfo->statusTime = simnfo->sampleTime;
347			simnfo->thermalTime = simnfo->sampleTime;
348			}
349
350			if( the_globals->haveStatusTime() ){
351			simnfo->statusTime = the_globals->getStatusTime();
352			}
353
354			if( the_globals->haveThermalTime() ){
355			simnfo->thermalTime = the_globals->getThermalTime();
356			}
357
358			// check for the temperature set flag
359
360			if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
361	chuckv	124
362
363	mmeineke	10	// make the longe range forces and the integrator
364	chuckv	124
365	mmeineke	10	new AllLong( simnfo );
366	chuckv	124
367	mmeineke	10	if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
368			if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
369			if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
370			}
371
372			void SimSetup::makeAtoms( void ){
373	chuckv	124
374	mmeineke	10	int i, j, k, index;
375			double ux, uy, uz, uSqr, u;
376			AtomStamp* current_atom;
377			DirectionalAtom* dAtom;
378	mmeineke	117	int molIndex, molStart, molEnd, nMemb;
379	mmeineke	10
380	chuckv	124
381	mmeineke	117	molIndex = 0;
382	mmeineke	10	index = 0;
383			for( i=0; i<n_components; i++ ){
384	chuckv	124
385	mmeineke	10	for( j=0; j<components_nmol[i]; j++ ){
386	chuckv	124
387	mmeineke	117	molStart = index;
388			nMemb = comp_stamps[i]->getNAtoms();
389	mmeineke	10	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
390	chuckv	124
391	mmeineke	10	current_atom = comp_stamps[i]->getAtom( k );
392	chuckv	124	if( current_atom->haveOrientation() ){
393	mmeineke	10
394			dAtom = new DirectionalAtom;
395			simnfo->n_oriented++;
396			the_atoms[index] = dAtom;
397	chuckv	124
398	mmeineke	10	ux = current_atom->getOrntX();
399			uy = current_atom->getOrntY();
400			uz = current_atom->getOrntZ();
401	chuckv	124
402	mmeineke	10	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
403	chuckv	124
404	mmeineke	10	u = sqrt( uSqr );
405			ux = ux / u;
406			uy = uy / u;
407			uz = uz / u;
408	chuckv	124
409	mmeineke	10	dAtom->setSUx( ux );
410			dAtom->setSUy( uy );
411			dAtom->setSUz( uz );
412			}
413			else{
414			the_atoms[index] = new GeneralAtom;
415			}
416			the_atoms[index]->setType( current_atom->getType() );
417			the_atoms[index]->setIndex( index );
418	chuckv	124
419	mmeineke	10	// increment the index and repeat;
420			index++;
421			}
422	chuckv	124
423	mmeineke	117	molEnd = index -1;
424			the_molecules[molIndex].setNMembers( nMemb );
425			the_molecules[molIndex].setStartAtom( molStart );
426			the_molecules[molIndex].setEndAtom( molEnd );
427			molIndex++;
428
429	mmeineke	10	}
430			}
431	chuckv	124
432	mmeineke	10	the_ff->initializeAtoms();
433			}
434
435			void SimSetup::makeBonds( void ){
436
437			int i, j, k, index, offset;
438			bond_pair* the_bonds;
439			BondStamp* current_bond;
440
441			the_bonds = new bond_pair[tot_bonds];
442			index = 0;
443			offset = 0;
444			for( i=0; i<n_components; i++ ){
445	chuckv	124
446	mmeineke	10	for( j=0; j<components_nmol[i]; j++ ){
447	chuckv	124
448	mmeineke	10	for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
449	chuckv	124
450	mmeineke	10	current_bond = comp_stamps[i]->getBond( k );
451			the_bonds[index].a = current_bond->getA() + offset;
452			the_bonds[index].b = current_bond->getB() + offset;
453
454			the_excludes[index].i = the_bonds[index].a;
455			the_excludes[index].j = the_bonds[index].b;
456
457			// increment the index and repeat;
458			index++;
459			}
460			offset += comp_stamps[i]->getNAtoms();
461			}
462			}
463	chuckv	124
464	mmeineke	10	the_ff->initializeBonds( the_bonds );
465			}
466
467			void SimSetup::makeBends( void ){
468
469			int i, j, k, index, offset;
470			bend_set* the_bends;
471			BendStamp* current_bend;
472
473			the_bends = new bend_set[tot_bends];
474			index = 0;
475			offset = 0;
476			for( i=0; i<n_components; i++ ){
477	chuckv	124
478	mmeineke	10	for( j=0; j<components_nmol[i]; j++ ){
479	chuckv	124
480	mmeineke	10	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
481	chuckv	124
482	mmeineke	10	current_bend = comp_stamps[i]->getBend( k );
483			the_bends[index].a = current_bend->getA() + offset;
484			the_bends[index].b = current_bend->getB() + offset;
485			the_bends[index].c = current_bend->getC() + offset;
486
487			the_excludes[index + tot_bonds].i = the_bends[index].a;
488			the_excludes[index + tot_bonds].j = the_bends[index].c;
489
490			// increment the index and repeat;
491			index++;
492			}
493			offset += comp_stamps[i]->getNAtoms();
494			}
495			}
496	chuckv	124
497	mmeineke	10	the_ff->initializeBends( the_bends );
498			}
499
500			void SimSetup::makeTorsions( void ){
501
502			int i, j, k, index, offset;
503			torsion_set* the_torsions;
504			TorsionStamp* current_torsion;
505
506			the_torsions = new torsion_set[tot_torsions];
507			index = 0;
508			offset = 0;
509			for( i=0; i<n_components; i++ ){
510	chuckv	124
511	mmeineke	10	for( j=0; j<components_nmol[i]; j++ ){
512	chuckv	124
513	mmeineke	10	for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
514	chuckv	124
515	mmeineke	10	current_torsion = comp_stamps[i]->getTorsion( k );
516			the_torsions[index].a = current_torsion->getA() + offset;
517			the_torsions[index].b = current_torsion->getB() + offset;
518			the_torsions[index].c = current_torsion->getC() + offset;
519			the_torsions[index].d = current_torsion->getD() + offset;
520
521			the_excludes[index + tot_bonds + tot_bends].i = the_torsions[index].a;
522			the_excludes[index + tot_bonds + tot_bends].j = the_torsions[index].d;
523
524			// increment the index and repeat;
525			index++;
526			}
527			offset += comp_stamps[i]->getNAtoms();
528			}
529			}
530	chuckv	124
531	mmeineke	10	the_ff->initializeTorsions( the_torsions );
532			}
533
534			void SimSetup::initFromBass( void ){
535
536			int i, j, k;
537			int n_cells;
538			double cellx, celly, cellz;
539			double temp1, temp2, temp3;
540			int n_per_extra;
541			int n_extra;
542			int have_extra, done;
543
544			temp1 = (double)tot_nmol / 4.0;
545			temp2 = pow( temp1, ( 1.0 / 3.0 ) );
546			temp3 = ceil( temp2 );
547
548			have_extra =0;
549			if( temp2 < temp3 ){ // we have a non-complete lattice
550			have_extra =1;
551
552			n_cells = (int)temp3 - 1;
553			cellx = simnfo->box_x / temp3;
554			celly = simnfo->box_y / temp3;
555			cellz = simnfo->box_z / temp3;
556			n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
557			temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
558			n_per_extra = (int)ceil( temp1 );
559
560			if( n_per_extra > 4){
561			std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
562			exit(8);
563			}
564			}
565			else{
566			n_cells = (int)temp3;
567			cellx = simnfo->box_x / temp3;
568			celly = simnfo->box_y / temp3;
569			cellz = simnfo->box_z / temp3;
570			}
571	chuckv	124
572	mmeineke	10	current_mol = 0;
573			current_comp_mol = 0;
574			current_comp = 0;
575			current_atom_ndx = 0;
576	chuckv	124
577	mmeineke	10	for( i=0; i < n_cells ; i++ ){
578			for( j=0; j < n_cells; j++ ){
579			for( k=0; k < n_cells; k++ ){
580	chuckv	124
581	mmeineke	10	makeElement( i * cellx,
582			j * celly,
583			k * cellz );
584	chuckv	124
585	mmeineke	10	makeElement( i * cellx + 0.5 * cellx,
586			j * celly + 0.5 * celly,
587			k * cellz );
588	chuckv	124
589	mmeineke	10	makeElement( i * cellx,
590			j * celly + 0.5 * celly,
591			k * cellz + 0.5 * cellz );
592	chuckv	124
593	mmeineke	10	makeElement( i * cellx + 0.5 * cellx,
594			j * celly,
595			k * cellz + 0.5 * cellz );
596			}
597			}
598			}
599
600			if( have_extra ){
601			done = 0;
602	chuckv	124
603	mmeineke	10	int start_ndx;
604			for( i=0; i < (n_cells+1) && !done; i++ ){
605			for( j=0; j < (n_cells+1) && !done; j++ ){
606	chuckv	124
607	mmeineke	10	if( i < n_cells ){
608	chuckv	124
609	mmeineke	10	if( j < n_cells ){
610			start_ndx = n_cells;
611			}
612			else start_ndx = 0;
613			}
614			else start_ndx = 0;
615	chuckv	124
616	mmeineke	10	for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
617	chuckv	124
618	mmeineke	10	makeElement( i * cellx,
619			j * celly,
620			k * cellz );
621			done = ( current_mol >= tot_nmol );
622	chuckv	124
623	mmeineke	10	if( !done && n_per_extra > 1 ){
624			makeElement( i * cellx + 0.5 * cellx,
625			j * celly + 0.5 * celly,
626			k * cellz );
627			done = ( current_mol >= tot_nmol );
628			}
629	chuckv	124
630	mmeineke	10	if( !done && n_per_extra > 2){
631			makeElement( i * cellx,
632			j * celly + 0.5 * celly,
633			k * cellz + 0.5 * cellz );
634			done = ( current_mol >= tot_nmol );
635			}
636	chuckv	124
637	mmeineke	10	if( !done && n_per_extra > 3){
638			makeElement( i * cellx + 0.5 * cellx,
639			j * celly,
640			k * cellz + 0.5 * cellz );
641			done = ( current_mol >= tot_nmol );
642			}
643			}
644			}
645			}
646			}
647	chuckv	124
648
649	mmeineke	10	for( i=0; i<simnfo->n_atoms; i++ ){
650			simnfo->atoms[i]->set_vx( 0.0 );
651			simnfo->atoms[i]->set_vy( 0.0 );
652			simnfo->atoms[i]->set_vz( 0.0 );
653			}
654			}
655
656			void SimSetup::makeElement( double x, double y, double z ){
657
658			int k;
659			AtomStamp* current_atom;
660			DirectionalAtom* dAtom;
661			double rotMat[3][3];
662
663			for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
664	chuckv	124
665	mmeineke	10	current_atom = comp_stamps[current_comp]->getAtom( k );
666			if( !current_atom->havePosition() ){
667			std::cerr << "Component " << comp_stamps[current_comp]->getID()
668			<< ", atom " << current_atom->getType()
669			<< " does not have a position specified.\n"
670			<< "The initialization routine is unable to give a start"
671			<< " position.\n";
672			exit(8);
673			}
674	chuckv	124
675	mmeineke	10	the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
676			the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
677			the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
678	chuckv	124
679	mmeineke	10	if( the_atoms[current_atom_ndx]->isDirectional() ){
680	chuckv	124
681	mmeineke	10	dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
682	chuckv	124
683	mmeineke	10	rotMat[0][0] = 1.0;
684			rotMat[0][1] = 0.0;
685			rotMat[0][2] = 0.0;
686
687			rotMat[1][0] = 0.0;
688			rotMat[1][1] = 1.0;
689			rotMat[1][2] = 0.0;
690
691			rotMat[2][0] = 0.0;
692			rotMat[2][1] = 0.0;
693			rotMat[2][2] = 1.0;
694
695			dAtom->setA( rotMat );
696			}
697
698			current_atom_ndx++;
699			}
700	chuckv	124
701	mmeineke	10	current_mol++;
702			current_comp_mol++;
703
704			if( current_comp_mol >= components_nmol[current_comp] ){
705	chuckv	124
706	mmeineke	10	current_comp_mol = 0;
707			current_comp++;
708			}
709			}