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 |
atoms = testInfo->atoms; |
517 |
|
518 |
// create the test box for initial water displacement |
519 |
|
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 |
|
524 |
double* waterX = new double[testWaters]; |
525 |
double* waterY = new double[testWaters]; |
526 |
double* waterZ = new double[testWaters]; |
527 |
|
528 |
double x0 = 0.0 - ( testBox * 0.5 ); |
529 |
double y0 = 0.0 - ( testBox * 0.5 ); |
530 |
double z0 = 0.0 - ( testBox * 0.5 ); |
531 |
|
532 |
|
533 |
// create an fcc lattice in the water box. |
534 |
|
535 |
int ndx = 0; |
536 |
for( i=0; i < nCells; i++ ){ |
537 |
for( j=0; j < nCells; j++ ){ |
538 |
for( k=0; k < nCells; k++ ){ |
539 |
|
540 |
myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k); |
541 |
for(l=0; l<4; l++){ |
542 |
waterX[ndx]=posX[l]; |
543 |
waterY[ndx]=posY[l]; |
544 |
waterZ[ndx]=posZ[l]; |
545 |
ndx++; |
546 |
} |
547 |
} |
548 |
} |
549 |
} |
550 |
|
551 |
// calculate the number of water's displaced by our lipid. |
552 |
|
553 |
testSite.rot[0][0] = 1.0; |
554 |
testSite.rot[0][1] = 0.0; |
555 |
testSite.rot[0][2] = 0.0; |
556 |
|
557 |
testSite.rot[1][0] = 0.0; |
558 |
testSite.rot[1][1] = 1.0; |
559 |
testSite.rot[1][2] = 0.0; |
560 |
|
561 |
testSite.rot[2][0] = 0.0; |
562 |
testSite.rot[2][1] = 0.0; |
563 |
testSite.rot[2][2] = 1.0; |
564 |
|
565 |
testSite.pos[0] = 0.0; |
566 |
testSite.pos[1] = 0.0; |
567 |
testSite.pos[2] = 0.0; |
568 |
|
569 |
lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig ); |
570 |
|
571 |
int *isActive = new int[testWaters]; |
572 |
for(i=0; i<testWaters; i++) isActive[i] = 1; |
573 |
|
574 |
int n_deleted = 0; |
575 |
double dx, dy, dz; |
576 |
double dx2, dy2, dz2, dSqr; |
577 |
double rCutSqr = water_padding * water_padding; |
578 |
|
579 |
for(i=0; ( (i<testWaters) && isActive[i] ); i++){ |
580 |
for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){ |
581 |
|
582 |
atoms[j]->getPos( pos ); |
583 |
|
584 |
dx = waterX[i] - pos[0]; |
585 |
dy = waterY[i] - pos[1]; |
586 |
dz = waterZ[i] - pos[2]; |
587 |
|
588 |
buildMap( dx, dy, dz, testBox, testBox, testBox ); |
589 |
|
590 |
dx2 = dx * dx; |
591 |
dy2 = dy * dy; |
592 |
dz2 = dz * dz; |
593 |
|
594 |
dSqr = dx2 + dy2 + dz2; |
595 |
if( dSqr < rCutSqr ){ |
596 |
isActive[i] = 0; |
597 |
n_deleted++; |
598 |
} |
599 |
} |
600 |
} |
601 |
|
602 |
int targetWaters = nWaters + n_deleted * nLipids; |
603 |
|
604 |
targetWaters = (int) ( targetWaters * 1.2 ); |
605 |
|
606 |
// find the best box size for the sim |
607 |
|
608 |
int nCellsX, nCellsY, nCellsZ; |
609 |
|
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); |
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 |
|
628 |
if( testTot >= targetWaters ) done = 1; |
629 |
} |
630 |
|
631 |
// create the new water box to the new specifications |
632 |
|
633 |
int newWaters = nCellsX * nCellsY * nCellsZ * 4; |
634 |
|
635 |
delete[] waterX; |
636 |
delete[] waterY; |
637 |
delete[] waterZ; |
638 |
|
639 |
coord* waterSites = new coord[newWaters]; |
640 |
|
641 |
double box_x = waterCell * nCellsX; |
642 |
double box_y = waterCell * nCellsY; |
643 |
double box_z = waterCell * nCellsZ; |
644 |
|
645 |
// create an fcc lattice in the water box. |
646 |
|
647 |
ndx = 0; |
648 |
for( i=0; i < nCellsX; i++ ){ |
649 |
for( j=0; j < nCellsY; j++ ){ |
650 |
for( k=0; k < nCellsZ; k++ ){ |
651 |
|
652 |
myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k); |
653 |
for(l=0; l<4; l++){ |
654 |
waterSites[ndx].pos[0] = posX[l]; |
655 |
waterSites[ndx].pos[1] = posY[l]; |
656 |
waterSites[ndx].pos[2] = posZ[l]; |
657 |
ndx++; |
658 |
} |
659 |
} |
660 |
} |
661 |
} |
662 |
|
663 |
coord* lipidSites = new coord[nLipids]; |
664 |
|
665 |
// start a 3D RSA for the for the lipid placements |
666 |
|
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 |
|
680 |
rCutSqr = lipid_spaceing * lipid_spaceing; |
681 |
|
682 |
for(i=0; i<nLipids; i++ ){ |
683 |
done = 0; |
684 |
while( !done ){ |
685 |
|
686 |
lipidSites[i].pos[0] = drand48() * box_x; |
687 |
lipidSites[i].pos[1] = drand48() * box_y; |
688 |
lipidSites[i].pos[2] = drand48() * box_z; |
689 |
|
690 |
getRandomRot( lipidSites[i].rot ); |
691 |
|
692 |
ndx = i * lipidNatoms; |
693 |
|
694 |
lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms, |
695 |
ndx, theConfig ); |
696 |
|
697 |
reject = 0; |
698 |
for( j=0; !reject && j<i; j++){ |
699 |
for(k=0; !reject && k<lipidNatoms; k++){ |
700 |
|
701 |
acceptedDX = j*lipidNatoms + k; |
702 |
for(l=0; !reject && l<lipidNatoms; l++){ |
703 |
|
704 |
testDX = ndx + l; |
705 |
|
706 |
atoms[testDX]->getPos( posA ); |
707 |
atoms[acceptedDX]->getPos( posB ); |
708 |
|
709 |
dx = posA[0] - posB[0]; |
710 |
dy = posA[1] - posB[1]; |
711 |
dz = posA[2] - posB[2]; |
712 |
|
713 |
buildMap( dx, dy, dz, box_x, box_y, box_z ); |
714 |
|
715 |
dx2 = dx * dx; |
716 |
dy2 = dy * dy; |
717 |
dz2 = dz * dz; |
718 |
|
719 |
dSqr = dx2 + dy2 + dz2; |
720 |
if( dSqr < rCutSqr ) reject = 1; |
721 |
} |
722 |
} |
723 |
} |
724 |
|
725 |
if( reject ){ |
726 |
|
727 |
for(j=0; j< lipidNatoms; j++) delete atoms[ndx+j]; |
728 |
} |
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 |
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]) ); |
744 |
|
745 |
// sort(zSortArray.begin(),zSortArray.end(),SortCond()); |
746 |
|
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 = "; |
753 |
|
754 |
// if(i<32) outFile << "60.0;\n"; |
755 |
// else outFile << "100.0;\n"; |
756 |
|
757 |
// outFile << " kRatio = 0.5;\n" |
758 |
// << "}\n"; |
759 |
// } |
760 |
|
761 |
// outFile.close(); |
762 |
|
763 |
|
764 |
// cut out the waters that overlap with the lipids. |
765 |
|
766 |
|
767 |
delete[] isActive; |
768 |
isActive = new int[newWaters]; |
769 |
for(i=0; i<newWaters; i++) isActive[i] = 1; |
770 |
int n_active = newWaters; |
771 |
rCutSqr = water_padding * water_padding; |
772 |
|
773 |
for(i=0; ( (i<newWaters) && isActive[i] ); i++){ |
774 |
for(j=0; ( (j<nAtoms) && isActive[i] ); j++){ |
775 |
|
776 |
atoms[j]->getPos( pos ); |
777 |
|
778 |
dx = waterSites[i].pos[0] - pos[0]; |
779 |
dy = waterSites[i].pos[1] - pos[1]; |
780 |
dz = waterSites[i].pos[2] - pos[2]; |
781 |
|
782 |
buildMap( dx, dy, dz, box_x, box_y, box_z ); |
783 |
|
784 |
dx2 = dx * dx; |
785 |
dy2 = dy * dy; |
786 |
dz2 = dz * dz; |
787 |
|
788 |
dSqr = dx2 + dy2 + dz2; |
789 |
if( dSqr < rCutSqr ){ |
790 |
isActive[i] = 0; |
791 |
n_active--; |
792 |
|
793 |
|
794 |
} |
795 |
} |
796 |
} |
797 |
|
798 |
|
799 |
|
800 |
|
801 |
if( n_active < nWaters ){ |
802 |
|
803 |
sprintf( painCave.errMsg, |
804 |
"Too many waters were removed, edit code and try again.\n" ); |
805 |
|
806 |
painCave.isFatal = 1; |
807 |
simError(); |
808 |
} |
809 |
|
810 |
int quickKill; |
811 |
while( n_active > nWaters ){ |
812 |
|
813 |
quickKill = (int)(drand48()*newWaters); |
814 |
|
815 |
if( isActive[quickKill] ){ |
816 |
isActive[quickKill] = 0; |
817 |
n_active--; |
818 |
|
819 |
} |
820 |
} |
821 |
|
822 |
if( n_active != nWaters ){ |
823 |
|
824 |
sprintf( painCave.errMsg, |
825 |
"QuickKill didn't work right. n_active = %d, and nWaters = %d\n", |
826 |
n_active, nWaters ); |
827 |
painCave.isFatal = 1; |
828 |
simError(); |
829 |
} |
830 |
|
831 |
// clean up our messes before building the final system. |
832 |
|
833 |
testInfo->getConfiguration()->destroyArrays(); |
834 |
testInfo2->getConfiguration()->destroyArrays(); |
835 |
|
836 |
// create the real Atom arrays |
837 |
|
838 |
nAtoms = 0; |
839 |
molIndex = 0; |
840 |
molStart = new int[nLipids + nWaters]; |
841 |
|
842 |
for(j=0; j<nLipids; j++){ |
843 |
molStart[molIndex] = nAtoms; |
844 |
molIndex++; |
845 |
nAtoms += lipidNatoms; |
846 |
} |
847 |
|
848 |
for(j=0; j<nWaters; j++){ |
849 |
molStart[molIndex] = nAtoms; |
850 |
molIndex++; |
851 |
nAtoms += waterNatoms; |
852 |
} |
853 |
|
854 |
theConfig = mainInfo->getConfiguration(); |
855 |
theConfig->createArrays( nAtoms ); |
856 |
mainInfo->atoms = new Atom*[nAtoms]; |
857 |
atoms = mainInfo->atoms; |
858 |
mainInfo->n_atoms = nAtoms; |
859 |
|
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, |
895 |
molStart[molIndex], theConfig ); |
896 |
molIndex++; |
897 |
} |
898 |
|
899 |
// initialize the water positions |
900 |
|
901 |
for(i=0; i<newWaters; i++){ |
902 |
|
903 |
if( isActive[i] ){ |
904 |
|
905 |
getRandomRot( waterSites[i].rot ); |
906 |
waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms, |
907 |
molStart[molIndex], theConfig ); |
908 |
molIndex++; |
909 |
} |
910 |
} |
911 |
|
912 |
|
913 |
|
914 |
// set up the writer and write out |
915 |
|
916 |
writer = new DumpWriter( mainInfo ); |
917 |
writer->writeFinal( 0.0 ); |
918 |
|
919 |
return 1; |
920 |
} |
921 |
|
922 |
void buildMap( double &x, double &y, double &z, |
923 |
double boxX, double boxY, double boxZ ){ |
924 |
|
925 |
if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) ); |
926 |
else x -= boxX * (double)( (int)( (x / boxX ) + 0.5)); |
927 |
|
928 |
if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) ); |
929 |
else y -= boxY * (double)( (int)( (y / boxY ) + 0.5)); |
930 |
|
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 |
} |