ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/bilayerSys.cpp
Revision: 536
Committed: Fri May 16 14:28:27 2003 UTC (21 years, 11 months ago) by mmeineke
File size: 12645 byte(s)
Log Message:
doing some work to overhaul sysbuild.

File Contents

# Content
1 #include <iostream>
2
3 #include <cstdlib>
4 #include <cstring>
5 #include <cmath>
6
7 #include "simError.h"
8 #include "SimInfo.hpp"
9 #include "ReadWrite.hpp"
10
11 #include "MoLocator.hpp"
12 #include "sysBuild.hpp"
13 #include "bilayerSys.hpp"
14
15
16
17 void map( double &x, double &y, double &z,
18 double boxX, double boxY, double boxZ );
19
20 int buildRandomBilayer( void );
21
22 void getRandomRot( double rot[3][3] );
23
24 int buildBilayer( int isRandom ){
25
26 if( isRandom ){
27 return buildRandomBilayer();
28 }
29 else{
30 sprintf( painCave.errMsg,
31 "Cannot currently create a non-random bilayer.\n" );
32 painCave.isFatal = 1;
33 simError();
34 return 0;
35 }
36 }
37
38
39 int buildRandomBilayer( void ){
40
41 typedef struct{
42 double rot[3][3];
43 double pos[3];
44 } coord;
45
46
47 const double waterRho = 0.0334; // number density per cubic angstrom
48 const double waterVol = 4.0 / water_rho; // volume occupied by 4 waters
49 const double waterCell = 4.929; // fcc unit cell length
50
51 const double water_padding = 2.5;
52 const double lipid_spaceing = 2.5;
53
54
55 int i,j,k;
56 int nAtoms, atomIndex, molIndex, molID;
57 int* molSeq;
58 int* molMap;
59 int* molStart;
60 int* cardDeck;
61 int deckSize;
62 int rSite, rCard;
63 double cell;
64 int nCells, nSites, siteIndex;
65
66 coord *siteArray;
67 coord testSite;
68
69 MoleculeStamp* lipidStamp;
70 MoleculeStamp* waterStamp;
71 MoLocator *lipidLocate;
72 MoLocator *waterLocate
73 int foundLipid, foundWater;
74 int nLipids, lipiNatoms, nWaters;
75 double testBox, maxLength;
76
77 srand48( RAND_SEED );
78
79
80 // set the the lipidStamp
81
82 foundLipid = 0;
83 for(i=0; i<bsInfo.nComponents; i++){
84 if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
85
86 foundlipid = 1;
87 lipidStamp = bsInfo.compStamps[i];
88 nLipids = bsInfo.componentsNmol[i];
89 }
90 if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
91
92 foundWater = 1;
93 waterStamp = bsInfo.compStamps[i];
94 nWaters = bsInfo.componentsNmol[i];
95 }
96 }
97 if( !foundLipid ){
98 sprintf(painCave.errMsg,
99 "Could not find lipid \"%s\" in the bass file.\n",
100 bsInfo.lipidName );
101 painCave.isFatal = 1;
102 simError();
103 }
104 if( !foundWater ){
105 sprintf(painCave.errMsg,
106 "Could not find water \"%s\" in the bass file.\n",
107 bsInfo.waterName );
108 painCave.isFatal = 1;
109 simError();
110 }
111
112 //create the temp Molocator and atom Arrays
113
114 lipidLocate = new MoLocator( lipidStamp );
115 lipidNatoms = lipidStamp->getNAtoms();
116 maxLength = lipidLocate->getMaxLength();
117
118 waterLocate = new MoLocator( waterStamp );
119
120 nAtoms = nLipids * lipidNatoms;
121
122 Atom::createArrays( nAtoms );
123 atoms = new Atom*[nAtoms];
124
125 // create the test box for initial water displacement
126
127 testBox = maxLength + waterCell * 4.0; // pad with 4 cells
128 int nCells = (int)( testBox / waterCell + 1.0 );
129 int testWaters = 4 * nCells * nCells * nCells;
130
131 double* waterX = new double[testWaters];
132 double* waterX = new double[testWaters];
133 double* waterX = new double[testWaters];
134
135 double x0 = 0.0 - ( testBox * 0.5 );
136 double y0 = 0.0 - ( testBox * 0.5 );
137 double z0 = 0.0 - ( testBox * 0.5 );
138
139
140 // create an fcc lattice in the water box.
141
142 int ndx = 0;
143 for( i=0; i < nCells; i++ ){
144 for( j=0; j < nCells; j++ ){
145 for( k=0; k < nCells; k++ ){
146
147 waterX[ndx] = i * waterCell + x0;
148 waterY[ndx] = j * waterCell + y0;
149 waterZ[ndx] = k * waterCell + z0;
150 ndx++;
151
152 waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
153 waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
154 waterZ[ndx] = k * waterCell + z0;
155 ndx++;
156
157 waterX[ndx] = i * waterCell + x0;
158 waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
159 waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
160 ndx++;
161
162 waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
163 waterY[ndx] = j * waterCell + y0;
164 waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
165 ndx++;
166 }
167 }
168 }
169
170 // calculate the number of water's displaced by our lipid.
171
172 testSite.rot[0][0] = 1.0;
173 testSite.rot[0][1] = 0.0;
174 testSite.rot[0][2] = 0.0;
175
176 testSite.rot[1][0] = 0.0;
177 testSite.rot[1][1] = 1.0;
178 testSite.rot[1][2] = 0.0;
179
180 testSite.rot[2][0] = 0.0;
181 testSite.rot[2][1] = 0.0;
182 testSite.rot[2][2] = 1.0;
183
184 testSite.pos[0] = 0.0;
185 testSite.pos[1] = 0.0;
186 testSite.pos[2] = 0.0;
187
188 lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
189
190 int *isActive = new int[testWaters];
191 for(i=0; i<testWaters; i++) isActive[i] = 1;
192
193 int n_deleted = 0;
194 double dx, dy, dz;
195 double dx2, dy2, dz2, dSqr;
196 double rCutSqr = water_padding * water_padding;
197
198 for(i=0; ( (i<testWaters) && isActive[i] ); i++){
199 for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
200
201 dx = waterX[i] - atoms[j]->getX();
202 dy = waterY[i] - atoms[j]->getY();
203 dz = waterZ[i] - atoms[j]->getZ();
204
205 map( dx, dy, dz, testBox, testBox, testBox );
206
207 dx2 = dx * dx;
208 dy2 = dy * dy;
209 dz2 = dz * dz;
210
211 dSqr = dx2 + dy2 + dz2;
212 if( dSqr < rCutSqr ){
213 isActive[i] = 0;
214 n_deleted++;
215 }
216 }
217 }
218
219 int targetWaters = nWaters + n_deleted * nLipids;
220
221 // find the best box size for the sim
222
223 int testTot;
224 int done = 0;
225 ndx = 0;
226 while( !done ){
227
228 ndx++;
229 testTot = 4 * ndx * ndx * ndx;
230
231 if( testTot >= targetWaters ) done = 1;
232 }
233
234 nCells = ndx;
235
236
237 // create the new water box to the new specifications
238
239 int newWaters = nCells * nCells * nCells * 4;
240
241 delete[] waterX;
242 delete[] waterY;
243 delete[] waterZ;
244
245 waterX = new double[newWater];
246 waterY = new double[newWater];
247 waterZ = new double[newWater];
248
249 double box_x = waterCell * nCells;
250 double box_y = waterCell * nCells;
251 double box_z = waterCell * nCells;
252
253 // create an fcc lattice in the water box.
254
255 ndx = 0;
256 for( i=0; i < nCells; i++ ){
257 for( j=0; j < nCells; j++ ){
258 for( k=0; k < nCells; k++ ){
259
260 waterX[ndx] = i * waterCell;
261 waterY[ndx] = j * waterCell;
262 waterZ[ndx] = k * waterCell;
263 ndx++;
264
265 waterX[ndx] = i * waterCell + 0.5 * waterCell;
266 waterY[ndx] = j * waterCell + 0.5 * waterCell;
267 waterZ[ndx] = k * waterCell;
268 ndx++;
269
270 waterX[ndx] = i * waterCell;
271 waterY[ndx] = j * waterCell + 0.5 * waterCell;
272 waterZ[ndx] = k * waterCell + 0.5 * waterCell;
273 ndx++;
274
275 waterX[ndx] = i * waterCell + 0.5 * waterCell;
276 waterY[ndx] = j * waterCell;
277 waterZ[ndx] = k * waterCell + 0.5 * waterCell;
278 ndx++;
279 }
280 }
281 }
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296 // create the real MoLocator and Atom arrays
297
298 nAtoms = 0;
299 molIndex = 0;
300 locate = new MoLocator*[bsInfo.nComponents];
301 molSeq = new int[bsInfo.totNmol];
302 molStart = new int[bsInfo.totNmol];
303 for(i=0; i<bsInfo.nComponents; i++){
304 locate[i] = new MoLocator( bsInfo.compStamps[i] );
305 for(j=0; j<bsInfo.componentsNmol[i]; j++){
306 molSeq[molIndex] = i;
307 molStart[molIndex] = nAtoms;
308 molIndex++;
309 nAtoms += bsInfo.compStamps[i]->getNAtoms();
310 }
311 }
312
313 Atom::createArrays( nAtoms );
314 atoms = new Atom*[nAtoms];
315
316
317 // find the width, height, and length of the molecule
318
319
320
321
322 }
323
324
325
326 int Old_buildRandomBilayer( void ){
327
328 int i,j,k;
329 int nAtoms, atomIndex, molIndex, molID;
330 int* molSeq;
331 int* molMap;
332 int* molStart;
333 int* cardDeck;
334 int deckSize;
335 int rSite, rCard;
336 double cell;
337 int nCells, nSites, siteIndex;
338 double rot[3][3];
339 double pos[3];
340
341 Atom** atoms;
342 SimInfo* simnfo;
343 DumpWriter* writer;
344 MoLocator** locate;
345
346 // initialize functions and variables
347
348 srand48( RAND_SEED );
349 molSeq = NULL;
350 molStart = NULL;
351 molMap = NULL;
352 cardDeck = NULL;
353 atoms = NULL;
354 locate = NULL;
355 simnfo = NULL;
356 writer = NULL;
357
358 // calculate the number of cells in the fcc box
359
360 nCells = 0;
361 nSites = 0;
362 while( nSites < bsInfo.totNmol ){
363 nCells++;
364 nSites = 4.0 * pow( (double)nCells, 3.0 );
365 }
366
367
368 // create the molMap and cardDeck arrays
369
370 molMap = new int[nSites];
371 cardDeck = new int[nSites];
372
373 for(i=0; i<nSites; i++){
374 molMap[i] = -1;
375 cardDeck[i] = i;
376 }
377
378 // randomly place the molecules on the sites
379
380 deckSize = nSites;
381 for(i=0; i<bsInfo.totNmol; i++){
382 rCard = (int)( deckSize * drand48() );
383 rSite = cardDeck[rCard];
384 molMap[rSite] = i;
385
386 // book keep the card deck;
387
388 deckSize--;
389 cardDeck[rCard] = cardDeck[deckSize];
390 }
391
392
393 // create the MoLocator and Atom arrays
394
395 nAtoms = 0;
396 molIndex = 0;
397 locate = new MoLocator*[bsInfo.nComponents];
398 molSeq = new int[bsInfo.totNmol];
399 molStart = new int[bsInfo.totNmol];
400 for(i=0; i<bsInfo.nComponents; i++){
401 locate[i] = new MoLocator( bsInfo.compStamps[i] );
402 for(j=0; j<bsInfo.componentsNmol[i]; j++){
403 molSeq[molIndex] = i;
404 molStart[molIndex] = nAtoms;
405 molIndex++;
406 nAtoms += bsInfo.compStamps[i]->getNAtoms();
407 }
408 }
409
410 Atom::createArrays( nAtoms );
411 atoms = new Atom*[nAtoms];
412
413
414 // place the molecules at each FCC site
415
416 cell = 5.0;
417 for(i=0; i<bsInfo.nComponents; i++){
418 if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
419 }
420 cell *= 1.2; // add a little buffer
421
422 cell *= M_SQRT2;
423
424 siteIndex = 0;
425 for(i=0; i<nCells; i++){
426 for(j=0; j<nCells; j++){
427 for(k=0; k<nCells; k++){
428
429 if( molMap[siteIndex] >= 0 ){
430 pos[0] = i * cell;
431 pos[1] = j * cell;
432 pos[2] = k * cell;
433
434 getRandomRot( rot );
435 molID = molSeq[molMap[siteIndex]];
436 atomIndex = molStart[ molMap[siteIndex] ];
437 locate[molID]->placeMol( pos, rot, atoms, atomIndex );
438 }
439 siteIndex++;
440
441 if( molMap[siteIndex] >= 0 ){
442 pos[0] = i * cell + (0.5 * cell);
443 pos[1] = j * cell;
444 pos[2] = k * cell + (0.5 * cell);
445
446 getRandomRot( rot );
447 molID = molSeq[molMap[siteIndex]];
448 atomIndex = molStart[ molMap[siteIndex] ];
449 locate[molID]->placeMol( pos, rot, atoms, atomIndex );
450 }
451 siteIndex++;
452
453 if( molMap[siteIndex] >= 0 ){
454 pos[0] = i * cell + (0.5 * cell);
455 pos[1] = j * cell + (0.5 * cell);
456 pos[2] = k * cell;
457
458 getRandomRot( rot );
459 molID = molSeq[molMap[siteIndex]];
460 atomIndex = molStart[ molMap[siteIndex] ];
461 locate[molID]->placeMol( pos, rot, atoms, atomIndex );
462 }
463 siteIndex++;
464
465 if( molMap[siteIndex] >= 0 ){
466 pos[0] = i * cell;
467 pos[1] = j * cell + (0.5 * cell);
468 pos[2] = k * cell + (0.5 * cell);
469
470 getRandomRot( rot );
471 molID = molSeq[molMap[siteIndex]];
472 atomIndex = molStart[ molMap[siteIndex] ];
473 locate[molID]->placeMol( pos, rot, atoms, atomIndex );
474 }
475 siteIndex++;
476 }
477 }
478 }
479
480 // set up the SimInfo object
481
482 bsInfo.boxX = nCells * cell;
483 bsInfo.boxY = nCells * cell;
484 bsInfo.boxZ = nCells * cell;
485
486 simnfo = new SimInfo();
487 simnfo->n_atoms = nAtoms;
488 simnfo->box_x = bsInfo.boxX;
489 simnfo->box_y = bsInfo.boxY;
490 simnfo->box_z = bsInfo.boxZ;
491
492 sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
493 sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
494
495 simnfo->atoms = atoms;
496
497 // set up the writer and write out
498
499 writer = new DumpWriter( simnfo );
500 writer->writeFinal();
501
502 // clean up the memory
503
504 if( molMap != NULL ) delete[] molMap;
505 if( cardDeck != NULL ) delete[] cardDeck;
506 if( locate != NULL ){
507 for(i=0; i<bsInfo.nComponents; i++){
508 delete locate[i];
509 }
510 delete[] locate;
511 }
512 if( atoms != NULL ){
513 for(i=0; i<nAtoms; i++){
514 delete atoms[i];
515 }
516 Atom::destroyArrays();
517 delete[] atoms;
518 }
519 if( molSeq != NULL ) delete[] molSeq;
520 if( simnfo != NULL ) delete simnfo;
521 if( writer != NULL ) delete writer;
522
523 return 1;
524 }
525
526
527 void getRandomRot( double rot[3][3] ){
528
529 double theta, phi, psi;
530 double cosTheta;
531
532 // select random phi, psi, and cosTheta
533
534 phi = 2.0 * M_PI * drand48();
535 psi = 2.0 * M_PI * drand48();
536 cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
537
538 theta = acos( cosTheta );
539
540 rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
541 rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
542 rot[0][2] = sin(theta) * sin(psi);
543
544 rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
545 rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
546 rot[1][2] = sin(theta) * cos(psi);
547
548 rot[2][0] = sin(phi) * sin(theta);
549 rot[2][1] = -cos(phi) * sin(theta);
550 rot[2][2] = cos(theta);
551 }
552
553
554
555 void map( double &x, double &y, double &z,
556 double boxX, double boxY, double boxZ ){
557
558 if(x < 0) x -= boxX * (double)( (int)( (x / boxX) - 0.5 ) );
559 else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
560
561 if(y < 0) y -= boxY * (double)( (int)( (y / boxY) - 0.5 ) );
562 else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
563
564 if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ) - 0.5 ) );
565 else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
566 }