ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-1.0/utils/sysbuilder/randomBilayer.cpp
Revision: 1417
Committed: Tue Jul 27 15:41:17 2004 UTC (20 years, 1 month ago) by tim
File size: 21405 byte(s)
Log Message:
BASS eradication project (part 1)

File Contents

# User Rev Content
1 gezelter 1334 #include <iostream>
2     #include <vector>
3     #include <algorithm>
4    
5     #include <stdio.h>
6     #include <stdlib.h>
7     #include <string.h>
8     #include <math.h>
9    
10     #include "simError.h"
11     #include "SimInfo.hpp"
12     #include "ReadWrite.hpp"
13     #include "SimSetup.hpp"
14    
15     #include "MoLocator.hpp"
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:
26     bool operator()(const pair<int, double>& p1, const pair<int, double>& p2){
27     return p1.second < p2.second;
28     }
29    
30    
31     };
32    
33    
34     void buildMap( double &x, double &y, double &z,
35     double boxX, double boxY, double boxZ );
36    
37     int buildRandomBilayer( double waterCell,
38     double waterBuffer,
39     double lipidBuffer,
40     char* waterName,
41     char* lipidName );
42    
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     int i,j; // loop counters
52    
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     }
365    
366     // if no output prefix is given default to "donkey".
367    
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    
423     Lattice myFCC( FCC_LATTICE_TYPE, waterCell );
424     double *posX, *posY, *posZ;
425     double pos[3], posA[3], posB[3];
426    
427     int i,j,k, l, m;
428     int nAtoms, atomIndex, molIndex, molID;
429     int* molSeq;
430     int* molMap;
431     int* molStart;
432     int* cardDeck;
433     int deckSize;
434     int rSite, rCard;
435     double cell;
436     int nCells, nSites, siteIndex;
437    
438     coord testSite;
439    
440     Atom** atoms;
441     SimInfo* simnfo;
442     SimState* theConfig;
443     DumpWriter* writer;
444    
445     MoleculeStamp* lipidStamp;
446     MoleculeStamp* waterStamp;
447     MoLocator *lipidLocate;
448     MoLocator *waterLocate;
449     int foundLipid, foundWater;
450     int nLipids, lipidNatoms, nWaters, waterNatoms;
451     double testBox, maxLength;
452    
453     double waterRho, waterVol;
454    
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<mainInfo->nComponents; i++){
468    
469     if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){
470    
471     foundLipid = 1;
472     lipidStamp = mainInfo->compStamps[i];
473     nLipids = mainInfo->componentsNmol[i];
474     }
475     if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){
476    
477     foundWater = 1;
478    
479     waterStamp = mainInfo->compStamps[i];
480     nWaters = mainInfo->componentsNmol[i];
481     }
482     }
483     if( !foundLipid ){
484     sprintf(painCave.errMsg,
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     "randomBilayer error: Could not find solvent \"%s\" in the bass file.\n",
493     waterName );
494     painCave.isFatal = 1;
495     simError();
496     }
497    
498     //create the temp Molocator and atom Arrays
499    
500     lipidLocate = new MoLocator( lipidStamp );
501     lipidNatoms = lipidStamp->getNAtoms();
502     maxLength = lipidLocate->getMaxLength();
503    
504     waterLocate = new MoLocator( waterStamp );
505     waterNatoms = waterStamp->getNAtoms();
506    
507     // create a test box
508    
509     SimInfo* testInfo;
510    
511     testInfo = new SimInfo();
512     testInfo->n_atoms = lipidNatoms;
513     theConfig = testInfo->getConfiguration();
514     theConfig->createArrays( lipidNatoms );
515     testInfo->atoms = new Atom*[lipidNatoms];
516 tim 1417 for (i=0; i < lipidNatoms; i++) {
517     testInfo->atoms[i] = new Atom(i, theConfig);
518     testInfo->atoms[i]->setCoords();
519     }
520 gezelter 1334 atoms = testInfo->atoms;
521    
522     // create the test box for initial water displacement
523    
524     testBox = maxLength + waterCell * 10.0; // pad with 4 cells
525     nCells = (int)( testBox / waterCell + 1.0 );
526     int testWaters = 4 * nCells * nCells * nCells;
527    
528     double* waterX = new double[testWaters];
529     double* waterY = new double[testWaters];
530     double* waterZ = new double[testWaters];
531    
532     double x0 = 0.0 - ( testBox * 0.5 );
533     double y0 = 0.0 - ( testBox * 0.5 );
534     double z0 = 0.0 - ( testBox * 0.5 );
535    
536    
537     // create an fcc lattice in the water box.
538    
539     int ndx = 0;
540     for( i=0; i < nCells; i++ ){
541     for( j=0; j < nCells; j++ ){
542     for( k=0; k < nCells; k++ ){
543    
544     myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
545     for(l=0; l<4; l++){
546     waterX[ndx]=posX[l];
547     waterY[ndx]=posY[l];
548     waterZ[ndx]=posZ[l];
549     ndx++;
550     }
551     }
552     }
553     }
554    
555     // calculate the number of water's displaced by our lipid.
556    
557     testSite.rot[0][0] = 1.0;
558     testSite.rot[0][1] = 0.0;
559     testSite.rot[0][2] = 0.0;
560    
561     testSite.rot[1][0] = 0.0;
562     testSite.rot[1][1] = 1.0;
563     testSite.rot[1][2] = 0.0;
564    
565     testSite.rot[2][0] = 0.0;
566     testSite.rot[2][1] = 0.0;
567     testSite.rot[2][2] = 1.0;
568    
569     testSite.pos[0] = 0.0;
570     testSite.pos[1] = 0.0;
571     testSite.pos[2] = 0.0;
572    
573     lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig );
574    
575     int *isActive = new int[testWaters];
576     for(i=0; i<testWaters; i++) isActive[i] = 1;
577    
578     int n_deleted = 0;
579     double dx, dy, dz;
580     double dx2, dy2, dz2, dSqr;
581     double rCutSqr = water_padding * water_padding;
582    
583     for(i=0; ( (i<testWaters) && isActive[i] ); i++){
584     for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
585    
586     atoms[j]->getPos( pos );
587    
588     dx = waterX[i] - pos[0];
589     dy = waterY[i] - pos[1];
590     dz = waterZ[i] - pos[2];
591    
592     buildMap( dx, dy, dz, testBox, testBox, testBox );
593    
594     dx2 = dx * dx;
595     dy2 = dy * dy;
596     dz2 = dz * dz;
597    
598     dSqr = dx2 + dy2 + dz2;
599     if( dSqr < rCutSqr ){
600     isActive[i] = 0;
601     n_deleted++;
602     }
603     }
604     }
605    
606     int targetWaters = nWaters + n_deleted * nLipids;
607    
608     targetWaters = (int) ( targetWaters * 1.2 );
609    
610     // find the best box size for the sim
611    
612     int nCellsX, nCellsY, nCellsZ;
613    
614     const double boxTargetX = 66.22752;
615     const double boxTargetY = 60.53088;
616    
617     // nCellsX = (int)ceil(boxTargetX / waterCell);
618     // nCellsY = (int)ceil(boxTargetY / waterCell);
619    
620     int testTot;
621     int done = 0;
622     nCellsX = 0;
623     nCellsY = 0;
624     nCellsZ = 0;
625     while( !done ){
626    
627     nCellsX++;
628     nCellsY++;
629     nCellsZ++;
630     testTot = 4 * nCellsX * nCellsY * nCellsZ;
631    
632     if( testTot >= targetWaters ) done = 1;
633     }
634    
635     // create the new water box to the new specifications
636    
637     int newWaters = nCellsX * nCellsY * nCellsZ * 4;
638    
639     delete[] waterX;
640     delete[] waterY;
641     delete[] waterZ;
642    
643     coord* waterSites = new coord[newWaters];
644    
645     double box_x = waterCell * nCellsX;
646     double box_y = waterCell * nCellsY;
647     double box_z = waterCell * nCellsZ;
648    
649     // create an fcc lattice in the water box.
650    
651     ndx = 0;
652     for( i=0; i < nCellsX; i++ ){
653     for( j=0; j < nCellsY; j++ ){
654     for( k=0; k < nCellsZ; k++ ){
655    
656     myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
657     for(l=0; l<4; l++){
658     waterSites[ndx].pos[0] = posX[l];
659     waterSites[ndx].pos[1] = posY[l];
660     waterSites[ndx].pos[2] = posZ[l];
661     ndx++;
662     }
663     }
664     }
665     }
666    
667     coord* lipidSites = new coord[nLipids];
668    
669     // start a 3D RSA for the for the lipid placements
670    
671    
672     int reject;
673     int testDX, acceptedDX;
674     SimInfo* testInfo2;
675    
676     nAtoms = nLipids * lipidNatoms;
677     testInfo2 = new SimInfo();
678     testInfo2->n_atoms = nAtoms;
679     theConfig = testInfo2->getConfiguration();
680     theConfig->createArrays( nAtoms );
681     testInfo2->atoms = new Atom*[nAtoms];
682 tim 1417 for (i=0; i < nAtoms; i++) {
683     testInfo2->atoms[i] = new Atom(i, theConfig);
684     testInfo2->atoms[i]->setCoords();
685     }
686 gezelter 1334 atoms = testInfo2->atoms;
687    
688     rCutSqr = lipid_spaceing * lipid_spaceing;
689    
690     for(i=0; i<nLipids; i++ ){
691     done = 0;
692     while( !done ){
693    
694     lipidSites[i].pos[0] = drand48() * box_x;
695     lipidSites[i].pos[1] = drand48() * box_y;
696     lipidSites[i].pos[2] = drand48() * box_z;
697    
698     getRandomRot( lipidSites[i].rot );
699    
700     ndx = i * lipidNatoms;
701    
702     lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
703     ndx, theConfig );
704    
705     reject = 0;
706     for( j=0; !reject && j<i; j++){
707     for(k=0; !reject && k<lipidNatoms; k++){
708    
709     acceptedDX = j*lipidNatoms + k;
710     for(l=0; !reject && l<lipidNatoms; l++){
711    
712     testDX = ndx + l;
713    
714     atoms[testDX]->getPos( posA );
715     atoms[acceptedDX]->getPos( posB );
716    
717     dx = posA[0] - posB[0];
718     dy = posA[1] - posB[1];
719     dz = posA[2] - posB[2];
720    
721     buildMap( dx, dy, dz, box_x, box_y, box_z );
722    
723     dx2 = dx * dx;
724     dy2 = dy * dy;
725     dz2 = dz * dz;
726    
727     dSqr = dx2 + dy2 + dz2;
728     if( dSqr < rCutSqr ) reject = 1;
729     }
730     }
731     }
732    
733     if( reject ){
734    
735     for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j];
736     }
737     else{
738     done = 1;
739     std::cout << (i+1) << " has been accepted\n";
740     std::cout.flush();
741     }
742     }
743     }
744    
745    
746     // // zSort of the lipid positions
747    
748    
749     // vector< pair<int,double> >zSortArray;
750     // for(i=0;i<nLipids;i++)
751     // zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) );
752    
753     // sort(zSortArray.begin(),zSortArray.end(),SortCond());
754    
755     // ofstream outFile( "./zipper.bass", ios::app);
756    
757     // for(i=0; i<nLipids; i++){
758     // outFile << "zConstraint[" << i << "]{\n"
759     // << " molIndex = " << zSortArray[i].first << ";\n"
760     // << " zPos = ";
761    
762     // if(i<32) outFile << "60.0;\n";
763     // else outFile << "100.0;\n";
764    
765     // outFile << " kRatio = 0.5;\n"
766     // << "}\n";
767     // }
768    
769     // outFile.close();
770    
771    
772     // cut out the waters that overlap with the lipids.
773    
774    
775     delete[] isActive;
776     isActive = new int[newWaters];
777     for(i=0; i<newWaters; i++) isActive[i] = 1;
778     int n_active = newWaters;
779     rCutSqr = water_padding * water_padding;
780    
781     for(i=0; ( (i<newWaters) && isActive[i] ); i++){
782     for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
783    
784     atoms[j]->getPos( pos );
785    
786     dx = waterSites[i].pos[0] - pos[0];
787     dy = waterSites[i].pos[1] - pos[1];
788     dz = waterSites[i].pos[2] - pos[2];
789    
790     buildMap( dx, dy, dz, box_x, box_y, box_z );
791    
792     dx2 = dx * dx;
793     dy2 = dy * dy;
794     dz2 = dz * dz;
795    
796     dSqr = dx2 + dy2 + dz2;
797     if( dSqr < rCutSqr ){
798     isActive[i] = 0;
799     n_active--;
800    
801    
802     }
803     }
804     }
805    
806    
807    
808    
809     if( n_active < nWaters ){
810    
811     sprintf( painCave.errMsg,
812     "Too many waters were removed, edit code and try again.\n" );
813    
814     painCave.isFatal = 1;
815     simError();
816     }
817    
818     int quickKill;
819     while( n_active > nWaters ){
820    
821     quickKill = (int)(drand48()*newWaters);
822    
823     if( isActive[quickKill] ){
824     isActive[quickKill] = 0;
825     n_active--;
826    
827     }
828     }
829    
830     if( n_active != nWaters ){
831    
832     sprintf( painCave.errMsg,
833     "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
834     n_active, nWaters );
835     painCave.isFatal = 1;
836     simError();
837     }
838    
839     // clean up our messes before building the final system.
840    
841     testInfo->getConfiguration()->destroyArrays();
842     testInfo2->getConfiguration()->destroyArrays();
843    
844     // create the real Atom arrays
845    
846     nAtoms = 0;
847     molIndex = 0;
848     molStart = new int[nLipids + nWaters];
849    
850     for(j=0; j<nLipids; j++){
851     molStart[molIndex] = nAtoms;
852     molIndex++;
853     nAtoms += lipidNatoms;
854     }
855    
856     for(j=0; j<nWaters; j++){
857     molStart[molIndex] = nAtoms;
858     molIndex++;
859     nAtoms += waterNatoms;
860     }
861    
862 tim 1417 std::cerr << "Flag-1\n";
863 gezelter 1334 theConfig = mainInfo->getConfiguration();
864 tim 1417 mainInfo->n_atoms = nAtoms;
865 gezelter 1334 theConfig->createArrays( nAtoms );
866     mainInfo->atoms = new Atom*[nAtoms];
867 tim 1417 for (i=0; i < nAtoms; i++) {
868     mainInfo->atoms[i] = new Atom(i, theConfig);
869     mainInfo->atoms[i]->setCoords();
870     }
871 gezelter 1334 atoms = mainInfo->atoms;
872 tim 1417 std::cerr << "Flag0\n";
873 gezelter 1334
874     // set up the SimInfo object
875    
876     double Hmat[3][3];
877    
878     Hmat[0][0] = box_x;
879     Hmat[0][1] = 0.0;
880     Hmat[0][2] = 0.0;
881    
882     Hmat[1][0] = 0.0;
883     Hmat[1][1] = box_y;
884     Hmat[1][2] = 0.0;
885    
886     Hmat[2][0] = 0.0;
887     Hmat[2][1] = 0.0;
888     Hmat[2][2] = box_z;
889    
890     mainInfo->setBoxM( Hmat );
891    
892     // center the system on (0,0,0)
893    
894     for(j=0;j<nLipids;j++){
895    
896     mainInfo->wrapVector( lipidSites[j].pos );
897     }
898    
899     for(j=0;j<newWaters;j++){
900    
901     mainInfo->wrapVector( waterSites[j].pos );
902     }
903    
904     // initialize lipid positions
905    
906 tim 1417 std::cerr << "Flag1\n";
907    
908 gezelter 1334 molIndex = 0;
909     for(i=0; i<nLipids; i++ ){
910     lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
911     molStart[molIndex], theConfig );
912     molIndex++;
913     }
914    
915     // initialize the water positions
916 tim 1417 std::cerr << "Flag2\n";
917 gezelter 1334
918     for(i=0; i<newWaters; i++){
919    
920     if( isActive[i] ){
921    
922     getRandomRot( waterSites[i].rot );
923     waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
924     molStart[molIndex], theConfig );
925     molIndex++;
926     }
927     }
928    
929 tim 1417 std::cerr << "Flag3\n";
930 gezelter 1334
931    
932     // set up the writer and write out
933    
934     writer = new DumpWriter( mainInfo );
935 tim 1417 std::cerr << "Flag4\n";
936 gezelter 1334 writer->writeFinal( 0.0 );
937 tim 1417 std::cerr << "Flag5\n";
938 gezelter 1334
939     return 1;
940     }
941    
942     void buildMap( double &x, double &y, double &z,
943     double boxX, double boxY, double boxZ ){
944    
945     if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
946     else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
947    
948     if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
949     else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
950    
951     if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
952     else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
953     }
954    
955     /***************************************************************************
956     * prints out the usage for the command line arguments, then exits.
957     ***************************************************************************/
958    
959     void usage(){
960     (void)fprintf(stdout,
961     "\n"
962     "The proper usage is: %s [options] <input_file>\n"
963     "\n"
964     "Options:\n"
965     "\n"
966     " short:\n"
967     " ------\n"
968     " -o <name> The output prefix\n"
969     " -l <lipidName> The name of the lipid molecule specified in the BASS file.\n"
970     " -w <waterName> The name of the water molecule specified in the BASS file.\n"
971    
972     "\n"
973     " long:\n"
974     " -----\n"
975    
976     " --version displays the version number\n"
977     " --help displays this help message.\n"
978     " --waterBuffer <#> sets the distance of closest approach of the water around the lipid\n"
979     " defaults to 5.0\n"
980     " --lipidBuffer <#> sets the distance of closest approach between two lipids\n"
981     " defaults to 6.0\n"
982     " --waterLattice <#> sets the water lattice spacing\n"
983     " defaults to 4.929 ( 1 g/cm^3 )\n"
984     "\n"
985     "\n",
986     programName);
987     }