mdtools/interface_implementation/SimSetup.cpp

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

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

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

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

void SimSetup::parseFile( char* fileName ){

  inFileName = fileName;
  set_interface_stamps( stamps, globals );
  yacc_BASS( fileName );
}

void SimSetup::createSim( void ){

  MakeStamps *the_stamps;
  Globals* the_globals;
  int i;

  // get the stamps and globals;
  the_stamps = stamps;
  the_globals = globals;

  // set the easy ones first
  simnfo->target_temp = the_globals->getTargetTemp();
  simnfo->dt = the_globals->getDt();
  simnfo->run_time = the_globals->getRunTime();

  // get the ones we know are there, yet still may need some work.
  n_components = the_globals->getNComponents();
  strcpy( force_field, the_globals->getForceField() );
  strcpy( ensemble, the_globals->getEnsemble() );
  
  if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
  else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
  else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
  else{
    std::cerr<< "SimSetup Error. Unrecognized force field -> " 
             << force_field << "\n";
    exit(8);
  }

  // get the components and calculate the tot_nMol and indvidual n_mol
  the_components = the_globals->getComponents();
  components_nmol = new int[n_components];
  comp_stamps = new MoleculeStamp*[n_components];
  
  if( !the_globals->haveNMol() ){
    // we don't have the total number of molecules, so we assume it is 
    // given in each component

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

      tot_nmol += the_components[i]->getNMol();
      components_nmol[i] = the_components[i]->getNMol();
    }
  }
  else{
    std::cerr << "NOT A SUPPORTED FEATURE\n";
    
//     tot_nmol = the_globals->getNMol();
    
//     //we have the total number of molecules, now we check for molfractions
//     for( i=0; i<n_components; i++ ){
      
//       if( !the_components[i]->haveMolFraction() ){
        
//      if( !the_components[i]->haveNMol() ){
//        //we have a problem
//        std::cerr << "SimSetup error. Neither molFraction nor "
//                  << " nMol was given in component

  }

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

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

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

  

  // caclulate the number of atoms, bonds, bends and torsions

  tot_atoms = 0;
  tot_bonds = 0;
  tot_bends = 0;
  tot_torsions = 0;
  for( i=0; i<n_components; i++ ){
    
    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
  }
  
  tot_SRI = tot_bonds + tot_bends + tot_torsions;
  
  simnfo->n_atoms = tot_atoms;
  simnfo->n_bonds = tot_bonds;
  simnfo->n_bends = tot_bends;
  simnfo->n_torsions = tot_torsions;
  simnfo->n_SRI = tot_SRI;

  // create the atom and short range interaction arrays
  
  the_atoms = new Atom*[tot_atoms];
  the_molecules = new Molecule[tot_nmol];
  
  
  if( tot_SRI ){
    the_sris = new SRI*[tot_SRI];
    the_excludes = new ex_pair[tot_SRI];
  }

  // set the arrays into the SimInfo object

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

  // initialize the arrays
  
  the_ff->setSimInfo( simnfo );
    
  makeAtoms();

  if( tot_bonds ){
    makeBonds();
  }

  if( tot_bends ){
    makeBends();
  }

  if( tot_torsions ){
    makeTorsions();
  }

  //  makeMolecules();

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

  if( simnfo->n_dipoles ){

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

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

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

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

    if( !the_globals->haveBoxZ() ){
      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
      exit(8);
    }
    simnfo->box_z = the_globals->getBoxZ();
  }
    
  if( the_globals->haveInitialConfig() ){
    InitializeFromFile* fileInit;
    fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
    
    fileInit->read_xyz( simnfo ); // default velocities on

    delete fileInit;    
  }
  else{
    initFromBass();
  }

  if( the_globals->haveFinalConfig() ){
    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
  }
  else{
    strcpy( simnfo->finalName, inFileName );
    char* endTest;
    int nameLength = strlen( simnfo->finalName );
    endTest = &(simnfo->finalName[nameLength - 5]);
    if( !strcmp( endTest, ".bass" ) ){
      strcpy( endTest, ".eor" );
    }
    else if( !strcmp( endTest, ".BASS" ) ){
      strcpy( endTest, ".eor" );
    }
    else{
      endTest = &(simnfo->finalName[nameLength - 4]);
      if( !strcmp( endTest, ".bss" ) ){
        strcpy( endTest, ".eor" );
      }
      else if( !strcmp( endTest, ".mdl" ) ){
        strcpy( endTest, ".eor" );
      }
      else{
        strcat( simnfo->finalName, ".eor" );
      }
    }
  }
    
  // make the sample and status out names

  strcpy( simnfo->sampleName, inFileName );
  char* endTest;
  int nameLength = strlen( simnfo->sampleName );
  endTest = &(simnfo->sampleName[nameLength - 5]);
  if( !strcmp( endTest, ".bass" ) ){
    strcpy( endTest, ".dump" );
  }
  else if( !strcmp( endTest, ".BASS" ) ){
    strcpy( endTest, ".dump" );
  }
  else{
    endTest = &(simnfo->sampleName[nameLength - 4]);
    if( !strcmp( endTest, ".bss" ) ){
      strcpy( endTest, ".dump" );
    }
    else if( !strcmp( endTest, ".mdl" ) ){
      strcpy( endTest, ".dump" );
    }
    else{
      strcat( simnfo->sampleName, ".dump" );
    }
  }
  
  strcpy( simnfo->statusName, inFileName );
  nameLength = strlen( simnfo->statusName );
  endTest = &(simnfo->statusName[nameLength - 5]);
  if( !strcmp( endTest, ".bass" ) ){
    strcpy( endTest, ".stat" );
  }
  else if( !strcmp( endTest, ".BASS" ) ){
    strcpy( endTest, ".stat" );
  }
  else{
    endTest = &(simnfo->statusName[nameLength - 4]);
    if( !strcmp( endTest, ".bss" ) ){
      strcpy( endTest, ".stat" );
    }
    else if( !strcmp( endTest, ".mdl" ) ){
      strcpy( endTest, ".stat" );
    }
    else{
      strcat( simnfo->statusName, ".stat" );
    }
  }
  
  // set the status, sample, and themal kick times

  if( the_globals->haveSampleTime() ){
    simnfo->sampleTime = the_globals->getSampleTime(); 
    simnfo->statusTime = simnfo->sampleTime;
    simnfo->thermalTime = simnfo->sampleTime;
  }
  else{
    simnfo->sampleTime = the_globals->getRunTime(); 
    simnfo->statusTime = simnfo->sampleTime;
    simnfo->thermalTime = simnfo->sampleTime;
  }

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

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

  // check for the temperature set flag

  if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
  
    
  // make the longe range forces and the integrator
  
  new AllLong( simnfo );
      
  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
  if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
  if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
}

