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 (19 years, 11 months ago) by tim
File size: 21405 byte(s)
Log Message:
BASS eradication project (part 1)

File Contents

# Content
1 #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 for (i=0; i < lipidNatoms; i++) {
517 testInfo->atoms[i] = new Atom(i, theConfig);
518 testInfo->atoms[i]->setCoords();
519 }
520 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 for (i=0; i < nAtoms; i++) {
683 testInfo2->atoms[i] = new Atom(i, theConfig);
684 testInfo2->atoms[i]->setCoords();
685 }
686 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 std::cerr << "Flag-1\n";
863 theConfig = mainInfo->getConfiguration();
864 mainInfo->n_atoms = nAtoms;
865 theConfig->createArrays( nAtoms );
866 mainInfo->atoms = new Atom*[nAtoms];
867 for (i=0; i < nAtoms; i++) {
868 mainInfo->atoms[i] = new Atom(i, theConfig);
869 mainInfo->atoms[i]->setCoords();
870 }
871 atoms = mainInfo->atoms;
872 std::cerr << "Flag0\n";
873
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 std::cerr << "Flag1\n";
907
908 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 std::cerr << "Flag2\n";
917
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 std::cerr << "Flag3\n";
930
931
932 // set up the writer and write out
933
934 writer = new DumpWriter( mainInfo );
935 std::cerr << "Flag4\n";
936 writer->writeFinal( 0.0 );
937 std::cerr << "Flag5\n";
938
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 }