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 |
|
|
} |