ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp
(Generate patch)

Comparing trunk/OOPSE/utils/sysbuilder/bilayerSys.cpp (file contents):
Revision 678 by chuckv, Mon Aug 11 22:12:31 2003 UTC vs.
Revision 817 by gezelter, Fri Oct 24 17:36:18 2003 UTC

# Line 1 | Line 1
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"
# Line 12 | Line 16
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 );
33 >               double boxX, double boxY, double boxZ );
34  
35   int buildRandomBilayer( void );
36 + int buildLatticeBilayer( int isHexLattice,
37 +                         double hexSpacing,
38 +                         double aLat,
39 +                         double bLat,
40 +                         int targetNlipid,
41 +                         double targetWaterLipidRatio,
42 +                         double leafSpacing);
43  
44   void getRandomRot( double rot[3][3] );
45 + void getEulerRot( double theta, double phi, double psi, double rot[3][3] );
46 + void getUnitRot( double unit[3], double rot[3][3] );
47  
48   int buildBilayer( int isRandom ){
49    
# Line 27 | Line 51 | int buildBilayer( int isRandom ){
51      return buildRandomBilayer();
52    }
53    else{
54 <    sprintf( painCave.errMsg,
55 <             "Cannot currently create a non-random bilayer.\n" );
56 <    painCave.isFatal = 1;
57 <    simError();
34 <    return 0;
54 >
55 >    
56 >
57 >    return buildLatticeBilayer();
58    }
59   }
60  
# Line 44 | Line 67 | int buildRandomBilayer( void ){
67    } coord;
68  
69  
70 +
71    const double waterRho = 0.0334; // number density per cubic angstrom
72    const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters
73    const double waterCell = 4.929; // fcc unit cell length
74  
75 +  Lattice myFCC( FCC_LATTICE_TYPE, waterCell );
76 +  double *posX, *posY, *posZ;
77 +  double pos[3], posA[3], posB[3];
78 +
79    const double water_padding = 6.0;
80    const double lipid_spaceing = 8.0;
81  
82  
83 <  int i,j,k, l;
83 >  int i,j,k, l, m;
84    int nAtoms, atomIndex, molIndex, molID;
85    int* molSeq;
86    int* molMap;
# Line 67 | Line 95 | int buildRandomBilayer( void ){
95  
96    Atom** atoms;
97    SimInfo* simnfo;
98 +  SimState* theConfig;
99    DumpWriter* writer;
100  
101    MoleculeStamp* lipidStamp;
# Line 91 | Line 120 | int buildRandomBilayer( void ){
120    foundWater = 0;
121    for(i=0; i<bsInfo.nComponents; i++){
122      if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
123 <      
123 >          
124        foundLipid = 1;
125        lipidStamp = bsInfo.compStamps[i];
126        nLipids = bsInfo.componentsNmol[i];
127      }
128      if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
129 <      
129 >          
130        foundWater = 1;
131 <      
131 >          
132        waterStamp = bsInfo.compStamps[i];
133        nWaters = bsInfo.componentsNmol[i];
134      }
# Line 112 | Line 141 | int buildRandomBilayer( void ){
141      simError();
142    }
143    if( !foundWater ){
144 <        sprintf(painCave.errMsg,
145 <                "Could not find solvent \"%s\" in the bass file.\n",
144 >    sprintf(painCave.errMsg,
145 >            "Could not find solvent \"%s\" in the bass file.\n",
146              bsInfo.waterName );
147      painCave.isFatal = 1;
148      simError();
# Line 128 | Line 157 | int buildRandomBilayer( void ){
157    waterLocate = new MoLocator( waterStamp );
158    waterNatoms = waterStamp->getNAtoms();
159    
160 <  nAtoms = nLipids * lipidNatoms;
160 >  nAtoms = lipidNatoms;
161  
162 < simnfo[0].n_atoms = nAtoms;
163 < simnfo[0].atoms=new Atom*[nAtoms];
162 >  simnfo[0].n_atoms = nAtoms;
163 >  simnfo[0].atoms=new Atom*[nAtoms];
164  
165 < (simnfo[0]->getConfiguration())->createArrays( simnfo[0].n_atoms );
166 < for(i=0; i<simnfo[0].n_atoms; i++) simnfo[0].atoms[i]->setCoords();
165 >  theConfig = simnfo[0].getConfiguration();
166 >  theConfig->createArrays( simnfo[0].n_atoms );
167  
168 < atoms=simnfo[0].atoms;
168 >  atoms=simnfo[0].atoms;
169          
170  
171    // create the test box for initial water displacement
# Line 160 | Line 189 | int buildRandomBilayer( void ){
189    for( i=0; i < nCells; i++ ){
190      for( j=0; j < nCells; j++ ){
191        for( k=0; k < nCells; k++ ){
192 <        
193 <        waterX[ndx] = i * waterCell + x0;
194 <        waterY[ndx] = j * waterCell + y0;
195 <        waterZ[ndx] = k * waterCell + z0;
196 <        ndx++;
197 <        
198 <        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
199 <        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
171 <        waterZ[ndx] = k * waterCell + z0;
172 <        ndx++;
173 <        
174 <        waterX[ndx] = i * waterCell + x0;
175 <        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
176 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
177 <        ndx++;
178 <        
179 <        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
180 <        waterY[ndx] = j * waterCell + y0;
181 <        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
182 <        ndx++;
192 >
193 >        myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
194 >        for(l=0; l<4; l++){
195 >          waterX[ndx]=posX[l];
196 >          waterY[ndx]=posY[l];
197 >          waterZ[ndx]=posZ[l];
198 >          ndx++;
199 >        }
200        }
201      }
202    }
# Line 202 | Line 219 | int buildRandomBilayer( void ){
219    testSite.pos[1] = 0.0;
220    testSite.pos[2] = 0.0;
221  
222 <  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
222 >  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig );
223  
224    int *isActive = new int[testWaters];
225    for(i=0; i<testWaters; i++) isActive[i] = 1;
# Line 214 | Line 231 | int buildRandomBilayer( void ){
231    
232    for(i=0; ( (i<testWaters) && isActive[i] ); i++){
233      for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
234 +      
235 +      atoms[j]->getPos( pos );
236  
237 <      dx = waterX[i] - atoms[j]->getX();
238 <      dy = waterY[i] - atoms[j]->getY();
239 <      dz = waterZ[i] - atoms[j]->getZ();
237 >      dx = waterX[i] - pos[0];
238 >      dy = waterY[i] - pos[1];
239 >      dz = waterZ[i] - pos[2];
240  
241        buildMap( dx, dy, dz, testBox, testBox, testBox );
242  
243        dx2 = dx * dx;
244        dy2 = dy * dy;
245        dz2 = dz * dz;
246 <      
246 >          
247        dSqr = dx2 + dy2 + dz2;
248        if( dSqr < rCutSqr ){
249          isActive[i] = 0;
# Line 239 | Line 258 | int buildRandomBilayer( void ){
258  
259    // find the best box size for the sim
260  
261 +  int nCellsX, nCellsY, nCellsZ;
262 +
263 +  const double boxTargetX = 66.22752;
264 +  const double boxTargetY = 60.53088;
265 +
266 +  nCellsX = (int)ceil(boxTargetX / waterCell);
267 +  nCellsY = (int)ceil(boxTargetY / waterCell);
268 +  
269    int testTot;
270    int done = 0;
271 <  ndx = 0;
271 >  nCellsZ = 0;
272    while( !done ){
273 <    
274 <    ndx++;
275 <    testTot = 4 * ndx * ndx * ndx;
276 <    
273 >        
274 >    nCellsZ++;
275 >    testTot = 4 * nCellsX * nCellsY * nCellsZ;
276 >        
277      if( testTot >= targetWaters ) done = 1;
278    }
279  
253 –  nCells = ndx;
254 –    
255 –    
280    // create the new water box to the new specifications
281    
282 <  int newWaters = nCells * nCells * nCells * 4;
282 >  int newWaters = nCellsX * nCellsY * nCellsZ * 4;
283    
284    delete[] waterX;
285    delete[] waterY;
# Line 263 | Line 287 | int buildRandomBilayer( void ){
287    
288    coord* waterSites = new coord[newWaters];
289    
290 <  double box_x = waterCell * nCells;
291 <  double box_y = waterCell * nCells;
292 <  double box_z = waterCell * nCells;
293 <    
290 >  double box_x = waterCell * nCellsX;
291 >  double box_y = waterCell * nCellsY;
292 >  double box_z = waterCell * nCellsZ;
293 >        
294    // create an fcc lattice in the water box.
295    
296    ndx = 0;
297 <  for( i=0; i < nCells; i++ ){
298 <    for( j=0; j < nCells; j++ ){
299 <      for( k=0; k < nCells; k++ ){
300 <        
301 <        waterSites[ndx].pos[0] = i * waterCell;
302 <        waterSites[ndx].pos[1] = j * waterCell;
303 <        waterSites[ndx].pos[2] = k * waterCell;
304 <        ndx++;
305 <        
306 <        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
307 <        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
284 <        waterSites[ndx].pos[2] = k * waterCell;
285 <        ndx++;
286 <        
287 <        waterSites[ndx].pos[0] = i * waterCell;
288 <        waterSites[ndx].pos[1] = j * waterCell + 0.5 * waterCell;
289 <        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
290 <        ndx++;
291 <        
292 <        waterSites[ndx].pos[0] = i * waterCell + 0.5 * waterCell;
293 <        waterSites[ndx].pos[1] = j * waterCell;
294 <        waterSites[ndx].pos[2] = k * waterCell + 0.5 * waterCell;
295 <        ndx++;
297 >  for( i=0; i < nCellsX; i++ ){
298 >    for( j=0; j < nCellsY; j++ ){
299 >      for( k=0; k < nCellsZ; k++ ){
300 >
301 >        myFCC.getLatticePoints(&posX, &posY, &posZ, i, j, k);
302 >        for(l=0; l<4; l++){
303 >          waterSites[ndx].pos[0] = posX[l];
304 >          waterSites[ndx].pos[1] = posY[l];
305 >          waterSites[ndx].pos[2] = posZ[l];
306 >          ndx++;
307 >        }
308        }
309      }
310    }  
# Line 305 | Line 317 | int buildRandomBilayer( void ){
317    int reject;
318    int testDX, acceptedDX;
319  
320 +  nAtoms = nLipids * lipidNatoms;
321 +
322 +  simnfo[1].n_atoms = nAtoms;
323 +  simnfo[1].atoms=new Atom*[nAtoms];
324 +
325 +  theConfig = simnfo[1].getConfiguration();
326 +  theConfig->createArrays( simnfo[1].n_atoms );
327 +
328 +  atoms=simnfo[1].atoms;
329 +
330    rCutSqr = lipid_spaceing * lipid_spaceing;
331  
332    for(i=0; i<nLipids; i++ ){
333      done = 0;
334      while( !done ){
335 <      
335 >          
336        lipidSites[i].pos[0] = drand48() * box_x;
337        lipidSites[i].pos[1] = drand48() * box_y;
338        lipidSites[i].pos[2] = drand48() * box_z;
339 <      
339 >          
340        getRandomRot( lipidSites[i].rot );
341 <      
341 >          
342        ndx = i * lipidNatoms;
343  
344        lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
345 <                             ndx );
346 <      
345 >                             ndx, theConfig );
346 >          
347        reject = 0;
348        for( j=0; !reject && j<i; j++){
349          for(k=0; !reject && k<lipidNatoms; k++){
350            
351            acceptedDX = j*lipidNatoms + k;
352            for(l=0; !reject && l<lipidNatoms; l++){
353 <            
353 >                
354              testDX = ndx + l;
355  
356 <            dx = atoms[testDX]->getX() - atoms[acceptedDX]->getX();
357 <            dy = atoms[testDX]->getY() - atoms[acceptedDX]->getY();
336 <            dz = atoms[testDX]->getZ() - atoms[acceptedDX]->getZ();
356 >            atoms[testDX]->getPos( posA );
357 >            atoms[acceptedDX]->getPos( posB );
358              
359 +            dx = posA[0] - posB[0];
360 +            dy = posA[1] - posB[1];
361 +            dz = posA[2] - posB[2];
362 +                
363              buildMap( dx, dy, dz, box_x, box_y, box_z );
364 <            
364 >                
365              dx2 = dx * dx;
366              dy2 = dy * dy;
367              dz2 = dz * dz;
368 <            
368 >                
369              dSqr = dx2 + dy2 + dz2;
370              if( dSqr < rCutSqr ) reject = 1;
371            }
# Line 353 | Line 378 | int buildRandomBilayer( void ){
378        }
379        else{
380          done = 1;
381 <        std::cout << i << " has been accepted\n";
381 >        std::cout << (i+1) << " has been accepted\n";
382        }
383      }
384    }
385          
386 +
387 +  // zSort of the lipid positions
388 +
389 +
390 +  vector< pair<int,double> >zSortArray;
391 +  for(i=0;i<nLipids;i++)
392 +    zSortArray.push_back( make_pair(i, lipidSites[i].pos[2]) );
393 +
394 +  sort(zSortArray.begin(),zSortArray.end(),SortCond());
395 +
396 +  ofstream outFile( "./zipper.bass", ios::app);
397 +
398 +  for(i=0; i<nLipids; i++){
399 +    outFile << "zConstraint[" << i << "]{\n"
400 +            << "  molIndex = " << zSortArray[i].first <<  ";\n"
401 +            << "  zPos = ";
402 +
403 +    if(i<32) outFile << "60.0;\n";
404 +    else outFile << "100.0;\n";
405 +
406 +    outFile << "  kRatio = 0.5;\n"
407 +            << "}\n";
408 +  }
409 +  
410 +  outFile.close();
411 +
412 +
413    // cut out the waters that overlap with the lipids.
414    
415 +
416    delete[] isActive;
417    isActive = new int[newWaters];
418    for(i=0; i<newWaters; i++) isActive[i] = 1;
# Line 369 | Line 422 | int buildRandomBilayer( void ){
422    for(i=0; ( (i<newWaters) && isActive[i] ); i++){
423      for(j=0; ( (j<nAtoms) && isActive[i] ); j++){
424  
425 <      dx = waterSites[i].pos[0] - atoms[j]->getX();
373 <      dy = waterSites[i].pos[1] - atoms[j]->getY();
374 <      dz = waterSites[i].pos[2] - atoms[j]->getZ();
425 >      atoms[j]->getPos( pos );
426  
427 +      dx = waterSites[i].pos[0] - pos[0];
428 +      dy = waterSites[i].pos[1] - pos[1];
429 +      dz = waterSites[i].pos[2] - pos[2];
430 +
431        buildMap( dx, dy, dz, box_x, box_y, box_z );
432  
433        dx2 = dx * dx;
434        dy2 = dy * dy;
435        dz2 = dz * dz;
436 <      
436 >          
437        dSqr = dx2 + dy2 + dz2;
438        if( dSqr < rCutSqr ){
439          isActive[i] = 0;
440          n_active--;
441 +
442 +
443        }
444      }
445    }
446  
447 +
448 +
449 +
450    if( n_active < nWaters ){
451 <    
451 >        
452      sprintf( painCave.errMsg,
453               "Too many waters were removed, edit code and try again.\n" );
454 <    
454 >        
455      painCave.isFatal = 1;
456      simError();
457    }
# Line 404 | Line 464 | int buildRandomBilayer( void ){
464      if( isActive[quickKill] ){
465        isActive[quickKill] = 0;
466        n_active--;
467 +
468      }
469    }
470  
471    if( n_active != nWaters ){
472 <    
472 >        
473      sprintf( painCave.errMsg,
474               "QuickKill didn't work right. n_active = %d, and nWaters = %d\n",
475               n_active, nWaters );
# Line 418 | Line 479 | int buildRandomBilayer( void ){
479  
480    // clean up our messes before building the final system.
481  
482 <  for(i=0; i<nAtoms; i++){
483 <    
423 <    delete atoms[i];
424 <  }
425 <  Atom::destroyArrays();
426 <  
482 >  simnfo[0].getConfiguration()->destroyArrays();
483 >  simnfo[1].getConfiguration()->destroyArrays();
484  
485    // create the real Atom arrays
486    
# Line 443 | Line 500 | int buildRandomBilayer( void ){
500      nAtoms += waterNatoms;
501    }
502    
503 +  theConfig = simnfo[2].getConfiguration();
504 +  theConfig->createArrays( nAtoms );
505 +  simnfo[2].atoms = new Atom*[nAtoms];
506 +  atoms = simnfo[2].atoms;
507 +  simnfo[2].n_atoms = nAtoms;
508    
447 –  Atom::createArrays( nAtoms );
448 –  atoms = new Atom*[nAtoms];
449 –
450 –  
509    // initialize lipid positions
510  
511    molIndex = 0;
512    for(i=0; i<nLipids; i++ ){
513      lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
514 <                           molStart[molIndex] );
514 >                           molStart[molIndex], theConfig );
515      molIndex++;
516    }
517  
518    // initialize the water positions
519  
520    for(i=0; i<newWaters; i++){
521 <    
521 >        
522      if( isActive[i] ){
523 <      
523 >          
524        getRandomRot( waterSites[i].rot );
525        waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
526 <                             molStart[molIndex] );
526 >                             molStart[molIndex], theConfig );
527        molIndex++;
528      }
529    }  
530  
531    // set up the SimInfo object
532    
533 +  double Hmat[3][3];
534 +
535 +  Hmat[0][0] = box_x;
536 +  Hmat[0][1] = 0.0;
537 +  Hmat[0][2] = 0.0;
538 +
539 +  Hmat[1][0] = 0.0;
540 +  Hmat[1][1] = box_y;
541 +  Hmat[1][2] = 0.0;
542 +
543 +  Hmat[2][0] = 0.0;
544 +  Hmat[2][1] = 0.0;
545 +  Hmat[2][2] = box_z;
546 +  
547 +
548    bsInfo.boxX = box_x;
549    bsInfo.boxY = box_y;
550    bsInfo.boxZ = box_z;
551    
552 <  double boxVector[3];
480 <
481 <  boxVector[0] = bsInfo.boxX;
482 <  boxVector[1] = bsInfo.boxY;
483 <  boxVector[2] = bsInfo.boxZ;
484 <  simnfo->setBox( boxVector );
552 >  simnfo[2].setBoxM( Hmat );
553  
554 <  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
555 <  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
554 >  sprintf( simnfo[2].sampleName, "%s.dump", bsInfo.outPrefix );
555 >  sprintf( simnfo[2].finalName, "%s.init", bsInfo.outPrefix );
556  
489 –  simnfo->atoms = atoms;
490 –  
557    // set up the writer and write out
558    
559 <  writer = new DumpWriter( simnfo );
559 >  writer = new DumpWriter( &simnfo[2] );
560    writer->writeFinal( 0.0 );
561 <    
561 >        
562    // clean up the memory
563    
564 < //   if( molMap != NULL )   delete[] molMap;
565 < //   if( cardDeck != NULL ) delete[] cardDeck;
566 < //   if( locate != NULL ){
567 < //     for(i=0; i<bsInfo.nComponents; i++){
568 < //       delete locate[i];
569 < //     }
570 < //     delete[] locate;
571 < //   }
572 < //   if( atoms != NULL ){
573 < //     for(i=0; i<nAtoms; i++){
574 < //       delete atoms[i];
575 < //     }
576 < //     Atom::destroyArrays();
577 < //     delete[] atoms;
578 < //   }
579 < //   if( molSeq != NULL ) delete[] molSeq;
580 < //   if( simnfo != NULL ) delete simnfo;
581 < //   if( writer != NULL ) delete writer;
564 >  //     if( molMap != NULL )   delete[] molMap;
565 >  //     if( cardDeck != NULL ) delete[] cardDeck;
566 >  //     if( locate != NULL ){
567 >  //       for(i=0; i<bsInfo.nComponents; i++){
568 >  //             delete locate[i];
569 >  //       }
570 >  //       delete[] locate;
571 >  //     }
572 >  //     if( atoms != NULL ){
573 >  //       for(i=0; i<nAtoms; i++){
574 >  //             delete atoms[i];
575 >  //       }
576 >  //       Atom::destroyArrays();
577 >  //       delete[] atoms;
578 >  //     }
579 >  //     if( molSeq != NULL ) delete[] molSeq;
580 >  //     if( simnfo != NULL ) delete simnfo;
581 >  //     if( writer != NULL ) delete writer;
582  
583    return 1;
584   }
585  
586 <  
586 > int buildLatticeBilayer(int isHexLattice,
587 >                        double hexSpacing,
588 >                        double aLat,
589 >                        double bLat,
590 >                        int targetNlipid,
591 >                        double targetWaterLipidRatio,
592 >                        double leafSpacing){
593  
594 < int Old_buildRandomBilayer( void ){
594 >  typedef struct{
595 >    double rot[3][3];
596 >    double pos[3];
597 >  } coord;
598  
599 <  int i,j,k;
600 <  int nAtoms, atomIndex, molIndex, molID;
526 <  int* molSeq;
527 <  int* molMap;
528 <  int* molStart;
529 <  int* cardDeck;
530 <  int deckSize;
531 <  int rSite, rCard;
532 <  double cell;
533 <  int nCells, nSites, siteIndex;
534 <  double rot[3][3];
535 <  double pos[3];
536 <
537 <  Atom** atoms;
538 <  SimInfo* simnfo;
539 <  DumpWriter* writer;
540 <  MoLocator** locate;
541 <
542 <  // initialize functions and variables
599 >  const double waterRho = 0.0334; // number density per cubic angstrom
600 >  const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters
601    
602 <  srand48( RAND_SEED );
545 <  molSeq = NULL;
546 <  molStart = NULL;
547 <  molMap = NULL;
548 <  cardDeck = NULL;
549 <  atoms = NULL;
550 <  locate = NULL;
551 <  simnfo = NULL;
552 <  writer = NULL;
602 >  double waterCell[3];
603  
604 <  // calculate the number of cells in the fcc box
604 >  double *posX, *posY, *posZ;
605 >  double pos[3], posA[3], posB[3];
606  
607 <  nCells = 0;
557 <  nSites = 0;
558 <  while( nSites < bsInfo.totNmol ){
559 <    nCells++;
560 <    nSites = 4.0 * pow( (double)nCells, 3.0 );
561 <  }
607 >  const double waterFudge = 5.0;
608  
609 <
610 <  // create the molMap and cardDeck arrays
609 >  int i,j,k,l;
610 >  int nAtoms, atomIndex, molIndex, molID;
611 >  int* molSeq;
612 >  int* molMap;
613 >  int* molStart;
614 >  int testTot, done;
615 >  int nCells, nCellsX, nCellsY, nCellsZ;
616 >  int nx, ny;
617 >  double boxX, boxY, boxZ;
618 >  double unitVector[3];
619 >  int which;
620 >  int targetWaters;
621    
566 –  molMap = new int[nSites];
567 –  cardDeck = new int[nSites];
622    
569 –  for(i=0; i<nSites; i++){
570 –    molMap[i] = -1;
571 –    cardDeck[i] = i;
572 –  }
623  
624 <  // randomly place the molecules on the sites
575 <  
576 <  deckSize = nSites;
577 <  for(i=0; i<bsInfo.totNmol; i++){
578 <    rCard = (int)( deckSize * drand48() );
579 <    rSite = cardDeck[rCard];
580 <    molMap[rSite] = i;
624 >  coord testSite;
625  
626 <    // book keep the card deck;
627 <    
628 <    deckSize--;
629 <    cardDeck[rCard] = cardDeck[deckSize];
630 <  }
631 <  
632 <  
633 <  // create the MoLocator and Atom arrays
634 <  
635 <  nAtoms = 0;
636 <  molIndex = 0;
637 <  locate = new MoLocator*[bsInfo.nComponents];
638 <  molSeq = new int[bsInfo.totNmol];
639 <  molStart = new int[bsInfo.totNmol];  
626 >  Atom** atoms;
627 >  SimInfo* simnfo;
628 >  SimState* theConfig;
629 >  DumpWriter* writer;
630 >
631 >  MoleculeStamp* lipidStamp;
632 >  MoleculeStamp* waterStamp;  
633 >  MoLocator *lipidLocate;
634 >  MoLocator *waterLocate;
635 >  int foundLipid, foundWater;
636 >  int nLipids, lipidNatoms, nWaters, waterNatoms;
637 >  
638 >  srand48( RAND_SEED );
639 >
640 >  // create the simInfo objects
641 >
642 >  simnfo = new SimInfo;
643 >
644 >  // set the the lipidStamp
645 >
646 >  foundLipid = 0;
647 >  foundWater = 0;
648    for(i=0; i<bsInfo.nComponents; i++){
649 <    locate[i] = new MoLocator( bsInfo.compStamps[i] );
650 <    for(j=0; j<bsInfo.componentsNmol[i]; j++){
651 <      molSeq[molIndex] = i;
652 <      molStart[molIndex] = nAtoms;
653 <      molIndex++;
654 <      nAtoms += bsInfo.compStamps[i]->getNAtoms();
649 >    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
650 >          
651 >      foundLipid = 1;
652 >      lipidStamp = bsInfo.compStamps[i];
653 >      nLipids = bsInfo.componentsNmol[i];
654 >      lipidNatoms = lipidStamp->getNAtoms();
655      }
656 +    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
657 +          
658 +      foundWater = 1;
659 +          
660 +      waterStamp = bsInfo.compStamps[i];
661 +      nWaters = bsInfo.componentsNmol[i];
662 +      waterNatoms = waterStamp->getNAtoms();
663 +    }
664    }
665 +  if( !foundLipid ){
666 +    sprintf(painCave.errMsg,
667 +            "Could not find lipid \"%s\" in the bass file.\n",
668 +            bsInfo.lipidName );
669 +    painCave.isFatal = 1;
670 +    simError();
671 +  }
672 +  if( !foundWater ){
673 +    sprintf(painCave.errMsg,
674 +            "Could not find solvent \"%s\" in the bass file.\n",
675 +            bsInfo.waterName );
676 +    painCave.isFatal = 1;
677 +    simError();
678 +  }
679    
680 <  Atom::createArrays( nAtoms );
607 <  atoms = new Atom*[nAtoms];
680 >  //create the Molocator arrays
681    
682 +  lipidLocate = new MoLocator( lipidStamp );
683 +  waterLocate = new MoLocator( waterStamp );
684 +
685 +
686 +  // set up the bilayer leaves
687    
688 <  // place the molecules at each FCC site
689 <  
690 <  cell = 5.0;
691 <  for(i=0; i<bsInfo.nComponents; i++){
614 <    if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
615 <  }
616 <  cell *= 1.2; // add a little buffer
688 >  if (isHexLattice) {
689 >    aLat = sqrt(3.0)*hexSpacing;
690 >    bLat = hexSpacing;
691 >  }
692  
693 <  cell *= M_SQRT2;
693 >  nCells = (int) sqrt( (double)targetNlipid * bLat / (4.0 * aLat) );
694  
695 <  siteIndex = 0;
696 <  for(i=0; i<nCells; i++){
622 <    for(j=0; j<nCells; j++){
623 <      for(k=0; k<nCells; k++){
624 <        
625 <        if( molMap[siteIndex] >= 0 ){
626 <          pos[0] = i * cell;
627 <          pos[1] = j * cell;
628 <          pos[2] = k * cell;
629 <          
630 <          getRandomRot( rot );
631 <          molID = molSeq[molMap[siteIndex]];
632 <          atomIndex = molStart[ molMap[siteIndex] ];
633 <          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
634 <        }
635 <        siteIndex++;
695 >  nx = nCells;
696 >  ny = (int) ((double)nCells  * aLat / bLat);
697  
698 <        if( molMap[siteIndex] >= 0 ){
699 <          pos[0] = i * cell + (0.5 * cell);
639 <          pos[1] = j * cell;
640 <          pos[2] = k * cell + (0.5 * cell);
698 >  boxX = nx * aLat;
699 >  boxY = ny * bLat;  
700  
701 <          getRandomRot( rot );
702 <          molID = molSeq[molMap[siteIndex]];
644 <          atomIndex = molStart[ molMap[siteIndex] ];
645 <          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
646 <        }
647 <        siteIndex++;
701 >  nLipids = 4 * nx * ny;
702 >  coord* lipidSites = new coord[nLipids];
703  
704 <        if( molMap[siteIndex] >= 0 ){
705 <          pos[0] = i * cell + (0.5 * cell);
651 <          pos[1] = j * cell + (0.5 * cell);
652 <          pos[2] = k * cell;
653 <          
654 <          getRandomRot( rot );
655 <          molID = molSeq[molMap[siteIndex]];
656 <          atomIndex = molStart[ molMap[siteIndex] ];
657 <          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
658 <        }
659 <        siteIndex++;
704 >  unitVector[0] = 0.0;
705 >  unitVector[1] = 0.0;
706  
707 <        if( molMap[siteIndex] >= 0 ){
662 <          pos[0] = i * cell;
663 <          pos[1] = j * cell + (0.5 * cell);
664 <          pos[2] = k * cell + (0.5 * cell);
707 >  which = 0;
708  
709 <          getRandomRot( rot );
710 <          molID = molSeq[molMap[siteIndex]];
711 <          atomIndex = molStart[ molMap[siteIndex] ];
712 <          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
709 >  for (i = 0; i < nx; i++) {
710 >    for (j = 0; j < ny; j++ ) {
711 >      for (k = 0; k < 2; k++) {
712 >        
713 >        lipidSites[which].pos[0] = (double)i * aLat;
714 >        lipidSites[which].pos[1] = (double)j * bLat;
715 >        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing / 2.0);
716 >        
717 >        unitVector[2] = 2.0 * (double)k  - 1.0;
718 >        
719 >        getUnitRot( unitVector, lipidSites[which].rot );
720 >
721 >        which++;
722 >
723 >        lipidSites[which].pos[0] = aLat * ((double)i + 0.5);
724 >        lipidSites[which].pos[1] = bLat * ((double)j + 0.5);
725 >        lipidSites[which].pos[2] = ((double)k - 0.5) * (leafSpacing / 2.0);
726 >        
727 >        unitVector[2] = 2.0 * (double)k  - 1.0;
728 >        
729 >        getUnitRot( unitVector, lipidSites[which].rot );
730 >        
731 >        which++;
732 >      }
733 >    }
734 >  }
735 >
736 >  targetWaters = targetWaterLipidRatio * nLipids;
737 >
738 >  // guess the size of the water box
739 >  
740 >  
741 >
742 >  nCellsX = (int)ceil(boxX / pow(waterVol, ( 1.0 / 3.0 )) );
743 >  nCellsY = (int)ceil(boxY / pow(waterVol, ( 1.0 / 3.0 )) );
744 >
745 >  done = 0;
746 >  nCellsZ = 0;
747 >  while( !done ){
748 >        
749 >    nCellsZ++;
750 >    testTot = 4 * nCellsX * nCellsY * nCellsZ;
751 >        
752 >    if( testTot >= targetWaters ) done = 1;
753 >  }
754 >  
755 >  nWaters = nCellsX * nCellsY * nCellsZ * 4;
756 >  
757 >  coord* waterSites = new coord[nWaters];
758 >
759 >  waterCell[0] = boxX / nCellsX;
760 >  waterCell[1] = boxY / nCellsY;
761 >  waterCell[2] = 4.0 / (waterRho * waterCell[0] * waterCell[1]);
762 >
763 >  Lattice *myORTHO;
764 >  myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell);
765 >  myORTHO->setStartZ( leafSpacing / 2.0 + waterFudge);
766 >
767 >  boxZ = waterCell[2] * nCellsZ;
768 >        
769 >  // create an fcc lattice in the water box.
770 >  
771 >  which = 0;
772 >  for( i=0; i < nCellsX; i++ ){
773 >    for( j=0; j < nCellsY; j++ ){
774 >      for( k=0; k < nCellsZ; k++ ){
775 >
776 >        myORTHO->getLatticePoints(&posX, &posY, &posZ, i, j, k);
777 >        for(l=0; l<4; l++){
778 >          waterSites[which].pos[0] = posX[l];
779 >          waterSites[which].pos[1] = posY[l];
780 >          waterSites[which].pos[2] = posZ[l];
781 >          which++;
782          }
671 –        siteIndex++;
783        }
784      }
785 +  }  
786 +
787 +  // create the real Atom arrays
788 +  
789 +  nAtoms = 0;
790 +  molIndex = 0;
791 +  molStart = new int[nLipids + nWaters];  
792 +  
793 +  for(j=0; j<nLipids; j++){
794 +    molStart[molIndex] = nAtoms;
795 +    molIndex++;
796 +    nAtoms += lipidNatoms;
797    }
798  
799 +  for(j=0; j<nWaters; j++){
800 +    molStart[molIndex] = nAtoms;
801 +    molIndex++;
802 +    nAtoms += waterNatoms;
803 +  }
804 +  
805 +  theConfig = simnfo->getConfiguration();
806 +  theConfig->createArrays( nAtoms );
807 +  simnfo->atoms = new Atom*[nAtoms];
808 +  atoms = simnfo->atoms;
809 +
810 +  // initialize lipid positions
811 +
812 +  molIndex = 0;
813 +  for(i=0; i<nLipids; i++ ){
814 +    lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms,
815 +                           molStart[molIndex], theConfig );
816 +    molIndex++;
817 +  }
818 +
819 +  // initialize the water positions
820 +
821 +  for(i=0; i<nWaters; i++){
822 +    
823 +    getRandomRot( waterSites[i].rot );
824 +    waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms,
825 +                           molStart[molIndex], theConfig );
826 +    molIndex++;
827 +  }
828 +
829    // set up the SimInfo object
830    
831 <  bsInfo.boxX = nCells * cell;
832 <  bsInfo.boxY = nCells * cell;
833 <  bsInfo.boxZ = nCells * cell;
831 >  double Hmat[3][3];
832 >
833 >  Hmat[0][0] = boxX;
834 >  Hmat[0][1] = 0.0;
835 >  Hmat[0][2] = 0.0;
836 >
837 >  Hmat[1][0] = 0.0;
838 >  Hmat[1][1] = boxY;
839 >  Hmat[1][2] = 0.0;
840 >
841 >  Hmat[2][0] = 0.0;
842 >  Hmat[2][1] = 0.0;
843 >  Hmat[2][2] = boxZ;
844    
682 –  double boxVector[3];
683 –  simnfo = new SimInfo();
684 –  simnfo->n_atoms = nAtoms;
685 –  boxVector[0] = bsInfo.boxX;
686 –  boxVector[1] = bsInfo.boxY;
687 –  boxVector[2] = bsInfo.boxZ;
688 –  simnfo->setBox( boxVector );
845  
846 +  bsInfo.boxX = boxX;
847 +  bsInfo.boxY = boxY;
848 +  bsInfo.boxZ = boxZ;
849 +  
850 +  simnfo->setBoxM( Hmat );
851 +
852    sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
853    sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
854  
693 –  simnfo->atoms = atoms;
694 –  
855    // set up the writer and write out
856    
857    writer = new DumpWriter( simnfo );
858 <  writer->writeFinal(0.0);
699 <    
700 <  // clean up the memory
701 <  
702 <  if( molMap != NULL )   delete[] molMap;
703 <  if( cardDeck != NULL ) delete[] cardDeck;
704 <  if( locate != NULL ){
705 <    for(i=0; i<bsInfo.nComponents; i++){
706 <      delete locate[i];
707 <    }
708 <    delete[] locate;
709 <  }
710 <  if( atoms != NULL ){
711 <    for(i=0; i<nAtoms; i++){
712 <      delete atoms[i];
713 <    }
714 <    Atom::destroyArrays();
715 <    delete[] atoms;
716 <  }
717 <  if( molSeq != NULL ) delete[] molSeq;
718 <  if( simnfo != NULL ) delete simnfo;
719 <  if( writer != NULL ) delete writer;
858 >  writer->writeFinal( 0.0 );
859  
860    return 1;
861   }
# Line 735 | Line 874 | void getRandomRot( double rot[3][3] ){
874  
875    theta = acos( cosTheta );
876  
877 +  getEulerRot( theta, phi, psi, rot );
878 + }
879 +
880 +
881 + void getEulerRot( double theta, double phi, double psi, double rot[3][3] ){
882 +
883    rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
884    rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
885    rot[0][2] = sin(theta) * sin(psi);
# Line 747 | Line 892 | void getRandomRot( double rot[3][3] ){
892    rot[2][1] = -cos(phi) * sin(theta);
893    rot[2][2] = cos(theta);  
894   }
895 +
896 +
897 + void getUnitRot( double u[3], double rot[3][3] ){
898 +
899 +  double theta, phi, psi;
900 +
901 +  theta = acos(u[2]);
902 +  phi = atan(u[1] / u[0]);
903 +  psi = 0.0;
904 +
905 +  getEulerRot( theta, phi, psi, rot );
906 + }
907                          
908  
909  
910   void buildMap( double &x, double &y, double &z,
911 <          double boxX, double boxY, double boxZ ){
911 >               double boxX, double boxY, double boxZ ){
912    
913    if(x < 0) x -= boxX * (double)( (int)( (x / boxX)  - 0.5 ) );
914    else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines