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