mdtools/interface_implementation/SimSetup.cpp

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

#include "SimSetup.hpp"
#include "parse_me.h"
#include "LRI.hpp"
#include "Integrator.hpp"
#include "simError.h"

#ifdef IS_MPI
#include "mpiBASS.h"
#include "bassDiag.hpp"
#endif

SimSetup::SimSetup(){
  stamps = new MakeStamps();
  globals = new Globals();
  
#ifdef IS_MPI
  strcpy( checkPointMsg, "SimSetup creation successful" );
  MPIcheckPoint();
#endif // IS_MPI
}

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

void SimSetup::parseFile( char* fileName ){

#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    inFileName = fileName;
    set_interface_stamps( stamps, globals );
    
#ifdef IS_MPI
    mpiEventInit();
#endif

    yacc_BASS( fileName );

#ifdef IS_MPI
    throwMPIEvent(NULL);
  }
  else receiveParse();
#endif

}

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

    set_interface_stamps( stamps, globals );
    mpiEventInit();
    MPIcheckPoint();
    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{
    sprintf( painCave.errMsg,
             "SimSetup Error. Unrecognized force field -> %s\n",
             force_field );
    painCave.isFatal = 1;
    simError();
  }

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