utils/sysbuilder/latticeBilayer.cpp

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "simError.h"
#include "SimInfo.hpp"
#include "ReadWrite.hpp"
#include "SimSetup.hpp"

#include "MoLocator.hpp"
#include "latticeBuilder.hpp"


#define VERSION_MAJOR 0
#define VERSION_MINOR 1

char *programName; /*the name of the program */
void usage(void);
int buildLatticeBilayer( double aLat, 
                         double bLat, 
                         double leafSpacing,
                         char* waterName,
                         char* lipidName);
using namespace std;
SimInfo* mainInfo;

int main(int argC,char* argV[]){
  
  int i,j; // loop counters

  char* outPrefix; // the output prefix

  char* conversionCheck;
  bool conversionError;
  bool optionError;

  char currentFlag; // used in parsing the flags
  bool done = false; // multipurpose boolean
  bool havePrefix; // boolean for the output prefix  

  char* lipidName;
  char* waterName;
  bool haveWaterName, haveLipidName;

  bool haveSpacing, haveAlat, haveBlat;
  double aLat, bLat, leafSpacing;
  
  char* inName;
  
  SimSetup* simInit;
  
  // first things first, all of the initializations

  fflush(stdout);
  srand48( 1337 );          // the random number generator.
  initSimError();           // the error handler

  outPrefix = NULL;
  inName = NULL;
  
  conversionError = false;
  optionError = false;
 
  havePrefix = false;
  haveAlat = false;
  haveBlat = false;
  haveSpacing = false;

  programName = argV[0]; /*save the program name in case we need it*/

  for( i = 1; i < argC; i++){
    
    if(argV[i][0] =='-'){

      // parse the option
      
      if(argV[i][1] == '-' ){

        // parse long word options
        
        if( !strcmp( argV[i], "--version") ){
          
          printf("\n"
                 "staticProps version %d.%d\n"
                 "\n",
                 VERSION_MAJOR, VERSION_MINOR );
          exit(0);
          
        }

        else if( !strcmp( argV[i], "--help") ){
          
          usage();
          exit(0);
        }
        
        // anything else is an error

        else{
          fprintf( stderr, 
                   "Invalid option \"%s\"\n", argV[i] );
          usage();
          exit(0);
        }
      }
      
      else{
        
        // parse single character options
        
        done =0;
        j = 1;
        currentFlag = argV[i][j];
        while( (currentFlag != '\0') && (!done) ){
          
          switch(currentFlag){

          case 'o':
            // -o <prefix> => the output prefix.

            j++;
            currentFlag = argV[i][j];
              
            if( currentFlag != '\0' ) optionError = true;
            
            if( optionError ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "The -o flag should end an option sequence.\n"
                       "   example: -r <outname> *NOT* -or <outname>\n" );
              usage();
              painCave.isFatal = 1;
              simError();
            }
              
            i++;
            if( i>=argC ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "not enough arguments for -o\n");
              usage();
              painCave.isFatal = 1;
              simError();
            }     
            
            outPrefix = argV[i];
            if( outPrefix[0] == '-' ) optionError = true;
                
            if( optionError ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "\"%s\" is not a valid out prefix/name.\n"
                       "Out prefix/name should not begin with a dash.\n",
                       outPrefix );
              usage();
              painCave.isFatal = 1;
              simError();
            }
            
            havePrefix = true;
            done = true;
            break;

          case 'l':
            // -l <lipidName> => the lipid name.

            j++;
            currentFlag = argV[i][j];
              
            if( currentFlag != '\0' ) optionError = true;
            
            if( optionError ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "The -l flag should end an option sequence.\n"
                       "   example: -rl <lipidName> *NOT* -lr <lipidName>\n" );
              usage();
              painCave.isFatal = 1;
              simError();
            }
              
            i++;
            if( i>=argC ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "not enough arguments for -l\n");
              usage();
              painCave.isFatal = 1;
              simError();
            }     
            
            lipidName = argV[i];
            if( lipidName[0] == '-' ) optionError = true;
                
            if( optionError ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "\"%s\" is not a valid lipidName.\n"
                       "lipidName should not begin with a dash.\n",
                       lipidName );
              usage();
              painCave.isFatal = 1;
              simError();
            }
            
            haveLipidName = true;
            done = true;
            break;

          case 'w':
            // -w <waterName> => the water name.

            j++;
            currentFlag = argV[i][j];
              
            if( currentFlag != '\0' ) optionError = true;
            
            if( optionError ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "The -w flag should end an option sequence.\n"
                       "   example: -rw <waterName> *NOT* -lw <waterName>\n" );
              usage();
              painCave.isFatal = 1;
              simError();
            }
              
            i++;
            if( i>=argC ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "not enough arguments for -w\n");
              usage();
              painCave.isFatal = 1;
              simError();
            }     
            
            waterName = argV[i];
            if( waterName[0] == '-' ) optionError = true;
                
            if( optionError ){
              sprintf( painCave.errMsg, 
                       "\n"
                       "\"%s\" is not a valid waterName.\n"
                       "waterName should not begin with a dash.\n",
                       waterName );
              usage();
              painCave.isFatal = 1;
              simError();
            }
            
            haveWaterName = true;
            done = true;
            break;


          case 'h':
            // -h <double>    set <double> to the hex lattice spacing

            haveAlat = true;
            haveBlat = true;
            j++;
            currentFlag = argV[i][j];
              
            if( currentFlag != '\0' ) optionError = true;
            
            if( optionError ){
              sprintf( painCave.errMsg,
                       "\n"
                       "The -h flag should end an option sequence.\n"
                       "   example: -sh <double> *NOT* -hs <double>\n" );
              usage();
              painCave.isFatal = 1;
              simError();
            }

            i++;
            if( i>=argC ){
              sprintf( painCave.errMsg,
                       "\n"
                       "not enough arguments for -h\n");
              usage();
              painCave.isFatal = 1;
              simError();
            }     


            bLat = atof( argV[i] );
