| 1 | 
– | 
 | 
| 1 | 
  | 
#include <iostream> | 
| 2 | 
  | 
#include <vector> | 
| 3 | 
  | 
#include <algorithm> | 
| 4 | 
  | 
 | 
| 5 | 
< | 
#include <cstdlib> | 
| 6 | 
< | 
#include <cstring> | 
| 7 | 
< | 
#include <cmath> | 
| 5 | 
> | 
#include <stdio.h> | 
| 6 | 
> | 
#include <stdlib.h> | 
| 7 | 
> | 
#include <string.h> | 
| 8 | 
> | 
#include <math.h> | 
| 9 | 
  | 
 | 
| 10 | 
– | 
 | 
| 10 | 
  | 
#include "simError.h" | 
| 11 | 
  | 
#include "SimInfo.hpp" | 
| 12 | 
  | 
#include "ReadWrite.hpp" | 
| 13 | 
+ | 
#include "SimSetup.hpp" | 
| 14 | 
  | 
 | 
| 15 | 
  | 
#include "MoLocator.hpp" | 
| 16 | 
– | 
#include "sysBuild.hpp" | 
| 17 | 
– | 
#include "bilayerSys.hpp" | 
| 18 | 
– | 
 | 
| 16 | 
  | 
#include "latticeBuilder.hpp" | 
| 17 | 
  | 
 | 
| 18 | 
+ | 
#define RAND_SEED 31337 // \/\/007! | 
| 19 | 
+ | 
 | 
| 20 | 
+ | 
#define VERSION_MAJOR 0 | 
| 21 | 
+ | 
#define VERSION_MINOR 1 | 
| 22 | 
+ | 
 | 
| 23 | 
  | 
