ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp
Revision: 700
Committed: Mon Aug 18 20:59:47 2003 UTC (21 years, 8 months ago) by chuckv
File size: 12538 byte(s)
Log Message:
Fixed sysBuild -bilayer works. Nanobuilder still broke.

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