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
        sprintf( painCave.errMsg,
                 "SimSetup Error. No global NMol or component NMol"
                 " given. Cannot calculate the number of atoms.\n" );
        painCave.isFatal = 1;
        simError();
      }

      tot_nmol += the_components[i]->getNMol();
      components_nmol[i] = the_components[i]->getNMol();
    }
  }
  else{
    sprintf( painCave.errMsg,
             "SimSetup error.\n"
             "\tSorry, the ability to specify total"
             " nMols and then give molfractions in the components\n"
             "\tis not currently supported."
             " Please give nMol in the components.\n" );
    painCave.isFatal = 1;
    simError();
    
    
    //     tot_nmol = the_globals->getNMol();
    
    //   //we have the total number of molecules, now we check for molfractions
    //     for( i=0; i<n_components; i++ ){
    
    //       if( !the_components[i]->haveMolFraction() ){
    
    //  if( !the_components[i]->haveNMol() ){
    //    //we have a problem
    //    std::cerr << "SimSetup error. Neither molFraction nor "
    //              << " nMol was given in component
    
  }

#ifdef IS_MPI
  strcpy( checkPointMsg, "Have the number of components" );
  MPIcheckPoint();
#endif // is_mpi

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

  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;
  simnfo->n_mol = tot_nmol;

  
#ifdef IS_MPI

  // divide the molecules among processors here.
  
  new mpiSimulation( simnfo );
  
  simnfo->mpiSim->divideLabor( n_components, comp_stamps, components_nmol );

#endif // is_mpi
  

  // create the atom and short range interaction arrays

  Atom::createArrays(simnfo->n_atoms);
  the_atoms = new Atom*[simnfo->n_atoms];
  the_molecules = new Molecule[simnfo->n_mol];


  if( simnfo->n_SRI ){
    the_sris = new SRI*[simnfo->n_SRI];
    the_excludes = new ex_pair[simnfo->n_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();
  }


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

  if( simnfo->n_dipoles ){

    if( !the_globals->haveRRF() ){
      sprintf( painCave.errMsg,
               "SimSetup Error, system has dipoles, but no rRF was set.\n");
      painCave.isFatal = 1;
      simError();
    }
    if( !the_globals->haveDielectric() ){
      sprintf( painCave.errMsg,
               "SimSetup Error, system has dipoles, but no"
               " dielectric was set.\n" );
      painCave.isFatal = 1;
      simError();
    }

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

#ifdef IS_MPI
  strcpy( checkPointMsg, "rRf and dielectric check out" );
  MPIcheckPoint();
#endif // is_mpi
  
  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() ){
      sprintf( painCave.errMsg,
               "SimSetup error, no periodic BoxX size given.\n" );
      painCave.isFatal = 1;
      simError();
    }
    simnfo->box_x = the_globals->getBoxX();

    if( !the_globals->haveBoxY() ){
      sprintf( painCave.errMsg,
               "SimSetup error, no periodic BoxY size given.\n" );
      painCave.isFatal = 1;
      simError();
    }
    simnfo->box_y = the_globals->getBoxY();

    if( !the_globals->haveBoxZ() ){
      sprintf( painCave.errMsg,
               "SimSetup error, no periodic BoxZ size given.\n" );
      painCave.isFatal = 1;
      simError();
    }
    simnfo->box_z = the_globals->getBoxZ();
  }

#ifdef IS_MPI
  strcpy( checkPointMsg, "Box size set up" );
  MPIcheckPoint();
#endif // is_mpi


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

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

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

#ifdef IS_MPI 

  // no init from bass
  
  sprintf( painCave.errMsg,
           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
  painCave.isFatal;
  simError();
  
#else

  initFromBass();

#endif // is_mpi
 
#ifdef IS_MPI
  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
  MPIcheckPoint();
#endif // is_mpi


  //   }
  
#ifdef IS_MPI
  if( worldRank == 0 ){
#endif // is_mpi
    
    if( the_globals->haveFinalConfig() ){
      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
    }
    else{
      strcpy( simnfo->finalName, inFileName );
      char* endTest;
      int nameLength = strlen( simnfo->finalName );
      endTest = &(simnfo->finalName[nameLength - 5]);
      if( !strcmp( endTest, ".bass" ) ){
        strcpy( endTest, ".eor" );
      }
      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" );
      }
    }
    
#ifdef IS_MPI
  }
#endif // is_mpi
  
  // set the status, sample, and themal kick times
  
  if( the_globals->haveSampleTime() ){
    simnfo->sampleTime = the_globals->getSampleTime();
    simnfo->statusTime = simnfo->sampleTime;
    simnfo->thermalTime = simnfo->sampleTime;
  }
  else{
    simnfo->sampleTime = the_globals->getRunTime();
    simnfo->statusTime = simnfo->sampleTime;
    simnfo->thermalTime = simnfo->sampleTime;
  }

  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){
      sprintf( painCave.errMsg,
               "SimSetup error. There has been an error in constructing"
               " the non-complete lattice.\n" );
      painCave.isFatal = 1;
      simError();
    }
  }
  else{
    n_cells = (int)temp3;
    cellx = simnfo->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() ){
      sprintf( painCave.errMsg,
               "SimSetup:initFromBass error.\n"
               "\tComponent %s, atom %s does not have a position specified.\n"
               "\tThe initialization routine is unable to give a start"
               " position.\n",
               comp_stamps[current_comp]->getID(),
               current_atom->getType() );
      painCave.isFatal = 1;
      simError();
    }

    the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
    the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
    the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );

    if( the_atoms[current_atom_ndx]->isDirectional() ){

      dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];

      rotMat[0][0] = 1.0;
      rotMat[0][1] = 0.0;
      rotMat[0][2] = 0.0;

      rotMat[1][0] = 0.0;
      rotMat[1][1] = 1.0;
      rotMat[1][2] = 0.0;

      rotMat[2][0] = 0.0;
      rotMat[2][1] = 0.0;
      rotMat[2][2] = 1.0;

      dAtom->setA( rotMat );
    }

    current_atom_ndx++;
  }

  current_mol++;
  current_comp_mol++;

  if( current_comp_mol >= components_nmol[current_comp] ){

    current_comp_mol = 0;
    current_comp++;
  }
}
Revision:	195
Committed:	Thu Dec 5 18:53:40 2002 UTC (22 years, 10 months ago) by chuckv
File size:	19251 byte(s)
Log Message:	Added divideLabor.