class SortCond{ | 
| 24 | 
  | 
   | 
| 25 | 
  | 
public: | 
| 34 | 
  | 
void buildMap( double &x, double &y, double &z, | 
| 35 | 
  | 
               double boxX, double boxY, double boxZ ); | 
| 36 | 
  | 
 | 
| 37 | 
< | 
int buildRandomBilayer( void ); | 
| 37 | 
> | 
int buildRandomBilayer( double waterCell, | 
| 38 | 
> | 
                        double waterBuffer, | 
| 39 | 
> | 
                        double lipidBuffer, | 
| 40 | 
> | 
                        char* waterName, | 
| 41 | 
> | 
                        char* lipidName ); | 
| 42 | 
  | 
 | 
| 43 | 
< | 
int buildBilayer( int isRandom ){ | 
| 43 | 
> | 
 | 
| 44 | 
> | 
char *programName; /*the name of the program */ | 
| 45 | 
> | 
void usage(void); | 
| 46 | 
> | 
using namespace std; | 
| 47 | 
> | 
SimInfo* mainInfo; | 
| 48 | 
> | 
 | 
| 49 | 
> | 
int main(int argC,char* argV[]){ | 
| 50 | 
  | 
   | 
| 51 | 
< | 
  if( isRandom ){ | 
| 40 | 
< | 
    return buildRandomBilayer(); | 
| 41 | 
< | 
  } | 
| 42 | 
< | 
  else{ | 
| 51 | 
> | 
  int i,j; // loop counters | 
| 52 | 
  | 
 | 
| 53 | 
< | 
    std::cerr << "unsupported feature\n"; | 
| 54 | 
< | 
    return 0; | 
| 53 | 
> | 
  char* outPrefix; // the output prefix | 
| 54 | 
> | 
 | 
| 55 | 
> | 
  char* conversionCheck; | 
| 56 | 
> | 
  bool conversionError; | 
| 57 | 
> | 
  bool optionError; | 
| 58 | 
> | 
 | 
| 59 | 
> | 
  char currentFlag; // used in parsing the flags | 
| 60 | 
> | 
  bool done = false; // multipurpose boolean | 
| 61 | 
> | 
  bool havePrefix; // boolean for the output prefix   | 
| 62 | 
> | 
 | 
| 63 | 
> | 
  char* lipidName; | 
| 64 | 
> | 
  char* waterName; | 
| 65 | 
> | 
  bool haveWaterName, haveLipidName; | 
| 66 | 
> | 
 | 
| 67 | 
> | 
  double waterLattice, waterBuffer, lipidBuffer; | 
| 68 | 
> | 
   | 
| 69 | 
> | 
  char* inName; | 
| 70 | 
> | 
   | 
| 71 | 
> | 
  SimSetup* simInit; | 
| 72 | 
> | 
   | 
| 73 | 
> | 
  // first things first, all of the initializations | 
| 74 | 
> | 
 | 
| 75 | 
> | 
  fflush(stdout); | 
| 76 | 
> | 
  srand48( 1337 );          // the random number generator. | 
| 77 | 
> | 
  initSimError();           // the error handler | 
| 78 | 
> | 
 | 
| 79 | 
> | 
  outPrefix = NULL; | 
| 80 | 
> | 
  inName = NULL; | 
| 81 | 
> | 
   | 
| 82 | 
> | 
  conversionError = false; | 
| 83 | 
> | 
  optionError = false; | 
| 84 | 
> | 
  | 
| 85 | 
> | 
  havePrefix = false; | 
| 86 | 
> | 
   | 
| 87 | 
> | 
  waterBuffer = 5.0; | 
| 88 | 
> | 
  lipidBuffer = 6.0; | 
| 89 | 
> | 
  waterLattice = 4.929; | 
| 90 | 
> | 
 | 
| 91 | 
> | 
  programName = argV[0]; /*save the program name in case we need it*/ | 
| 92 | 
> | 
 | 
| 93 | 
> | 
  for( i = 1; i < argC; i++){ | 
| 94 | 
> | 
     | 
| 95 | 
> | 
    if(argV[i][0] =='-'){ | 
| 96 | 
> | 
 | 
| 97 | 
> | 
      // parse the option | 
| 98 | 
> | 
       | 
| 99 | 
> | 
      if(argV[i][1] == '-' ){ | 
| 100 | 
> | 
 | 
| 101 | 
> | 
        // parse long word options | 
| 102 | 
> | 
         | 
| 103 | 
> | 
        if( !strcmp( argV[i], "--version") ){ | 
| 104 | 
> | 
           | 
| 105 | 
> | 
          printf("\n" | 
| 106 | 
> | 
                 "randomBilayer version %d.%d\n" | 
| 107 | 
> | 
                 "\n", | 
| 108 | 
> | 
                 VERSION_MAJOR, VERSION_MINOR ); | 
| 109 | 
> | 
          exit(0); | 
| 110 | 
> | 
           | 
| 111 | 
> | 
        } | 
| 112 | 
> | 
 | 
| 113 | 
> | 
        else if( !strcmp( argV[i], "--help") ){ | 
| 114 | 
> | 
           | 
| 115 | 
> | 
          usage(); | 
| 116 | 
> | 
          exit(0); | 
| 117 | 
> | 
        } | 
| 118 | 
> | 
 | 
| 119 | 
> | 
        else if (!strcmp( argV[i], "--lipidBuffer" )){ | 
| 120 | 
> | 
           | 
| 121 | 
> | 
          i++; | 
| 122 | 
> | 
          if( i>=argC ){ | 
| 123 | 
> | 
            sprintf( painCave.errMsg, | 
| 124 | 
> | 
                     "\n" | 
| 125 | 
> | 
                     "not enough arguments for -lipidBuffer\n"); | 
| 126 | 
> | 
            usage(); | 
| 127 | 
> | 
            painCave.isFatal = 1; | 
| 128 | 
> | 
            simError(); | 
| 129 | 
> | 
          }        | 
| 130 | 
> | 
           | 
| 131 | 
> | 
          lipidBuffer = atof( argV[i] ); | 
| 132 | 
> | 
        } | 
| 133 | 
> | 
 | 
| 134 | 
> | 
        else if (!strcmp( argV[i], "--waterBuffer" )){ | 
| 135 | 
> | 
           | 
| 136 | 
> | 
          i++; | 
| 137 | 
> | 
          if( i>=argC ){ | 
| 138 | 
> | 
            sprintf( painCave.errMsg, | 
| 139 | 
> | 
                     "\n" | 
| 140 | 
> | 
                     "not enough arguments for --waterBuffer\n"); | 
| 141 | 
> | 
            usage(); | 
| 142 | 
> | 
            painCave.isFatal = 1; | 
| 143 | 
> | 
            simError(); | 
| 144 | 
> | 
          }        | 
| 145 | 
> | 
           | 
| 146 | 
> | 
          waterBuffer = atof( argV[i] ); | 
| 147 | 
> | 
        } | 
| 148 | 
> | 
 | 
| 149 | 
> | 
        else if (!strcmp( argV[i], "--waterLattice" )){ | 
| 150 | 
> | 
           | 
| 151 | 
> | 
          i++; | 
| 152 | 
> | 
          if( i>=argC ){ | 
| 153 | 
> | 
            sprintf( painCave.errMsg, | 
| 154 | 
> | 
                     "\n" | 
| 155 | 
> | 
                     "not enough arguments for -waterLattice\n"); | 
| 156 | 
> | 
            usage(); | 
| 157 | 
> | 
            painCave.isFatal = 1; | 
| 158 | 
> | 
            simError(); | 
| 159 | 
> | 
          }        | 
| 160 | 
> | 
           | 
| 161 | 
> | 
          waterLattice = atof( argV[i] ); | 
| 162 | 
> | 
        } | 
| 163 | 
> | 
           | 
| 164 | 
> | 
           | 
| 165 | 
> | 
         | 
| 166 | 
> | 
        // anything else is an error | 
| 167 | 
> | 
 | 
| 168 | 
> | 
        else{ | 
| 169 | 
> | 
          fprintf( stderr,  | 
| 170 | 
> | 
                   "Invalid option \"%s\"\n", argV[i] ); | 
| 171 | 
> | 
          usage(); | 
| 172 | 
> | 
          exit(0); | 
| 173 | 
> | 
        } | 
| 174 | 
> | 
      } | 
| 175 | 
> | 
       | 
| 176 | 
> | 
      else{ | 
| 177 | 
> | 
         | 
| 178 | 
> | 
        // parse single character options | 
| 179 | 
> | 
         | 
| 180 | 
> | 
        done =0; | 
| 181 | 
> | 
        j = 1; | 
| 182 | 
> | 
        currentFlag = argV[i][j]; | 
| 183 | 
> | 
        while( (currentFlag != '\0') && (!done) ){ | 
| 184 | 
> | 
           | 
| 185 | 
> | 
          switch(currentFlag){ | 
| 186 | 
> | 
 | 
| 187 | 
> | 
          case 'o': | 
| 188 | 
> | 
            // -o <prefix> => the output prefix. | 
| 189 | 
> | 
 | 
| 190 | 
> | 
            j++; | 
| 191 | 
> | 
            currentFlag = argV[i][j]; | 
| 192 | 
> | 
               | 
| 193 | 
> | 
            if( currentFlag != '\0' ) optionError = true; | 
| 194 | 
> | 
             | 
| 195 | 
> | 
            if( optionError ){ | 
| 196 | 
> | 
              sprintf( painCave.errMsg,  | 
| 197 | 
> | 
                       "\n" | 
| 198 | 
> | 
                       "The -o flag should end an option sequence.\n" | 
| 199 | 
> | 
                       "   example: -r <outname> *NOT* -or <outname>\n" ); | 
| 200 | 
> | 
              usage(); | 
| 201 | 
> | 
              painCave.isFatal = 1; | 
| 202 | 
> | 
              simError(); | 
| 203 | 
> | 
            } | 
| 204 | 
> | 
               | 
| 205 | 
> | 
            i++; | 
| 206 | 
> | 
            if( i>=argC ){ | 
| 207 | 
> | 
              sprintf( painCave.errMsg,  | 
| 208 | 
> | 
                       "\n" | 
| 209 | 
> | 
                       "not enough arguments for -o\n"); | 
| 210 | 
> | 
              usage(); | 
| 211 | 
> | 
              painCave.isFatal = 1; | 
| 212 | 
> | 
              simError(); | 
| 213 | 
> | 
            }      | 
| 214 | 
> | 
             | 
| 215 | 
> | 
            outPrefix = argV[i]; | 
| 216 | 
> | 
            if( outPrefix[0] == '-' ) optionError = true; | 
| 217 | 
> | 
                 | 
| 218 | 
> | 
            if( optionError ){ | 
| 219 | 
> | 
              sprintf( painCave.errMsg,  | 
| 220 | 
> | 
                       "\n" | 
| 221 | 
> | 
                       "\"%s\" is not a valid out prefix/name.\n" | 
| 222 | 
> | 
                       "Out prefix/name should not begin with a dash.\n", | 
| 223 | 
> | 
                       outPrefix ); | 
| 224 | 
> | 
              usage(); | 
| 225 | 
> | 
              painCave.isFatal = 1; | 
| 226 | 
> | 
              simError(); | 
| 227 | 
> | 
            } | 
| 228 | 
> | 
             | 
| 229 | 
> | 
            havePrefix = true; | 
| 230 | 
> | 
            done = true; | 
| 231 | 
> | 
            break; | 
| 232 | 
> | 
 | 
| 233 | 
> | 
          case 'l': | 
| 234 | 
> | 
            // -l <lipidName> => the lipid name. | 
| 235 | 
> | 
 | 
| 236 | 
> | 
            j++; | 
| 237 | 
> | 
            currentFlag = argV[i][j]; | 
| 238 | 
> | 
               | 
| 239 | 
> | 
            if( currentFlag != '\0' ) optionError = true; | 
| 240 | 
> | 
             | 
| 241 | 
> | 
            if( optionError ){ | 
| 242 | 
> | 
              sprintf( painCave.errMsg,  | 
| 243 | 
> | 
                       "\n" | 
| 244 | 
> | 
                       "The -l flag should end an option sequence.\n" | 
| 245 | 
> | 
                       "   example: -rl <lipidName> *NOT* -lr <lipidName>\n" ); | 
| 246 | 
> | 
              usage(); | 
| 247 | 
> | 
              painCave.isFatal = 1; | 
| 248 | 
> | 
              simError(); | 
| 249 | 
> | 
            } | 
| 250 | 
> | 
               | 
| 251 | 
> | 
            i++; | 
| 252 | 
> | 
            if( i>=argC ){ | 
| 253 | 
> | 
              sprintf( painCave.errMsg,  | 
| 254 | 
> | 
                       "\n" | 
| 255 | 
> | 
                       "not enough arguments for -l\n"); | 
| 256 | 
> | 
              usage(); | 
| 257 | 
> | 
              painCave.isFatal = 1; | 
| 258 | 
> | 
              simError(); | 
| 259 | 
> | 
            }      | 
| 260 | 
> | 
             | 
| 261 | 
> | 
            lipidName = argV[i]; | 
| 262 | 
> | 
            if( lipidName[0] == '-' ) optionError = true; | 
| 263 | 
> | 
                 | 
| 264 | 
> | 
            if( optionError ){ | 
| 265 | 
> | 
              sprintf( painCave.errMsg,  | 
| 266 | 
> | 
                       "\n" | 
| 267 | 
> | 
                       "\"%s\" is not a valid lipidName.\n" | 
| 268 | 
> | 
                       "lipidName should not begin with a dash.\n", | 
| 269 | 
> | 
                       lipidName ); | 
| 270 | 
> | 
              usage(); | 
| 271 | 
> | 
              painCave.isFatal = 1; | 
| 272 | 
> | 
              simError(); | 
| 273 | 
> | 
            } | 
| 274 | 
> | 
             | 
| 275 | 
> | 
            haveLipidName = true; | 
| 276 | 
> | 
            done = true; | 
| 277 | 
> | 
            break; | 
| 278 | 
> | 
 | 
| 279 | 
> | 
          case 'w': | 
| 280 | 
> | 
            // -w <waterName> => the water name. | 
| 281 | 
> | 
 | 
| 282 | 
> | 
            j++; | 
| 283 | 
> | 
            currentFlag = argV[i][j]; | 
| 284 | 
> | 
               | 
| 285 | 
> | 
            if( currentFlag != '\0' ) optionError = true; | 
| 286 | 
> | 
             | 
| 287 | 
> | 
            if( optionError ){ | 
| 288 | 
> | 
              sprintf( painCave.errMsg,  | 
| 289 | 
> | 
                       "\n" | 
| 290 | 
> | 
                       "The -w flag should end an option sequence.\n" | 
| 291 | 
> | 
                       "   example: -rw <waterName> *NOT* -lw <waterName>\n" ); | 
| 292 | 
> | 
              usage(); | 
| 293 | 
> | 
              painCave.isFatal = 1; | 
| 294 | 
> | 
              simError(); | 
| 295 | 
> | 
            } | 
| 296 | 
> | 
               | 
| 297 | 
> | 
            i++; | 
| 298 | 
> | 
            if( i>=argC ){ | 
| 299 | 
> | 
              sprintf( painCave.errMsg,  | 
| 300 | 
> | 
                       "\n" | 
| 301 | 
> | 
                       "not enough arguments for -w\n"); | 
| 302 | 
> | 
              usage(); | 
| 303 | 
> | 
              painCave.isFatal = 1; | 
| 304 | 
> | 
              simError(); | 
| 305 | 
> | 
            }      | 
| 306 | 
> | 
             | 
| 307 | 
> | 
            waterName = argV[i]; | 
| 308 | 
> | 
            if( waterName[0] == '-' ) optionError = true; | 
| 309 | 
> | 
                 | 
| 310 | 
> | 
            if( optionError ){ | 
| 311 | 
> | 
              sprintf( painCave.errMsg,  | 
| 312 | 
> | 
                       "\n" | 
| 313 | 
> | 
                       "\"%s\" is not a valid waterName.\n" | 
| 314 | 
> | 
                       "waterName should not begin with a dash.\n", | 
| 315 | 
> | 
                       waterName ); | 
| 316 | 
> | 
              usage(); | 
| 317 | 
> | 
              painCave.isFatal = 1; | 
| 318 | 
> | 
              simError(); | 
| 319 | 
> | 
            } | 
| 320 | 
> | 
             | 
| 321 | 
> | 
            haveWaterName = true; | 
| 322 | 
> | 
            done = true; | 
| 323 | 
> | 
            break; | 
| 324 | 
> | 
 | 
| 325 | 
> | 
          default: | 
| 326 | 
> | 
 | 
| 327 | 
> | 
            sprintf(painCave.errMsg,  | 
| 328 | 
> | 
                    "\n" | 
| 329 | 
> | 
                    "Bad option \"-%c\"\n", currentFlag); | 
| 330 | 
> | 
            usage(); | 
| 331 | 
> | 
            painCave.isFatal = 1; | 
| 332 | 
> | 
            simError(); | 
| 333 | 
> | 
          } | 
| 334 | 
> | 
          j++; | 
| 335 | 
> | 
          currentFlag = argV[i][j]; | 
| 336 | 
> | 
        } | 
| 337 | 
> | 
      } | 
| 338 | 
> | 
    } | 
| 339 | 
> | 
 | 
| 340 | 
> | 
    else{ | 
| 341 | 
> | 
       | 
| 342 | 
> | 
      if( inName != NULL ){ | 
| 343 | 
> | 
        sprintf( painCave.errMsg, | 
| 344 | 
> | 
                 "Error at \"%s\", program does not currently support\n" | 
| 345 | 
> | 
                 "more than one input bass file.\n" | 
| 346 | 
> | 
                 "\n", | 
| 347 | 
> | 
                 argV[i]); | 
| 348 | 
> | 
        usage(); | 
| 349 | 
> | 
        painCave.isFatal = 1; | 
| 350 | 
> | 
        simError(); | 
| 351 | 
> | 
      } | 
| 352 | 
> | 
       | 
| 353 | 
> | 
      inName = argV[i]; | 
| 354 | 
> | 
       | 
| 355 | 
> | 
    } | 
| 356 | 
> | 
  }   | 
| 357 | 
> | 
   | 
| 358 | 
> | 
  if( inName == NULL ){ | 
| 359 | 
> | 
    sprintf( painCave.errMsg, | 
| 360 | 
> | 
             "Error, bass file is needed to run.\n" ); | 
| 361 | 
> | 
    usage(); | 
| 362 | 
> | 
    painCave.isFatal = 1; | 
| 363 | 
> | 
    simError(); | 
| 364 | 
  | 
  } | 
| 47 | 
– | 
} | 
| 365 | 
  | 
 | 
| 366 | 
+ | 
  // if no output prefix is given default to "donkey". | 
| 367 | 
  | 
 | 
| 368 | 
< | 
int buildRandomBilayer( void ){ | 
| 368 | 
> | 
  if( !havePrefix ){ | 
| 369 | 
> | 
    outPrefix = strdup( "donkey" ); | 
| 370 | 
> | 
  } | 
| 371 | 
  | 
 | 
| 372 | 
+ | 
  if( !haveWaterName ){ | 
| 373 | 
+ | 
    sprintf( painCave.errMsg, | 
| 374 | 
+ | 
             "Error, the water name is needed to run.\n" | 
| 375 | 
+ | 
             ); | 
| 376 | 
+ | 
    usage(); | 
| 377 | 
+ | 
    painCave.isFatal = 1; | 
| 378 | 
+ | 
    simError(); | 
| 379 | 
+ | 
  }   | 
| 380 | 
+ | 
   | 
| 381 | 
+ | 
  if( !haveLipidName ){ | 
| 382 | 
+ | 
    sprintf( painCave.errMsg, | 
| 383 | 
+ | 
             "Error, the lipid name is needed to run.\n" | 
| 384 | 
+ | 
             ); | 
| 385 | 
+ | 
    usage(); | 
| 386 | 
+ | 
    painCave.isFatal = 1; | 
| 387 | 
+ | 
    simError(); | 
| 388 | 
+ | 
  } | 
| 389 | 
+ | 
   | 
| 390 | 
+ | 
   | 
| 391 | 
+ | 
  // create and initialize the info object | 
| 392 | 
+ | 
 | 
| 393 | 
+ | 
  mainInfo = new SimInfo(); | 
| 394 | 
+ | 
  simInit = new SimSetup(); | 
| 395 | 
+ | 
  simInit->setSimInfo( mainInfo ); | 
| 396 | 
+ | 
  simInit->suspendInit(); | 
| 397 | 
+ | 
  simInit->parseFile( inName ); | 
| 398 | 
+ | 
  simInit->createSim(); | 
| 399 | 
+ | 
 | 
| 400 | 
+ | 
  delete simInit; | 
| 401 | 
+ | 
 | 
| 402 | 
+ | 
  sprintf( mainInfo->statusName, "%s.stat", outPrefix ); | 
| 403 | 
+ | 
  sprintf( mainInfo->sampleName, "%s.dump", outPrefix ); | 
| 404 | 
+ | 
  sprintf( mainInfo->finalName, "%s.init", outPrefix ); | 
| 405 | 
+ | 
 | 
| 406 | 
+ | 
  buildRandomBilayer( waterLattice, waterBuffer, lipidBuffer,  | 
| 407 | 
+ | 
                      waterName, lipidName ); | 
| 408 | 
+ | 
 | 
| 409 | 
+ | 
  return 0;   | 
| 410 | 
+ | 
} | 
| 411 | 
+ | 
 | 
| 412 | 
+ | 
int buildRandomBilayer( double waterCell, | 
| 413 | 
+ | 
                        double water_padding, | 
| 414 | 
+ | 
                        double lipid_spaceing, | 
| 415 | 
+ | 
                        char* waterName, | 
| 416 | 
+ | 
                        char* lipidName ){ | 
| 417 | 
+ | 
 | 
| 418 | 
  | 
  typedef struct{ | 
| 419 | 
  | 
    double rot[3][3]; | 
| 420 | 
  | 
    double pos[3]; | 
| 421 | 
  | 
  } coord; | 
| 422 | 
  | 
 | 
| 57 | 
– | 
 | 
| 58 | 
– | 
 | 
| 59 | 
– | 
  const double waterRho = 0.0334; // number density per cubic angstrom | 
| 60 | 
– | 
  const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters | 
| 61 | 
– | 
  const double waterCell = 4.929; // fcc unit cell length | 
| 62 | 
– | 
 | 
| 423 | 
  | 
  Lattice myFCC( FCC_LATTICE_TYPE, waterCell ); | 
| 424 | 
  | 
  double *posX, *posY, *posZ; | 
| 425 | 
  | 
  double pos[3], posA[3], posB[3]; | 
| 426 | 
  | 
 | 
| 67 | 
– | 
  const double water_padding = 6.0; | 
| 68 | 
– | 
  const double lipid_spaceing = 8.0; | 
| 69 | 
– | 
 | 
| 70 | 
– | 
 | 
| 427 | 
  | 
  int i,j,k, l, m; | 
| 428 | 
  | 
  int nAtoms, atomIndex, molIndex, molID; | 
| 429 | 
  | 
  int* molSeq; | 
| 449 | 
  | 
  int foundLipid, foundWater; | 
| 450 | 
  | 
  int nLipids, lipidNatoms, nWaters, waterNatoms; | 
| 451 | 
  | 
  double testBox, maxLength; | 
| 96 | 
– | 
   | 
| 97 | 
– | 
  srand48( RAND_SEED ); | 
| 452 | 
  | 
 | 
| 453 | 
+ | 
  double waterRho, waterVol; | 
| 454 | 
  | 
 | 
| 455 | 
< | 
  // create the simInfo objects | 
| 456 | 
< | 
 | 
| 102 | 
< | 
  simnfo = new SimInfo[3]; | 
| 455 | 
> | 
   | 
| 456 | 
> | 
  srand48( RAND_SEED ); | 
| 457 | 
  | 
 | 
| 458 | 
+ | 
  // calculate the water parameters | 
| 459 | 
+ | 
   | 
| 460 | 
+ | 
  waterVol = waterCell * waterCell * waterCell; | 
| 461 | 
+ | 
  waterRho = 4.0 / waterVol; | 
| 462 | 
  | 
 | 
| 463 | 
  | 
  // set the the lipidStamp | 
| 464 | 
  | 
 | 
| 465 | 
  | 
  foundLipid = 0; | 
| 466 | 
  | 
  foundWater = 0; | 
| 467 | 
< | 
  for(i=0; i<bsInfo.nComponents; i++){ | 
| 468 | 
< | 
    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){ | 
| 467 | 
> | 
  for(i=0; i<mainInfo->nComponents; i++){ | 
| 468 | 
> | 
     | 
| 469 | 
> | 
    if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){ | 
| 470 | 
  | 
           | 
| 471 | 
  | 
      foundLipid = 1; | 
| 472 | 
< | 
      lipidStamp = bsInfo.compStamps[i]; | 
| 473 | 
< | 
      nLipids = bsInfo.componentsNmol[i]; | 
| 472 | 
> | 
      lipidStamp = mainInfo->compStamps[i]; | 
| 473 | 
> | 
      nLipids = mainInfo->componentsNmol[i]; | 
| 474 | 
  | 
    } | 
| 475 | 
< | 
    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){ | 
| 475 | 
> | 
    if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){ | 
| 476 | 
  | 
           | 
| 477 | 
  | 
      foundWater = 1; | 
| 478 | 
  | 
           | 
| 479 | 
< | 
      waterStamp = bsInfo.compStamps[i]; | 
| 480 | 
< | 
      nWaters = bsInfo.componentsNmol[i]; | 
| 479 | 
> | 
      waterStamp = mainInfo->compStamps[i]; | 
| 480 | 
> | 
      nWaters = mainInfo->componentsNmol[i]; | 
| 481 | 
  | 
    } | 
| 482 | 
  | 
  } | 
| 483 | 
  | 
  if( !foundLipid ){ | 
| 484 | 
  | 
    sprintf(painCave.errMsg, | 
| 485 | 
< | 
            "Could not find lipid \"%s\" in the bass file.\n", | 
| 486 | 
< | 
            bsInfo.lipidName ); | 
| 485 | 
> | 
            "randomBilayer error: Could not find lipid \"%s\" in the bass file.\n", | 
| 486 | 
> | 
            lipidName ); | 
| 487 | 
  | 
    painCave.isFatal = 1; | 
| 488 | 
  | 
    simError(); | 
| 489 | 
  | 
  } | 
| 490 | 
  | 
  if( !foundWater ){ | 
| 491 | 
  | 
    sprintf(painCave.errMsg, | 
| 492 | 
< | 
            "Could not find solvent \"%s\" in the bass file.\n", | 
| 493 | 
< | 
            bsInfo.waterName ); | 
| 492 | 
> | 
            "randomBilayer error: Could not find solvent \"%s\" in the bass file.\n", | 
| 493 | 
> | 
            waterName ); | 
| 494 | 
  | 
    painCave.isFatal = 1; | 
| 495 | 
  | 
    simError(); | 
| 496 | 
  | 
  } | 
| 504 | 
  | 
  waterLocate = new MoLocator( waterStamp ); | 
| 505 | 
  | 
  waterNatoms = waterStamp->getNAtoms(); | 
| 506 | 
  | 
   | 
| 507 | 
< | 
  nAtoms = lipidNatoms; | 
| 507 | 
> | 
  // create a test box | 
| 508 | 
  | 
 | 
| 509 | 
< | 
  simnfo[0].n_atoms = nAtoms; | 
| 151 | 
< | 
  simnfo[0].atoms=new Atom*[nAtoms]; | 
| 152 | 
< | 
  | 
| 153 | 
< | 
  theConfig = simnfo[0].getConfiguration(); | 
| 154 | 
< | 
  theConfig->createArrays( simnfo[0].n_atoms ); | 
| 509 | 
> | 
  SimInfo* testInfo; | 
| 510 | 
  | 
 | 
| 511 | 
< | 
  atoms=simnfo[0].atoms; | 
| 512 | 
< | 
         | 
| 511 | 
> | 
  testInfo = new SimInfo(); | 
| 512 | 
> | 
  testInfo->n_atoms = lipidNatoms; | 
| 513 | 
> | 
  theConfig = testInfo->getConfiguration(); | 
| 514 | 
> | 
  theConfig->createArrays( lipidNatoms ); | 
| 515 | 
> | 
  testInfo->atoms = new Atom*[lipidNatoms]; | 
| 516 | 
> | 
  atoms = testInfo->atoms; | 
| 517 | 
  | 
 | 
| 518 | 
  | 
  // create the test box for initial water displacement | 
| 519 | 
  | 
 | 
| 520 | 
< | 
  testBox = maxLength + waterCell * 4.0; // pad with 4 cells | 
| 520 | 
> | 
  testBox = maxLength + waterCell * 10.0; // pad with 4 cells | 
| 521 | 
  | 
  nCells = (int)( testBox / waterCell + 1.0 ); | 
| 522 | 
  | 
  int testWaters = 4 * nCells * nCells * nCells; | 
| 523 | 
  | 
   | 
| 610 | 
  | 
  const double boxTargetX = 66.22752; | 
| 611 | 
  | 
  const double boxTargetY = 60.53088; | 
| 612 | 
  | 
 | 
| 613 | 
< | 
  nCellsX = (int)ceil(boxTargetX / waterCell); | 
| 614 | 
< | 
  nCellsY = (int)ceil(boxTargetY / waterCell); | 
| 613 | 
> | 
//   nCellsX = (int)ceil(boxTargetX / waterCell); | 
| 614 | 
> | 
//   nCellsY = (int)ceil(boxTargetY / waterCell); | 
| 615 | 
  | 
   | 
| 616 | 
  | 
  int testTot; | 
| 617 | 
  | 
  int done = 0; | 
| 618 | 
+ | 
  nCellsX = 0; | 
| 619 | 
+ | 
  nCellsY = 0; | 
| 620 | 
  | 
  nCellsZ = 0; | 
| 621 | 
  | 
  while( !done ){ | 
| 622 | 
  | 
         | 
| 623 | 
+ | 
    nCellsX++; | 
| 624 | 
+ | 
    nCellsY++; | 
| 625 | 
  | 
    nCellsZ++; | 
| 626 | 
  | 
    testTot = 4 * nCellsX * nCellsY * nCellsZ; | 
| 627 | 
  | 
         | 
| 667 | 
  | 
   | 
| 668 | 
  | 
  int reject; | 
| 669 | 
  | 
  int testDX, acceptedDX; | 
| 670 | 
+ | 
  SimInfo* testInfo2; | 
| 671 | 
  | 
 | 
| 672 | 
  | 
  nAtoms = nLipids * lipidNatoms; | 
| 673 | 
+ | 
  testInfo2 = new SimInfo(); | 
| 674 | 
+ | 
  testInfo2->n_atoms = nAtoms; | 
| 675 | 
+ | 
  theConfig = testInfo2->getConfiguration(); | 
| 676 | 
+ | 
  theConfig->createArrays( nAtoms ); | 
| 677 | 
+ | 
  testInfo2->atoms = new Atom*[nAtoms]; | 
| 678 | 
+ | 
  atoms = testInfo2->atoms; | 
| 679 | 
  | 
 | 
| 310 | 
– | 
  simnfo[1].n_atoms = nAtoms; | 
| 311 | 
– | 
  simnfo[1].atoms=new Atom*[nAtoms]; | 
| 312 | 
– | 
  | 
| 313 | 
– | 
  theConfig = simnfo[1].getConfiguration(); | 
| 314 | 
– | 
  theConfig->createArrays( simnfo[1].n_atoms ); | 
| 315 | 
– | 
 | 
| 316 | 
– | 
  atoms=simnfo[1].atoms; | 
| 317 | 
– | 
 | 
| 680 | 
  | 
  rCutSqr = lipid_spaceing * lipid_spaceing; | 
| 681 | 
  | 
 | 
| 682 | 
  | 
  for(i=0; i<nLipids; i++ ){ | 
| 729 | 
  | 
      else{ | 
| 730 | 
  | 
        done = 1; | 
| 731 | 
  | 
        std::cout << (i+1) << " has been accepted\n"; | 
| 732 | 
+ | 
        std::cout.flush(); | 
| 733 | 
  | 
      } | 
| 734 | 
  | 
    } | 
| 735 | 
  | 
  } | 
| 736 | 
  | 
         | 
| 737 | 
  | 
 | 
| 738 | 
< | 
  // zSort of the lipid positions | 
| 738 | 
> | 
//   // zSort of the lipid positions | 
| 739 | 
  | 
 | 
| 740 | 
  | 
 | 
| 741 | 
< | 
  vector< pair<int,double> >zSortArray; | 
| 742 | 
< | 
  for(i=0;i<nLipids;i++)  | 
| 743 | 
< | 
    zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) ); | 