void SimSetup::makeAtoms( void ){
  
  int i, j, k, index;
  double ux, uy, uz, uSqr, u;
  AtomStamp* current_atom;
  DirectionalAtom* dAtom;

  index = 0;
  for( i=0; i<n_components; i++ ){
    
    for( j=0; j<components_nmol[i]; j++ ){
      
      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
        
        current_atom = comp_stamps[i]->getAtom( k );
        if( current_atom->haveOrientation() ){ 

          dAtom = new DirectionalAtom;
          simnfo->n_oriented++;
          the_atoms[index] = dAtom;
      
          ux = current_atom->getOrntX();
          uy = current_atom->getOrntY();
          uz = current_atom->getOrntZ();
          
          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
          
          u = sqrt( uSqr );
          ux = ux / u;
          uy = uy / u;
          uz = uz / u;
          
          dAtom->setSUx( ux );
          dAtom->setSUy( uy );
          dAtom->setSUz( uz );
        }
        else{
          the_atoms[index] = new GeneralAtom;
        }
        the_atoms[index]->setType( current_atom->getType() );
        the_atoms[index]->setIndex( index );
        
        // increment the index and repeat;
        index++;
      }
    }
  }
  
  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::makeMolecules( void ){

  int i,j,k;
  
  for( i=0; i<n_components; i++ ){
    
    for( j=0; j<components_nmol[i]; j++ ){
      
      

}

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