ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp
Revision: 821
Committed: Fri Oct 24 22:17:45 2003 UTC (21 years, 6 months ago) by mmeineke
File size: 12623 byte(s)
Log Message:
put QuickBass, MoLocator, and latticeBuilder into a Builder Library
overhauled latticeBilayer into its own program. Removed sysBuild from the Makefile

File Contents

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