| 741 | 
> | 
//   vector< pair<int,double> >zSortArray; | 
| 742 | 
> | 
//   for(i=0;i<nLipids;i++)  | 
| 743 | 
> | 
//     zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) ); | 
| 744 | 
  | 
 | 
| 745 | 
< | 
  sort(zSortArray.begin(),zSortArray.end(),SortCond()); | 
| 745 | 
> | 
//   sort(zSortArray.begin(),zSortArray.end(),SortCond()); | 
| 746 | 
  | 
 | 
| 747 | 
< | 
  ofstream outFile( "./zipper.bass", ios::app); | 
| 747 | 
> | 
//   ofstream outFile( "./zipper.bass", ios::app); | 
| 748 | 
  | 
 | 
| 749 | 
< | 
  for(i=0; i<nLipids; i++){ | 
| 750 | 
< | 
    outFile << "zConstraint[" << i << "]{\n" | 
| 751 | 
< | 
            << "  molIndex = " << zSortArray[i].first <<  ";\n" | 
| 752 | 
< | 
            << "  zPos = "; | 
| 749 | 
> | 
//   for(i=0; i<nLipids; i++){ | 
| 750 | 
> | 
//     outFile << "zConstraint[" << i << "]{\n" | 
| 751 | 
> | 
//          << "  molIndex = " << zSortArray[i].first <<  ";\n" | 
| 752 | 
> | 
//          << "  zPos = "; | 
| 753 | 
  | 
 | 
| 754 | 
< | 
    if(i<32) outFile << "60.0;\n"; | 
| 755 | 
< | 
    else outFile << "100.0;\n"; | 
| 754 | 
> | 
//     if(i<32) outFile << "60.0;\n"; | 
| 755 | 
> | 
//     else outFile << "100.0;\n"; | 
| 756 | 
  | 
 | 
| 757 | 
< | 
    outFile << "  kRatio = 0.5;\n" | 
| 758 | 
< | 
            << "}\n"; | 
| 759 | 
< | 
  } | 
| 757 | 
> | 
//     outFile << "  kRatio = 0.5;\n" | 
| 758 | 
> | 
//          << "}\n"; | 
| 759 | 
> | 
//   } | 
| 760 | 
  | 
   | 
| 761 | 
< | 
  outFile.close(); | 
| 761 | 
> | 
//   outFile.close(); | 
| 762 | 
  | 
 | 
| 763 | 
  | 
 | 
| 764 | 
  | 
  // cut out the waters that overlap with the lipids. | 
| 830 | 
  | 
 | 
| 831 | 
  | 
  // clean up our messes before building the final system. | 
| 832 | 
  | 
 | 
| 833 | 
< | 
  simnfo[0].getConfiguration()->destroyArrays(); | 
| 834 | 
< | 
  simnfo[1].getConfiguration()->destroyArrays(); | 
| 833 | 
> | 
  testInfo->getConfiguration()->destroyArrays(); | 
| 834 | 
> | 
  testInfo2->getConfiguration()->destroyArrays(); | 
| 835 | 
  | 
 | 
| 836 | 
  | 
  // create the real Atom arrays | 
| 837 | 
  | 
   | 
| 851 | 
  | 
    nAtoms += waterNatoms; | 
| 852 | 
  | 
  } | 
| 853 | 
  | 
   | 
