--- trunk/OOPSE/utils/sysbuilder/randomBilayer.cpp	2003/12/10 16:52:46	876
+++ trunk/OOPSE/utils/sysbuilder/randomBilayer.cpp	2003/12/19 15:12:38	886
@@ -1,23 +1,25 @@
-
 #include <iostream>
 #include <vector>
 #include <algorithm>
 
-#include <cstdlib>
-#include <cstring>
-#include <cmath>
+#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 "sysBuild.hpp"
-#include "bilayerSys.hpp"
-
 #include "latticeBuilder.hpp"
 
+#define RAND_SEED 31337 // \/\/007!
+
+#define VERSION_MAJOR 0
+#define VERSION_MINOR 1
+
 class SortCond{
   
 public:
@@ -32,42 +34,396 @@ void buildMap( double &x, double &y, double &z,
 void buildMap( double &x, double &y, double &z,
 	       double boxX, double boxY, double boxZ );
 
-int buildRandomBilayer( void );
+int buildRandomBilayer( double waterCell,
+			double waterBuffer,
+			double lipidBuffer,
+			char* waterName,
+			char* lipidName );
 
-int buildBilayer( int isRandom ){
+
+char *programName; /*the name of the program */
+void usage(void);
+using namespace std;
+SimInfo* mainInfo;
+
+int main(int argC,char* argV[]){
   
-  if( isRandom ){
-    return buildRandomBilayer();
-  }
-  else{
+  int i,j; // loop counters
 
-    std::cerr << "unsupported feature\n";
-    return 0;
+  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;
+
+  double waterLattice, waterBuffer, lipidBuffer;
+  
+  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;
+  
+  waterBuffer = 5.0;
+  lipidBuffer = 6.0;
+  waterLattice = 4.929;
+
+  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"
+		 "randomBilayer version %d.%d\n"
+		 "\n",
+		 VERSION_MAJOR, VERSION_MINOR );
+	  exit(0);
+	  
+	}
+
+	else if( !strcmp( argV[i], "--help") ){
+	  
+	  usage();
+	  exit(0);
+	}
+
+	else if (!strcmp( argV[i], "--lipidBuffer" )){
+	  
+	  i++;
+	  if( i>=argC ){
+	    sprintf( painCave.errMsg,
+		     "\n"
+		     "not enough arguments for -lipidBuffer\n");
+	    usage();
+	    painCave.isFatal = 1;
+	    simError();
+	  }	  
+	  
+	  lipidBuffer = atof( argV[i] );
+	}
+
+	else if (!strcmp( argV[i], "--waterBuffer" )){
+	  
+	  i++;
+	  if( i>=argC ){
+	    sprintf( painCave.errMsg,
+		     "\n"
+		     "not enough arguments for --waterBuffer\n");
+	    usage();
+	    painCave.isFatal = 1;
+	    simError();
+	  }	  
+	  
+	  waterBuffer = atof( argV[i] );
+	}
+
+	else if (!strcmp( argV[i], "--waterLattice" )){
+	  
+	  i++;
+	  if( i>=argC ){
+	    sprintf( painCave.errMsg,
+		     "\n"
+		     "not enough arguments for -waterLattice\n");
+	    usage();
+	    painCave.isFatal = 1;
+	    simError();
+	  }	  
+	  
+	  waterLattice = atof( argV[i] );
+	}
+	  
+	  
+	
+	// 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;
+
+	  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".
 
-int buildRandomBilayer( void ){
+  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();
+  }
+  
+  
+  // create and initialize the info object
+
+  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 );
+
+  buildRandomBilayer( waterLattice, waterBuffer, lipidBuffer, 
+		      waterName, lipidName );
+
+  return 0;  
+}
+
+int buildRandomBilayer( double waterCell,
+			double water_padding,
+			double lipid_spaceing,
+			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
-  const double waterCell = 4.929; // fcc unit cell length
-
   Lattice myFCC( FCC_LATTICE_TYPE, waterCell );
   double *posX, *posY, *posZ;
   double pos[3], posA[3], posB[3];
 
-  const double water_padding = 6.0;
-  const double lipid_spaceing = 8.0;
-
-
   int i,j,k, l, m;
   int nAtoms, atomIndex, molIndex, molID;
   int* molSeq;
@@ -93,45 +449,48 @@ int buildRandomBilayer( void ){
   int foundLipid, foundWater;
   int nLipids, lipidNatoms, nWaters, waterNatoms;
   double testBox, maxLength;
-  
-  srand48( RAND_SEED );
 
+  double waterRho, waterVol;
 
-  // create the simInfo objects
-
-  simnfo = new SimInfo[3];
+  
+  srand48( RAND_SEED );
 
+  // calculate the water parameters
+  
+  waterVol = waterCell * waterCell * waterCell;
+  waterRho = 4.0 / waterVol;
 
   // set the the lipidStamp
 
   foundLipid = 0;
   foundWater = 0;
-  for(i=0; i<bsInfo.nComponents; i++){
-    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
+  for(i=0; i<mainInfo->nComponents; i++){
+    
+    if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){
 	  
       foundLipid = 1;
-      lipidStamp = bsInfo.compStamps[i];
-      nLipids = bsInfo.componentsNmol[i];
+      lipidStamp = mainInfo->compStamps[i];
+      nLipids = mainInfo->componentsNmol[i];
     }
-    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
+    if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){
 	  
       foundWater = 1;
 	  
-      waterStamp = bsInfo.compStamps[i];
-      nWaters = bsInfo.componentsNmol[i];
+      waterStamp = mainInfo->compStamps[i];
+      nWaters = mainInfo->componentsNmol[i];
     }
   }
   if( !foundLipid ){
     sprintf(painCave.errMsg,
-	    "Could not find lipid \"%s\" in the bass file.\n",
-	    bsInfo.lipidName );
+	    "randomBilayer error: Could not find lipid \"%s\" in the bass file.\n",
+	    lipidName );
     painCave.isFatal = 1;
     simError();
   }
   if( !foundWater ){
     sprintf(painCave.errMsg,
-	    "Could not find solvent \"%s\" in the bass file.\n",
-	    bsInfo.waterName );
+	    "randomBilayer error: Could not find solvent \"%s\" in the bass file.\n",
+	    waterName );
     painCave.isFatal = 1;
     simError();
   }