//          bLat = strtod( argV[i], &conversionCheck);
//          if( conversionCheck == argV[i] ) conversionError = true;
//          if( *conversionCheck != '\0' ) conversionError = true;
            
            if( conversionError ){
              sprintf( painCave.errMsg,
                       "Error converting \"%s\" to a double for \"h\".\n",
                       argV[i] );
              usage();
              painCave.isFatal = 1;
              simError();
            }
            
            aLat = sqrt( 3.0 ) * bLat;
            
            done = true;

            break;

          case 'a':
            // -a <double>    set <double> to the aLat

            haveAlat = true;
            j++;
            currentFlag = argV[i][j];
              
            if( currentFlag != '\0' ) optionError = true;
            
            if( optionError ){
              sprintf( painCave.errMsg,
                       "\n"
                       "The -a flag should end an option sequence.\n"
                       "   example: -sa <double> *NOT* -as <double>\n" );
              usage();
              painCave.isFatal = 1;
              simError();
            }

            i++;
            if( i>=argC ){
              sprintf( painCave.errMsg,
                       "\n"
                       "not enough arguments for -a\n");
              usage();
              painCave.isFatal = 1;
              simError();
            }     

            aLat = atof( argV[i] );
//          aLat = strtod( argV[i], &conversionCheck);
//          if( conversionCheck == argV[i] ) conversionError = true;
//          if( *conversionCheck != '\0' ) conversionError = true;
            
            if( conversionError ){
              sprintf( painCave.errMsg,
                       "Error converting \"%s\" to a double for \"a\".\n",
                       argV[i] );
              usage();
              painCave.isFatal = 1;
              simError();
            }
            
            done = true;

            break;

          case 'b':
            // -b <double>    set <double> to the bLat

            haveBlat = true;
            j++;
            currentFlag = argV[i][j];
              
            if( currentFlag != '\0' ) optionError = true;
            
            if( optionError ){
              sprintf( painCave.errMsg,
                       "\n"
                       "The -b flag should end an option sequence.\n"
                       "   example: -sb <double> *NOT* -bs <double>\n" );
              usage();
              painCave.isFatal = 1;
              simError();
            }

            i++;
            if( i>=argC ){
              sprintf( painCave.errMsg,
                       "\n"
                       "not enough arguments for -b\n");
              usage();
              painCave.isFatal = 1;
              simError();
            }     

            bLat = atof( argV[i] );