| 854 | 
< | 
  theConfig = simnfo[2].getConfiguration(); | 
| 854 | 
> | 
  theConfig = mainInfo->getConfiguration(); | 
| 855 | 
  | 
  theConfig->createArrays( nAtoms ); | 
| 856 | 
< | 
  simnfo[2].atoms = new Atom*[nAtoms]; | 
| 857 | 
< | 
  atoms = simnfo[2].atoms; | 
| 858 | 
< | 
  simnfo[2].n_atoms = nAtoms; | 
| 856 | 
> | 
  mainInfo->atoms = new Atom*[nAtoms]; | 
| 857 | 
> | 
  atoms = mainInfo->atoms; | 
| 858 | 
> | 
  mainInfo->n_atoms = nAtoms; | 
| 859 | 
  | 
   | 
| 860 | 
< | 
  // initialize lipid positions | 
| 860 | 
> | 
  // set up the SimInfo object | 
| 861 | 
> | 
   | 
| 862 | 
> | 
  double Hmat[3][3]; | 
| 863 | 
  | 
 | 
| 864 | 
+ | 
  Hmat[0][0] = box_x; | 
| 865 | 
+ | 
  Hmat[0][1] = 0.0; | 
| 866 | 
+ | 
  Hmat[0][2] = 0.0; | 