@@ -145,20 +504,20 @@ int buildRandomBilayer( void ){
   waterLocate = new MoLocator( waterStamp );
   waterNatoms = waterStamp->getNAtoms();
   
-  nAtoms = lipidNatoms;
+  // create a test box
 
-  simnfo[0].n_atoms = nAtoms;
-  simnfo[0].atoms=new Atom*[nAtoms];
- 
-  theConfig = simnfo[0].getConfiguration();
-  theConfig->createArrays( simnfo[0].n_atoms );
+  SimInfo* testInfo;
 
-  atoms=simnfo[0].atoms;
-	
+  testInfo = new SimInfo();
+  testInfo->n_atoms = lipidNatoms;
+  theConfig = testInfo->getConfiguration();
+  theConfig->createArrays( lipidNatoms );
+  testInfo->atoms = new Atom*[lipidNatoms];
+  atoms = testInfo->atoms;
 
   // create the test box for initial water displacement
 
-  testBox = maxLength + waterCell * 4.0; // pad with 4 cells
+  testBox = maxLength + waterCell * 10.0; // pad with 4 cells
   nCells = (int)( testBox / waterCell + 1.0 );
   int testWaters = 4 * nCells * nCells * nCells;
   
@@ -251,14 +610,18 @@ int buildRandomBilayer( void ){
   const double boxTargetX = 66.22752;
   const double boxTargetY = 60.53088;
 
-  nCellsX = (int)ceil(boxTargetX / waterCell);
-  nCellsY = (int)ceil(boxTargetY / waterCell);
+//   nCellsX = (int)ceil(boxTargetX / waterCell);
+//   nCellsY = (int)ceil(boxTargetY / waterCell);
   
   int testTot;
   int done = 0;
+  nCellsX = 0;
+  nCellsY = 0;
   nCellsZ = 0;
   while( !done ){
 	
+    nCellsX++;
+    nCellsY++;
     nCellsZ++;
     testTot = 4 * nCellsX * nCellsY * nCellsZ;
 	
@@ -304,17 +667,16 @@ int buildRandomBilayer( void ){
   
   int reject;
   int testDX, acceptedDX;
+  SimInfo* testInfo2;
 
   nAtoms = nLipids * lipidNatoms;
+  testInfo2 = new SimInfo();
+  testInfo2->n_atoms = nAtoms;
+  theConfig = testInfo2->getConfiguration();
+  theConfig->createArrays( nAtoms );
+  testInfo2->atoms = new Atom*[nAtoms];
+  atoms = testInfo2->atoms;
 
-  simnfo[1].n_atoms = nAtoms;
-  simnfo[1].atoms=new Atom*[nAtoms];
- 
-  theConfig = simnfo[1].getConfiguration();
-  theConfig->createArrays( simnfo[1].n_atoms );
-
-  atoms=simnfo[1].atoms;
-
   rCutSqr = lipid_spaceing * lipid_spaceing;
 
   for(i=0; i<nLipids; i++ ){
@@ -367,35 +729,36 @@ int buildRandomBilayer( void ){
       else{
 	done = 1;
 	std::cout << (i+1) << " has been accepted\n";
+	std::cout.flush();
       }
     }
   }
 	
 
-  // zSort of the lipid positions
+//   // zSort of the lipid positions
 
 
-  vector< pair<int,double> >zSortArray;
-  for(i=0;i<nLipids;i++) 
-    zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) );
+//   vector< pair<int,double> >zSortArray;
+//   for(i=0;i<nLipids;i++) 
+//     zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) );
 
-  sort(zSortArray.begin(),zSortArray.end(),SortCond());
+//   sort(zSortArray.begin(),zSortArray.end(),SortCond());
 
-  ofstream outFile( "./zipper.bass", ios::app);
+//   ofstream outFile( "./zipper.bass", ios::app);
 
-  for(i=0; i<nLipids; i++){
-    outFile << "zConstraint[" << i << "]{\n"
-	    << "  molIndex = " << zSortArray[i].first <<  ";\n"
-	    << "  zPos = ";
+//   for(i=0; i<nLipids; i++){
+//     outFile << "zConstraint[" << i << "]{\n"
+// 	    << "  molIndex = " << zSortArray[i].first <<  ";\n"
+// 	    << "  zPos = ";
 
-    if(i<32) outFile << "60.0;\n";
-    else outFile << "100.0;\n";
+//     if(i<32) outFile << "60.0;\n";
+//     else outFile << "100.0;\n";
 
-    outFile << "  kRatio = 0.5;\n"
-	    << "}\n";
-  }
+//     outFile << "  kRatio = 0.5;\n"
+// 	    << "}\n";
+//   }
   
-  outFile.close();
+//   outFile.close();
 
 
   // cut out the waters that overlap with the lipids.
@@ -467,8 +830,8 @@ int buildRandomBilayer( void ){
 
   // clean up our messes before building the final system.
 
-  simnfo[0].getConfiguration()->destroyArrays();
-  simnfo[1].getConfiguration()->destroyArrays();
+  testInfo->getConfiguration()->destroyArrays();
+  testInfo2->getConfiguration()->destroyArrays();
 
   // create the real Atom arrays
   
@@ -488,14 +851,44 @@ int buildRandomBilayer( void ){
     nAtoms += waterNatoms;
   }
   
-  theConfig = simnfo[2].getConfiguration();
+  theConfig = mainInfo->getConfiguration();
   theConfig->createArrays( nAtoms );
-  simnfo[2].atoms = new Atom*[nAtoms];
-  atoms = simnfo[2].atoms;
-  simnfo[2].n_atoms = nAtoms;
+  mainInfo->atoms = new Atom*[nAtoms];
+  atoms = mainInfo->atoms;
+  mainInfo->n_atoms = nAtoms;
   
-  // initialize lipid positions
+  // set up the SimInfo object
+  
+  double Hmat[3][3];
 
+  Hmat[0][0] = box_x;
+  Hmat[0][1] = 0.0;
+  Hmat[0][2] = 0.0;
+
+  Hmat[1][0] = 0.0;
+  Hmat[1][1] = box_y;
+  Hmat[1][2] = 0.0;
+
+  Hmat[2][0] = 0.0;
+  Hmat[2][1] = 0.0;
+  Hmat[2][2] = box_z;
+    
+  mainInfo->setBoxM( Hmat );
+
+  // center the system on (0,0,0)
+
+  for(j=0;j<nLipids;j++){
+
+    mainInfo->wrapVector( lipidSites[j].pos );
+  }
+
+  for(j=0;j<newWaters;j++){
+
+    mainInfo->wrapVector( waterSites[j].pos );
+  }
+  
+  // initialize lipid positions
+  
   molIndex = 0;
   for(i=0; i<nLipids; i++ ){
     lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
@@ -516,58 +909,13 @@ int buildRandomBilayer( void ){
     }
   }  
 
-  // set up the SimInfo object
-  
-  double Hmat[3][3];
 
-  Hmat[0][0] = box_x;
-  Hmat[0][1] = 0.0;
-  Hmat[0][2] = 0.0;
 
-  Hmat[1][0] = 0.0;
-  Hmat[1][1] = box_y;
-  Hmat[1][2] = 0.0;
-
-  Hmat[2][0] = 0.0;
-  Hmat[2][1] = 0.0;
-  Hmat[2][2] = box_z;
-  
-
-  bsInfo.boxX = box_x;
-  bsInfo.boxY = box_y;
-  bsInfo.boxZ = box_z;
-  
-  simnfo[2].setBoxM( Hmat );
-
-  sprintf( simnfo[2].sampleName, "%s.dump", bsInfo.outPrefix );
-  sprintf( simnfo[2].finalName, "%s.init", bsInfo.outPrefix );
-
   // set up the writer and write out
   
-  writer = new DumpWriter( &simnfo[2] );
+  writer = new DumpWriter( mainInfo );
   writer->writeFinal( 0.0 );
 	
-  // clean up the memory
-  
-  //	 if( molMap != NULL )	delete[] molMap;
-  //	 if( cardDeck != NULL ) delete[] cardDeck;
-  //	 if( locate != NULL ){
-  //	   for(i=0; i<bsInfo.nComponents; i++){
-  //		 delete locate[i];
-  //	   }
-  //	   delete[] locate;
-  //	 }
-  //	 if( atoms != NULL ){
-  //	   for(i=0; i<nAtoms; i++){
-  //		 delete atoms[i];
-  //	   }
-  //	   Atom::destroyArrays();
-  //	   delete[] atoms;
-  //	 }
-  //	 if( molSeq != NULL ) delete[] molSeq;
-  //	 if( simnfo != NULL ) delete simnfo;
-  //	 if( writer != NULL ) delete writer;
-
   return 1;
 }
 
@@ -583,3 +931,37 @@ void buildMap( double &x, double &y, double &z, 
   if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ)  - 0.5 ) );
   else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
 }
+
+/***************************************************************************
+ * 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"
+
+		"\n"
+		"   long:\n"
+		"   -----\n"
+
+		"   --version          displays the version number\n"
+		"   --help             displays this help message.\n"
+		"   --waterBuffer <#>  sets the distance of closest approach of the water around the lipid\n"
+		"                        defaults to 5.0\n"
+		"   --lipidBuffer <#>  sets the distance of closest approach between two lipids\n"
+		"                        defaults to 6.0\n"
+		"   --waterLattice <#> sets the water lattice spacing\n"
+		"                        defaults to 4.929 ( 1 g/cm^3 )\n"
+		"\n"
+		"\n",
+		programName);
+}