//          bLat = strtod( argV[i], &conversionCheck);
//          if( conversionCheck == argV[i] ) conversionError = true;
//          if( *conversionCheck != '\0' ) conversionError = true;
            
            if( conversionError ){
              sprintf( painCave.errMsg,
                       "Error converting \"%s\" to a double for \"a\".\n",
                       argV[i] );
              usage();
              painCave.isFatal = 1;
              simError();
            }
            
            done = true;

            break;

          case 's':
            // -s <double>    set <double> to the leafSpacing

            haveSpacing = true;
            j++;
            currentFlag = argV[i][j];
              
            if( currentFlag != '\0' ) optionError = true;
            
            if( optionError ){
              sprintf( painCave.errMsg,
                       "\n"
                       "The -s flag should end an option sequence.\n"
                       "   example: -rs <double> *NOT* -sr <double>\n" );
              usage();
              painCave.isFatal = 1;
              simError();
            }

            i++;
            if( i>=argC ){
              sprintf( painCave.errMsg,
                       "\n"
                       "not enough arguments for -s\n");
              usage();
              painCave.isFatal = 1;
              simError();
            }     

            leafSpacing = atof( argV[i] );

//          leafSpacing = strtod( argV[i], &conversionCheck);
//          if( conversionCheck == argV[i] ) conversionError = true;
//          if( *conversionCheck != '\0' ) conversionError = true;
            
            if( conversionError ){
              sprintf( painCave.errMsg,
                       "Error converting \"%s\" to a double for \"s\".\n",
                       argV[i] );
              usage();
              painCave.isFatal = 1;
              simError();
            }
            
            done = true;

            break;

          default:

            sprintf(painCave.errMsg, 
                    "\n"
                    "Bad option \"-%c\"\n", currentFlag);
            usage();
            painCave.isFatal = 1;
            simError();
          }
          j++;
          currentFlag = argV[i][j];
        }
      }
    }

    else{
      
      if( inName != NULL ){
        sprintf( painCave.errMsg,
                 "Error at \"%s\", program does not currently support\n"
                 "more than one input bass file.\n"
                 "\n",
                 argV[i]);
        usage();
        painCave.isFatal = 1;
        simError();
      }
      
      inName = argV[i];
      
    }
  }  
  
  if( inName == NULL ){
    sprintf( painCave.errMsg,
             "Error, bass file is needed to run.\n" );
    usage();
    painCave.isFatal = 1;
    simError();
  }

  // if no output prefix is given default to "donkey".

  if( !havePrefix ){
    outPrefix = strdup( "donkey" );
  }


  if( !haveWaterName ){
    sprintf( painCave.errMsg,
             "Error, the water name is needed to run.\n"
             );
    usage();
    painCave.isFatal = 1;
    simError();
  }  
  
  if( !haveLipidName ){
    sprintf( painCave.errMsg,
             "Error, the lipid name is needed to run.\n"
             );
    usage();
    painCave.isFatal = 1;
    simError();
  }
  
  if( !haveAlat ){
    sprintf( painCave.errMsg,
             "Error, the hexagonal lattice parameter, a, is needed to run.\n"
             );
    usage();
    painCave.isFatal = 1;
    simError();
  }
  
  if( !haveBlat ){
    sprintf( painCave.errMsg,
             "Error, the hexagonal lattice parameter, b, is needed to run.\n"
             );
    usage();
    painCave.isFatal = 1;
    simError();
  }


  if( !haveSpacing ){
    sprintf( painCave.errMsg,
             "Error, the leaf spaceing is needed to run.\n"
             );
    usage();
    painCave.isFatal = 1;
    simError();
  }

  mainInfo = new SimInfo();
  simInit = new SimSetup();
  simInit->setSimInfo( mainInfo );
  simInit->suspendInit();
  simInit->parseFile( inName );
   simInit->createSim();

  delete simInit;

  sprintf( mainInfo->statusName, "%s.stat", outPrefix );
  sprintf( mainInfo->sampleName, "%s.dump", outPrefix );
  sprintf( mainInfo->finalName, "%s.init", outPrefix );

  buildLatticeBilayer( aLat, bLat, leafSpacing, waterName, lipidName );

  return 0;  
}

int buildLatticeBilayer(double aLat, 
                        double bLat, 
                        double leafSpacing,
                        char* waterName,
                        char* lipidName){

  typedef struct{
    double rot[3][3];
    double pos[3];
  } coord;

  const double waterRho = 0.0334; // number density per cubic angstrom
  const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters
  
  double waterCell[3];

  double *posX, *posY, *posZ;
  double pos[3], posA[3], posB[3];

  const double waterFudge = 5.0;

  int i,j,k,l;
  int nAtoms, atomIndex, molIndex, molID;
  int* molSeq;
  int* molMap;
  int* molStart;
  int testTot, done;
  int nCells, nCellsX, nCellsY, nCellsZ; 
  int nx, ny;
  double boxX, boxY, boxZ;
  double theta, phi, psi;
  int which;
  int targetWaters;
  
  Atom** atoms;
  SimInfo* simnfo;
  SimInfo* testInfo;
  coord testSite;
  SimState* theConfig;
  DumpWriter* writer;

  MoleculeStamp* lipidStamp;
  MoleculeStamp* waterStamp;  
  MoLocator *lipidLocate;
  MoLocator *waterLocate;
  int foundLipid, foundWater;
  int nLipids, lipidNatoms, nWaters, waterNatoms;
  int targetNlipids, targetNwaters;
  double targetWaterLipidRatio;
  double maxZ, minZ, zHeight;
  double maxY, minY;
  double maxX, minX;
   
  molStart = NULL;

  // create the simInfo objects

  simnfo = new SimInfo;

  // set the the lipidStamp

  foundLipid = 0;
  foundWater = 0;
 
  for(i=0; i<mainInfo->nComponents; i++){
     
    if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){
          
      foundLipid = 1;
      lipidStamp = mainInfo->compStamps[i];
      targetNlipids = mainInfo->componentsNmol[i];
      lipidNatoms = lipidStamp->getNAtoms();
    }
    if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){
          
      foundWater = 1;
          
      waterStamp = mainInfo->compStamps[i];
      targetNwaters = mainInfo->componentsNmol[i];
      waterNatoms = waterStamp->getNAtoms();
    }
  }
  if( !foundLipid ){
    sprintf(painCave.errMsg,
            "latticeBilayer error: Could not find lipid \"%s\" in the bass file.\n",
            lipidName );
    painCave.isFatal = 1;
    simError();
  }
  if( !foundWater ){
    sprintf(painCave.errMsg,
            "latticeBilayer error: Could not find solvent \"%s\" in the bass file.\n",
            waterName );
    painCave.isFatal = 1;
    simError();
  }
  
  //create the Molocator arrays
  
  lipidLocate = new MoLocator( lipidStamp );
  waterLocate = new MoLocator( waterStamp );

  // gather info about the lipid

  testInfo = new SimInfo();
  testInfo->n_atoms = lipidNatoms;
  theConfig = testInfo->getConfiguration();
  theConfig->createArrays( lipidNatoms );
  testInfo->atoms = new Atom*[lipidNatoms];
  atoms = testInfo->atoms;

  testSite.pos[0] = 0.0;
  testSite.pos[1] = 0.0;
  testSite.pos[2] = 0.0;
  
  theta = 0.0;
  phi = 0.0;
  psi = 0.0;

  getEulerRot(theta, phi, psi, testSite.rot );
  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig );

  minZ = 0.0;
  maxZ = 0.0;
  double myPos[3];
  for(i=0;i<lipidNatoms;i++){
    atoms[i]->getPos( myPos );
    minZ = (minZ > myPos[2]) ? myPos[2] : minZ;
    maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ;
  }
  zHeight = maxZ - minZ;
    
  nCells = (int) sqrt( (double)targetNlipids * bLat / (4.0 * aLat) );

  nx = nCells;
  ny = (int) ((double)nCells  * aLat / bLat);

  boxX = nx * aLat;
  boxY = ny * bLat;  

  nLipids = 4 * nx * ny;
  coord* lipidSites = new coord[nLipids];

  phi = 0.0;
  psi = 0.0;

  which = 0;

  for (i = 0; i < nx; i++) {
    for (j = 0; j < ny; j++ ) {
      for (k = 0; k < 2; k++) {
        
        lipidSites[which].pos[0] = (double)i * aLat;
        lipidSites[which].pos[1] = (double)j * bLat;
        lipidSites[which].pos[2] = (2.0* (double)k - 1.0) * 
          ((leafSpacing / 2.0) - maxZ);
        
        theta = (1.0 - (double)k) * M_PI;
        
        getEulerRot( theta, phi, psi, lipidSites[which].rot );

        which++;

        lipidSites[which].pos[0] = aLat * ((double)i + 0.5);
        lipidSites[which].pos[1] = bLat * ((double)j + 0.5);
        lipidSites[which].pos[2] = (2.0* (double)k - 1.0) * 
          ((leafSpacing / 2.0) - maxZ);
        
        theta = (1.0 - (double)k) * M_PI;
        
        getEulerRot( theta, phi, psi, lipidSites[which].rot );
        
        which++;
      }
    }
  }

  targetWaterLipidRatio = (double)targetNwaters / (double)targetNlipids;
 
  targetWaters = targetWaterLipidRatio * nLipids;

  // guess the size of the water box
  
  nCellsX = (int)ceil(boxX / pow(waterVol, ( 1.0 / 3.0 )) );
  nCellsY = (int)ceil(boxY / pow(waterVol, ( 1.0 / 3.0 )) );

  done = 0;
  nCellsZ = 0;
  while( !done ){
        
    nCellsZ++;
    testTot = 4 * nCellsX * nCellsY * nCellsZ;
        
    if( testTot >= targetWaters ) done = 1;
  }
  
  nWaters = nCellsX * nCellsY * nCellsZ * 4;
  
  // calc current system size;

  nAtoms = 0;
  molIndex = 0;
  if(molStart != NULL ) delete[] molStart;
  molStart = new int[nLipids];  
  
  for(j=0; j<nLipids; j++){
    molStart[molIndex] = nAtoms;
    molIndex++;
    nAtoms += lipidNatoms;
  }

  testInfo->n_atoms = nAtoms;
  theConfig = testInfo->getConfiguration();
  theConfig->destroyArrays();
  theConfig->createArrays( nAtoms );
  testInfo->atoms = new Atom*[nAtoms];
  atoms = testInfo->atoms;

  molIndex = 0;
  for(i=0; i<nLipids; i++ ){
    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
                           molStart[molIndex], theConfig );
    molIndex++;
  }

  atoms[0]->getPos( myPos );

  maxX = myPos[0];
  minX = myPos[0];

  maxY = myPos[1];
  minY = myPos[1];

  maxZ = myPos[2];
  minZ = myPos[2];

  for(i=0;i<nAtoms;i++){
    atoms[i]->getPos( myPos );
    minX = (minX > myPos[0]) ? myPos[0] : minX;
    maxX = (maxX < myPos[0]) ? myPos[0] : maxX;

    minY = (minY > myPos[1]) ? myPos[1] : minY;
    maxY = (maxY < myPos[1]) ? myPos[1] : maxY;

    minZ = (minZ > myPos[2]) ? myPos[2] : minZ;
    maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ;
  }

  boxX = (maxX - minX)+2.0;
  boxY = (maxY - minY)+2.0;
  boxZ = (maxZ - minZ)+2.0;

  double centerX, centerY, centerZ;
  
  centerX = ((maxX - minX) / 2.0) + minX;
  centerY = ((maxY - minY) / 2.0) + minY;
  centerZ = ((maxZ - minZ) / 2.0) + minZ;

  // set up water coordinates

  coord* waterSites = new coord[nWaters];

  waterCell[0] = boxX / nCellsX;
  waterCell[1] = boxY / nCellsY;
  waterCell[2] = 4.0 / (waterRho * waterCell[0] * waterCell[1]);

  Lattice *myORTHO;
  myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell);
  myORTHO->setStartZ( maxZ + waterFudge);

  // create an fcc lattice in the water box.
  
  which = 0;
  for( i=0; i < nCellsX; i++ ){
    for( j=0; j < nCellsY; j++ ){
      for( k=0; k < nCellsZ; k++ ){

        myORTHO->getLatticePoints(&posX, &posY, &posZ, i, j, k);
        for(l=0; l<4; l++){
          waterSites[which].pos[0] = posX[l];
          waterSites[which].pos[1] = posY[l];
          waterSites[which].pos[2] = posZ[l];
          getRandomRot( waterSites[which].rot );
          which++;
        }
      }
    }
  }  

  // calc the system sizes
  
  nAtoms = 0;
  molIndex = 0;
  if(molStart != NULL ) delete[] molStart;
  molStart = new int[nLipids + nWaters];  
  
  for(j=0; j<nLipids; j++){
    molStart[molIndex] = nAtoms;
    molIndex++;
    nAtoms += lipidNatoms;
  }

  for(j=0; j<nWaters; j++){
    molStart[molIndex] = nAtoms;
    molIndex++;
    nAtoms += waterNatoms;
  }
  
  testInfo->n_atoms = nAtoms;
  theConfig = testInfo->getConfiguration();
  theConfig->destroyArrays();
  theConfig->createArrays( nAtoms );
  testInfo->atoms = new Atom*[nAtoms];
  atoms = testInfo->atoms;

  molIndex = 0;
  for(i=0; i<nLipids; i++ ){
    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
                           molStart[molIndex], theConfig );
    molIndex++;
  }

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

    getRandomRot( waterSites[i].rot );
    waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
                           molStart[molIndex], theConfig );
    molIndex++;
  }

  atoms[0]->getPos( myPos );
  maxZ = myPos[2];
  minZ = myPos[2];
  for(i=0;i<nAtoms;i++){
    atoms[i]->getPos( myPos );
    minZ = (minZ > myPos[2]) ? myPos[2] : minZ;
    maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ;
  }
  boxZ = (maxZ - (minZ - waterFudge / 2.0));
  
  // create the real Atom arrays 

  theConfig = simnfo->getConfiguration();
  theConfig->createArrays( nAtoms );
  simnfo->atoms = new Atom*[nAtoms];
  simnfo->n_atoms = nAtoms;
  atoms = simnfo->atoms;

  // wrap back to <0,0,0> as center

  double Hmat[3][3];

  Hmat[0][0] = boxX;
  Hmat[0][1] = 0.0;
  Hmat[0][2] = 0.0;

  Hmat[1][0] = 0.0;
  Hmat[1][1] = boxY;
  Hmat[1][2] = 0.0;

  Hmat[2][0] = 0.0;
  Hmat[2][1] = 0.0;
  Hmat[2][2] = boxZ;
  
  mainInfo->setBoxM( Hmat );
  simnfo->setBoxM( Hmat );

  for(j=0;j<nLipids;j++){

    lipidSites[j].pos[0] -= centerX;
    lipidSites[j].pos[1] -= centerY;
    lipidSites[j].pos[2] -= centerZ;

    simnfo->wrapVector( lipidSites[j].pos );
  }

  for(j=0;j<nWaters;j++){

    waterSites[j].pos[0] -= centerX;
    waterSites[j].pos[1] -= centerY;
    waterSites[j].pos[2] -= centerZ;

    simnfo->wrapVector( waterSites[j].pos );
  }

  // initialize lipid positions

  molIndex = 0;
  for(i=0; i<nLipids; i++ ){
    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
                           molStart[molIndex], theConfig );
    molIndex++;
  }

  // initialize the water positions

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

    getRandomRot( waterSites[i].rot );
    waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
                           molStart[molIndex], theConfig );
    molIndex++;
  }

  strcpy( simnfo->sampleName, mainInfo->sampleName );
  strcpy( simnfo->finalName, mainInfo->finalName );

  // set up the writer and write out
  
  writer = new DumpWriter( simnfo );
  writer->writeFinal( 0.0 );

  std::cout << "\n"
            << "----------------------------------------\n"
            << "\n"
            << " final nLipids = " << nLipids << "\n"
            << " final nWaters = " << nWaters << "\n"
            << "\n";
       
  return 1;
}


