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 32 by mmeineke, Wed Jul 17 20:33:56 2002 UTC

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

Diff Legend

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