ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/makeLipid/makeRandomLipid.cpp
(Generate patch)

Comparing trunk/makeLipid/makeRandomLipid.cpp (file contents):
Revision 29 by mmeineke, Tue Jul 16 01:22:04 2002 UTC vs.
Revision 33 by mmeineke, Wed Jul 17 20:36:25 2002 UTC

# Line 12 | Line 12
12   #include "ReadWrite.hpp"
13  
14  
15 < void map( double *x, double *y, double *z, double centerX, double centerY,
16 <          double centerZ, double boxX, double boxY, double boxZ );
15 > void map( double &x, double &y, double &z,
16 >          double boxX, double boxY, double boxZ );
17  
18 + void rotate( double &x, double &y, double &z,
19 +             double theta, double phi, double psi );
20 +
21   char* program_name;
22   using namespace std;
23  
# Line 45 | Line 48 | int main(int argc,char* argv[]){
48  
49    double water_shell = 10.0;
50    double water_padding = 2.5;
51 <  double lipid_spaceing = 4.0;
51 >  double lipid_spaceing = 2.5;
52    
53    srand48( 1337 ); // initialize the random number generator.
54  
# Line 152 | Line 155 | int main(int argc,char* argv[]){
155    // create an fcc lattice in the water box.
156  
157    
158 <  index = 0;
158 >  int ndx = 0;
159    for( i=0; i < n_cellX; i++ ){
160      for( j=0; j < n_cellY; j++ ){
161        for( k=0; k < n_cellZ; k++ ){
162          
163 <        waterX[index] = i * water_cell + x0;
164 <        waterY[index] = j * water_cell + y0;
165 <        waterZ[index] = k * water_cell + z0;
166 <        index++;
163 >        waterX[ndx] = i * water_cell + x0;
164 >        waterY[ndx] = j * water_cell + y0;
165 >        waterZ[ndx] = k * water_cell + z0;
166 >        ndx++;
167          
168 <        waterX[index] = i * water_cell + 0.5 * water_cell + x0;
169 <        waterY[index] = j * water_cell + 0.5 * water_cell + y0;
170 <        waterZ[index] = k * water_cell + z0;
171 <        index++;
168 >        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
169 >        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
170 >        waterZ[ndx] =   k * water_cell + z0;
171 >        ndx++;
172          
173 <        waterX[index] = i * water_cell + x0;
174 <        waterY[index] = j * water_cell + 0.5 * water_cell + y0;
175 <        waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
176 <        index++;
173 >        waterX[ndx] = i * water_cell + x0;
174 >        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
175 >        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
176 >        ndx++;
177          
178 <        waterX[index] = i * water_cell + 0.5 * water_cell + x0;
179 <        waterY[index] = j * water_cell + y0;
180 <        waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
181 <        index++;
178 >        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
179 >        waterY[ndx] = j * water_cell + y0;
180 >        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
181 >        ndx++;
182        }
183      }
184    }
# Line 198 | Line 201 | int main(int argc,char* argv[]){
201        dy = waterY[i] - lipidAtoms[j]->getY();
202        dz = waterZ[i] - lipidAtoms[j]->getZ();
203  
204 <      map( &dx, &dy, &dz, cx, cy, cz, box_x, box_y, box_z );
204 >      map( dx, dy, dz, box_x, box_y, box_z );
205  
206        dx2 = dx * dx;
207        dy2 = dy * dy;
# Line 211 | Line 214 | int main(int argc,char* argv[]){
214        }
215      }
216    }
217 <
217 >  
218    n_h2o_target += n_deleted * n_lipids;
219  
220    // find a box size that best suits the number of waters we need.
221  
222    int done = 0;
223    
224 <  if( n_waters < n_h2o_target ){
224 >  if( n_water < n_h2o_target ){
225      
226      int n_generated = n_cellX;
227      int n_test, nx, ny, nz;
# Line 226 | Line 229 | int main(int argc,char* argv[]){
229      ny = n_cellY;
230      nz = n_cellZ;
231      
232 +    n_test = 4 * nx * ny * nz;
233 +
234      while( n_test < n_h2o_target ){
235        
236        nz++;
# Line 234 | Line 239 | int main(int argc,char* argv[]){
239      
240      int n_diff, goodX, goodY, goodZ;
241      
242 <    n_diff = ntest - n_h2o_target;
242 >    n_diff = n_test - n_h2o_target;
243      goodX = nx;
244      goodY = ny;
245      goodZ = nz;
246      
247      int test_diff;
248 <    int n_limit = n_z;
249 <    n_z = n_cellZ;
248 >    int n_limit = nz;
249 >    nz = n_cellZ;
250      
251      for( i=n_generated; i<=n_limit; i++ ){
252        for( j=i; j<=n_limit; j++ ){
# Line 270 | Line 275 | int main(int argc,char* argv[]){
275      n_cellZ = goodZ;
276    }
277  
278 <  // we now have the best box size for the simulation.
279 <    
280 <    
278 >  // we now have the best box size for the simulation. Next we
279 >  // recreate the water box to the new specifications.
280 >
281 >  n_water = n_cellX * n_cellY * n_cellZ * 4;
282 >  
283 >  delete[] waterX;
284 >  delete[] waterY;
285 >  delete[] waterZ;
286 >  
287 >  waterX = new double[n_water];
288 >  waterY = new double[n_water];
289 >  waterZ = new double[n_water];
290 >
291 >  box_x = water_cell * n_cellX;
292 >  box_y = water_cell * n_cellY;
293 >  box_z = water_cell * n_cellZ;
294 >    
295 >  x0 = 0.0;
296 >  y0 = 0.0;
297 >  z0 = 0.0;
298 >
299 >  cx = ( box_x * 0.5 );
300 >  cy = ( box_y * 0.5 );
301 >  cz = ( box_z * 0.5 );
302 >
303 >  // create an fcc lattice in the water box.
304 >  
305 >  ndx = 0;
306 >  for( i=0; i < n_cellX; i++ ){
307 >    for( j=0; j < n_cellY; j++ ){
308 >      for( k=0; k < n_cellZ; k++ ){
309 >        
310 >        waterX[ndx] = i * water_cell + x0;
311 >        waterY[ndx] = j * water_cell + y0;
312 >        waterZ[ndx] = k * water_cell + z0;
313 >        ndx++;
314 >        
315 >        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
316 >        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
317 >        waterZ[ndx] =   k * water_cell + z0;
318 >        ndx++;
319 >        
320 >        waterX[ndx] = i * water_cell + x0;
321 >        waterY[ndx] =   j * water_cell + 0.5 * water_cell + y0;
322 >        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
323 >        ndx++;
324 >        
325 >        waterX[ndx] = i * water_cell + 0.5 * water_cell + x0;
326 >        waterY[ndx] = j * water_cell + y0;
327 >        waterZ[ndx] = k * water_cell + 0.5 * water_cell + z0;
328 >        ndx++;
329 >      }
330 >    }
331 >  }
332 >  
333 >  // **************************************************************  
334 >  
335 >
336 >  
337 >  // start a 3D RSA for the for the lipid placements
338 >  
339 >  srand48( 1337 );
340 >
341 >  int rsaNAtoms = n_lipids * lipidNAtoms;
342 >  Atom** rsaAtoms = new Atom*[rsaNAtoms];
343            
344 +  DirectionalAtom* dAtom;
345 +  DirectionalAtom* dAtomNew;
346  
347 +  double rotMat[3][3];
348 +  double unitRotMat[3][3];
349 +  
350 +  unitRotMat[0][0] = 1.0;
351 +  unitRotMat[0][1] = 0.0;
352 +  unitRotMat[0][2] = 0.0;
353 +  
354 +  unitRotMat[1][0] = 0.0;
355 +  unitRotMat[1][1] = 1.0;
356 +  unitRotMat[1][2] = 0.0;
357 +  
358 +  unitRotMat[2][0] = 0.0;
359 +  unitRotMat[2][1] = 0.0;
360 +  unitRotMat[2][2] = 1.0;
361  
362 +  ndx = 0;
363 +  for(i=0; i<n_lipids; i++ ){
364 +    for(j=0; j<lipidNAtoms; j++){
365 +      
366 +      if( lipidAtoms[j]->isDirectional() ){
367 +        dAtom = (DirectionalAtom *)lipidAtoms[j];
368 +        
369 +        dAtomNew = new DirectionalAtom();
370 +        dAtomNew->setSUx( dAtom->getSUx() );
371 +        dAtomNew->setSUx( dAtom->getSUx() );
372 +        dAtomNew->setSUx( dAtom->getSUx() );
373 +        
374 +        dAtom->getA( rotMat );
375 +        dAtomNew->setA( rotMat );
376 +        
377 +        rsaAtoms[ndx] = dAtomNew;
378 +      }
379 +      else{
380 +        
381 +        rsaAtoms[ndx] = new GeneralAtom();
382 +      }
383 +      
384 +      rsaAtoms[ndx]->setType( lipidAtoms[j]->getType() );
385 +      
386 +      ndx++;
387 +    }
388 +  }
389 +  
390 +  double testX, testY, testZ;
391 +  double theta, phi, psi;
392 +  double tempX, tempY, tempZ;
393 +  int reject;
394 +  int testDX, acceptedDX;
395 +  
396 +  rCutSqr = lipid_spaceing * lipid_spaceing;
397 +  
398 +  for(i=0; i<n_lipids; i++ ){
399 +    done = 0;
400 +    while( !done ){
401 +      
402 +      testX = drand48() * box_x;
403 +      testY = drand48() * box_y;
404 +      testZ = drand48() * box_z;
405 +      
406 +      theta = drand48() * 2.0 * M_PI;
407 +      phi   = drand48() * 2.0 * M_PI;
408 +      psi   = drand48() * 2.0 * M_PI;
409 +      
410 +      ndx = i * lipidNAtoms;
411 +      for(j=0; j<lipidNAtoms; j++){
412 +      
413 +        tempX = lipidAtoms[j]->getX();
414 +        tempY = lipidAtoms[j]->getY();
415 +        tempZ = lipidAtoms[j]->getZ();
416  
417 +        rotate( tempX, tempY, tempZ, theta, phi, psi );
418 +        
419 +        rsaAtoms[ndx + j]->setX( tempX + testX );
420 +        rsaAtoms[ndx + j]->setY( tempY + testY );
421 +        rsaAtoms[ndx + j]->setZ( tempZ + testZ );
422 +      }
423 +      
424 +      reject = 0;
425 +      for( j=0; !reject && j<i; j++){
426 +        for(k=0; !reject && k<lipidNAtoms; k++){
427 +          
428 +          acceptedDX = j*lipidNAtoms + k;
429 +          for(l=0; !reject && l<lipidNAtoms; l++){
430 +            
431 +            testDX = ndx + l;
432  
433 +            dx = rsaAtoms[testDX]->getX() - rsaAtoms[acceptedDX]->getX();
434 +            dy = rsaAtoms[testDX]->getY() - rsaAtoms[acceptedDX]->getY();
435 +            dz = rsaAtoms[testDX]->getZ() - rsaAtoms[acceptedDX]->getZ();
436 +            
437 +            map( dx, dy, dz, box_x, box_y, box_z );
438 +            
439 +            dx2 = dx * dx;
440 +            dy2 = dy * dy;
441 +            dz2 = dz * dz;
442 +            
443 +            dSqr = dx2 + dy2 + dz2;
444 +            if( dSqr < rCutSqr ) reject = 1;
445 +          }
446 +        }
447 +      }
448  
449 +      if( !reject ){
450 +        done = 1;
451 +        std::cerr << i << " has been accepted\n";
452 +      }
453 +    }
454 +  }
455 +        
456 +  // cut out the waters that overlap with the lipids.
457 +  
458 +  delete[] isActive;
459 +  isActive = new int[n_water];
460 +  for(i=0; i<n_water; i++) isActive[i] = 1;
461 +  int n_active = n_water;
462 +  rCutSqr = water_padding * water_padding;
463 +  
464 +  for(i=0; ( (i<n_water) && isActive[i] ); i++){
465 +    for(j=0; ( (j<rsaNAtoms) && isActive[i] ); j++){
466  
467 +      dx = waterX[i] - rsaAtoms[j]->getX();
468 +      dy = waterY[i] - rsaAtoms[j]->getY();
469 +      dz = waterZ[i] - rsaAtoms[j]->getZ();
470  
471 <  int new_nAtoms = group_nAtoms + n_active;
471 >      map( dx, dy, dz, box_x, box_y, box_z );
472 >
473 >      dx2 = dx * dx;
474 >      dy2 = dy * dy;
475 >      dz2 = dz * dz;
476 >      
477 >      dSqr = dx2 + dy2 + dz2;
478 >      if( dSqr < rCutSqr ){
479 >        isActive[i] = 0;
480 >        n_active--;
481 >      }
482 >    }
483 >  }
484 >
485 >  std::cerr << "final n_waters = " << n_active << "\n";
486 >
487 >  // place all of the waters and lipids into one new array
488 >
489 >  int new_nAtoms = rsaNAtoms + n_active;
490    Atom** new_atoms = new Atom*[new_nAtoms];
491    
492 <  index = 0;
493 <  for(i=0; i<group_nAtoms; i++ ){
492 >  ndx = 0;
493 >  for(i=0; i<rsaNAtoms; i++ ){
494      
495 <    if( group_atoms[i]->isDirectional() ){
496 <      dAtom = (DirectionalAtom *)group_atoms[i];
495 >    if( rsaAtoms[i]->isDirectional() ){
496 >      dAtom = (DirectionalAtom *)rsaAtoms[i];
497        
498        dAtomNew = new DirectionalAtom();
499        dAtomNew->setSUx( dAtom->getSUx() );
500        dAtomNew->setSUx( dAtom->getSUx() );
501        dAtomNew->setSUx( dAtom->getSUx() );
502        
503 +      dAtom->getA( rotMat );
504        dAtomNew->setA( rotMat );
505        
506 <      new_atoms[index] = dAtomNew;
506 >      new_atoms[ndx] = dAtomNew;
507      }
508      else{
509        
510 <      new_atoms[index] = new GeneralAtom();
510 >      new_atoms[ndx] = new GeneralAtom();
511      }
512  
513 <    new_atoms[index]->setType( group_atoms[i]->getType() );    
513 >    new_atoms[ndx]->setType( rsaAtoms[i]->getType() );    
514      
515 <    new_atoms[index]->setX( group_atoms[i]->getX() );    
516 <    new_atoms[index]->setY( group_atoms[i]->getY() );
517 <    new_atoms[index]->setZ( group_atoms[i]->getZ() );
515 >    new_atoms[ndx]->setX( rsaAtoms[i]->getX() );    
516 >    new_atoms[ndx]->setY( rsaAtoms[i]->getY() );
517 >    new_atoms[ndx]->setZ( rsaAtoms[i]->getZ() );
518  
519 <    new_atoms[index]->set_vx( 0.0 );
520 <    new_atoms[index]->set_vy( 0.0 );
521 <    new_atoms[index]->set_vz( 0.0 );
519 >    new_atoms[ndx]->set_vx( 0.0 );
520 >    new_atoms[ndx]->set_vy( 0.0 );
521 >    new_atoms[ndx]->set_vz( 0.0 );
522      
523 <    index++;
523 >    ndx++;
524    }
525  
320 –
321 –
322 –  
526    for(i=0; i<n_water; i++){
527      if(isActive[i]){
528        
529 <      new_atoms[index] = new DirectionalAtom();
530 <      new_atoms[index]->setType( "SSD" );
529 >      new_atoms[ndx] = new DirectionalAtom();
530 >      new_atoms[ndx]->setType( "SSD" );
531  
532 <      new_atoms[index]->setX( waterX[i] );    
533 <      new_atoms[index]->setY( waterY[i] );
534 <      new_atoms[index]->setZ( waterZ[i] );
532 >      new_atoms[ndx]->setX( waterX[i] );    
533 >      new_atoms[ndx]->setY( waterY[i] );
534 >      new_atoms[ndx]->setZ( waterZ[i] );
535        
536 <      new_atoms[index]->set_vx( 0.0 );
537 <      new_atoms[index]->set_vy( 0.0 );
538 <      new_atoms[index]->set_vz( 0.0 );
536 >      new_atoms[ndx]->set_vx( 0.0 );
537 >      new_atoms[ndx]->set_vy( 0.0 );
538 >      new_atoms[ndx]->set_vz( 0.0 );
539        
540 <      dAtom = (DirectionalAtom *) new_atoms[index];
540 >      dAtom = (DirectionalAtom *) new_atoms[ndx];
541  
542        dAtom->setSUx( 0.0 );
543        dAtom->setSUy( 0.0 );
544        dAtom->setSUz( 1.0 );
545        
546 <      dAtom->setA( rotMat );
546 >      dAtom->setA( unitRotMat );
547        
548 <      index++;
548 >      ndx++;
549      }
550    }
551  
# Line 390 | Line 593 | int main(int argc,char* argv[]){
593   }
594  
595  
596 < void map( x, y, z, centerX, centerY, centerZ, boxX, boxY, boxZ )
597 <     double *x, *y, *z;
598 <     double centerX, centerY, centerZ;
599 <     double boxX, boxY, boxZ;
600 < {
596 > void map( double &x, double &y, double &z,
597 >          double boxX, double boxY, double boxZ ){
598 >  
599 >  if(x < 0) x -= boxX * (double)( (int)( (x / boxX)  - 0.5 ) );
600 >  else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
601  
602 <  *x -= centerX;
603 <  *y -= centerY;
604 <  *z -= centerZ;
602 >  if(y < 0) y -= boxY * (double)( (int)( (y / boxY)  - 0.5 ) );
603 >  else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
604 >  
605 >  if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ)  - 0.5 ) );
606 >  else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
607 > }
608  
403 –  if(*x < 0) *x -= boxX * (double)( (int)( (*x / boxX)  - 0.5 ) );
404 –  else *x -= boxX * (double)( (int)( (*x / boxX ) + 0.5));
609  
610 <  if(*y < 0) *y -= boxY * (double)( (int)( (*y / boxY)  - 0.5 ) );
611 <  else *y -= boxY * (double)( (int)( (*y / boxY ) + 0.5));
610 > void rotate( double &x, double &y, double &z,
611 >             double theta, double phi, double psi ){
612 >
613 >  double newX, newY, newZ;
614 >
615 >  double A[3][3];
616 >
617 >  A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
618 >  A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
619 >  A[0][2] = sin(theta) * sin(psi);
620    
621 <  if(*z < 0) *z -= boxZ * (double)( (int)( (*z / boxZ)  - 0.5 ) );
622 <  else *z -= boxZ * (double)( (int)( (*z / boxZ ) + 0.5));
621 >  A[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
622 >  A[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
623 >  A[1][2] = sin(theta) * cos(psi);
624 >
625 >  A[2][0] = sin(phi) * sin(theta);
626 >  A[2][1] = -cos(phi) * sin(theta);
627 >  A[2][2] = cos(theta);
628 >
629 >  newX = (x * A[0][0]) + (y * A[0][1]) + (z * A[0][2]);
630 >  newY = (x * A[1][0]) + (y * A[1][1]) + (z * A[1][2]);
631 >  newZ = (x * A[2][0]) + (y * A[2][1]) + (z * A[2][2]);
632 >
633 >  x = newX;
634 >  y = newY;
635 >  z = newZ;
636   }

Diff Legend

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