/***************************************************************************
 * prints out the usage for the command line arguments, then exits.
 ***************************************************************************/

void usage(){
  (void)fprintf(stdout, 
                "\n"
                "The proper usage is: %s [options] <input_file>\n"
                "\n"
                "Options:\n"
                "\n"
                "   short:\n"
                "   ------\n"
                "   -o <name>       The output prefix\n"
                "   -l <lipidName>  The name of the lipid molecule specified in the BASS file.\n"
                "   -w <waterName>  The name of the water molecule specified in the BASS file.\n"
                "   -s <double>     The leaf spaceing of the bilayer.\n"
                "   -a <double>     The hexagonal lattice constant, a, for the bilayer leaf.\n"
                "   -b <double>     The hexagonal lattice constant, b, for the bilayer leaf.\n"
                "   -h <double>     Use to set a and b for a normal hex lattice with spacing <double>\n"

                "\n"
                "   long:\n"
                "   -----\n"

                "   --version       displays the version number\n"
                "   --help          displays this help message.\n"
                "\n"
                "\n",
                programName);
}
Revision:	855
Committed:	Thu Nov 6 22:01:37 2003 UTC (22 years ago) by mmeineke
File size:	23053 byte(s)
Log Message:	added the following parameters to BASS: * useInitialExtendedSystemState * orthoBoxTolerance * useIntiTime => useInitialTime