| 867 | 
+ | 
 | 
| 868 | 
+ | 
  Hmat[1][0] = 0.0; | 
| 869 | 
+ | 
  Hmat[1][1] = box_y; | 
| 870 | 
+ | 
  Hmat[1][2] = 0.0; | 
| 871 | 
+ | 
 | 
| 872 | 
+ | 
  Hmat[2][0] = 0.0; | 
| 873 | 
+ | 
  Hmat[2][1] = 0.0; | 
| 874 | 
+ | 
  Hmat[2][2] = box_z; | 
| 875 | 
+ | 
     | 
| 876 | 
+ | 
  mainInfo->setBoxM( Hmat ); | 
| 877 | 
+ | 
 | 
| 878 | 
+ | 
  // center the system on (0,0,0) | 
| 879 | 
+ | 
 | 
| 880 | 
+ | 
  for(j=0;j<nLipids;j++){ | 
| 881 | 
+ | 
 | 
| 882 | 
+ | 
    mainInfo->wrapVector( lipidSites[j].pos ); | 
| 883 | 
+ | 
  } | 
| 884 | 
+ | 
 | 
| 885 | 
+ | 
  for(j=0;j<newWaters;j++){ | 
| 886 | 
+ | 
 | 
| 887 | 
+ | 
    mainInfo->wrapVector( waterSites[j].pos ); | 
| 888 | 
+ | 
  } | 
| 889 | 
+ | 
   | 
| 890 | 
+ | 
  // initialize lipid positions | 
| 891 | 
+ | 
   | 
| 892 | 
  | 
  molIndex = 0; | 
| 893 | 
  | 
  for(i=0; i<nLipids; i++ ){ | 
| 894 | 
  | 
    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms, | 
| 909 | 
  | 
    } | 
| 910 | 
  | 
  }   | 
| 911 | 
  | 
 | 
| 519 | 
– | 
  // set up the SimInfo object | 
| 520 | 
– | 
   | 
| 521 | 
– | 
  double Hmat[3][3]; | 
| 912 | 
  | 
 | 
| 523 | 
– | 
  Hmat[0][0] = box_x; | 
| 524 | 
– | 
  Hmat[0][1] = 0.0; | 
| 525 | 
– | 
  Hmat[0][2] = 0.0; | 
| 913 | 
  | 
 | 
| 527 | 
– | 
  Hmat[1][0] = 0.0; | 
| 528 | 
– | 
  Hmat[1][1] = box_y; | 
| 529 | 
– | 
  Hmat[1][2] = 0.0; | 
| 530 | 
– | 
 | 
| 531 | 
– | 
  Hmat[2][0] = 0.0; | 
| 532 | 
– | 
  Hmat[2][1] = 0.0; | 
| 533 | 
– | 
  Hmat[2][2] = box_z; | 
| 534 | 
– | 
   | 
| 535 | 
– | 
 | 
| 536 | 
– | 
  bsInfo.boxX = box_x; | 
| 537 | 
– | 
  bsInfo.boxY = box_y; | 
| 538 | 
– | 
  bsInfo.boxZ = box_z; | 
| 539 | 
– | 
   | 
| 540 | 
– | 
  simnfo[2].setBoxM( Hmat ); | 
| 541 | 
– | 
 | 
| 542 | 
– | 
  sprintf( simnfo[2].sampleName, "%s.dump", bsInfo.outPrefix ); | 
| 543 | 
– | 
  sprintf( simnfo[2].finalName, "%s.init", bsInfo.outPrefix ); | 
| 544 | 
– | 
 | 
| 914 | 
  | 
  // set up the writer and write out | 
| 915 | 
  | 
   | 
| 916 | 
< | 
  writer = new DumpWriter( &simnfo[2] ); | 
| 916 | 
> | 
  writer = new DumpWriter( mainInfo ); | 
| 917 | 
  | 
  writer->writeFinal( 0.0 ); | 
| 918 | 
  | 
         | 
| 550 | 
– | 
  // clean up the memory | 
| 551 | 
– | 
   | 
| 552 | 
– | 
  //     if( molMap != NULL )   delete[] molMap; | 
| 553 | 
– | 
  //     if( cardDeck != NULL ) delete[] cardDeck; | 
| 554 | 
– | 
  //     if( locate != NULL ){ | 
| 555 | 
– | 
  //       for(i=0; i<bsInfo.nComponents; i++){ | 
| 556 | 
– | 
  //             delete locate[i]; | 
| 557 | 
– | 
  //       } | 
| 558 | 
– | 
  //       delete[] locate; | 
| 559 | 
– | 
  //     } | 
| 560 | 
– | 
  //     if( atoms != NULL ){ | 
| 561 | 
– | 
  //       for(i=0; i<nAtoms; i++){ | 
| 562 | 
– | 
  //             delete atoms[i]; | 
| 563 | 
– | 
  //       } | 
| 564 | 
– | 
  //       Atom::destroyArrays(); | 
| 565 | 
– | 
  //       delete[] atoms; | 
| 566 | 
– | 
  //     } | 
| 567 | 
– | 
  //     if( molSeq != NULL ) delete[] molSeq; | 
| 568 | 
– | 
  //     if( simnfo != NULL ) delete simnfo; | 
| 569 | 
– | 
  //     if( writer != NULL ) delete writer; | 
| 570 | 
– | 
 | 
| 919 | 
  | 
  return 1; | 
| 920 | 
  | 
} | 
| 921 | 
  | 
 | 
| 931 | 
  | 
  if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ)  - 0.5 ) ); | 
| 932 | 
  | 
  else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5)); | 
| 933 | 
  | 
} | 
| 934 | 
+ | 
 | 
| 935 | 
+ | 
/*************************************************************************** | 
| 936 | 
+ | 
 * prints out the usage for the command line arguments, then exits. | 
| 937 | 
+ | 
 ***************************************************************************/ | 
| 938 | 
+ | 
 | 
| 939 | 
+ | 
void usage(){ | 
| 940 | 
+ | 
  (void)fprintf(stdout,  | 
| 941 | 
+ | 
                "\n" | 
| 942 | 
+ | 
                "The proper usage is: %s [options] <input_file>\n" | 
| 943 | 
+ | 
                "\n" | 
| 944 | 
+ | 
                "Options:\n" | 
| 945 | 
+ | 
                "\n" | 
| 946 | 
+ | 
                "   short:\n" | 
| 947 | 
+ | 
                "   ------\n" | 
| 948 | 
+ | 
                "   -o <name>       The output prefix\n" | 
| 949 | 
+ | 
                "   -l <lipidName>  The name of the lipid molecule specified in the BASS file.\n" | 
| 950 | 
+ | 
                "   -w <waterName>  The name of the water molecule specified in the BASS file.\n" | 
| 951 | 
+ | 
 | 
| 952 | 
+ | 
                "\n" | 
| 953 | 
+ | 
                "   long:\n" | 
| 954 | 
+ | 
                "   -----\n" | 
| 955 | 
+ | 
 | 
| 956 | 
+ | 
                "   --version          displays the version number\n" | 
| 957 | 
+ | 
                "   --help             displays this help message.\n" | 
| 958 | 
+ | 
                "   --waterBuffer <#>  sets the distance of closest approach of the water around the lipid\n" | 
| 959 | 
+ | 
                "                        defaults to 5.0\n" | 
| 960 | 
+ | 
                "   --lipidBuffer <#>  sets the distance of closest approach between two lipids\n" | 
| 961 | 
+ | 
                "                        defaults to 6.0\n" | 
| 962 | 
+ | 
                "   --waterLattice <#> sets the water lattice spacing\n" | 
| 963 | 
+ | 
                "                        defaults to 4.929 ( 1 g/cm^3 )\n" | 
| 964 | 
+ | 
                "\n" | 
| 965 | 
+ | 
                "\n", | 
| 966 | 
+ | 
                programName); | 
| 967 | 
+ | 
} |