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 28 by mmeineke, Mon Jul 15 20:28:33 2002 UTC vs.
Revision 29 by mmeineke, Tue Jul 16 01:22:04 2002 UTC

# Line 2 | Line 2
2   #include <cstdlib>
3   #include <cstring>
4   #include <cstdio>
5 + #include <cmath>
6  
7   #include "SimSetup.hpp"
8   #include "SimInfo.hpp"
# Line 11 | 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 );
17 +
18   char* program_name;
19   using namespace std;
20  
# Line 33 | Line 37 | int main(int argc,char* argv[]){
37    const double water_vol = 4.0 / water_rho; // volume occupied by 4 waters
38    const double water_cell = 4.929; // fcc unit cell length
39  
40 <  int n_lipidsX = 5;
41 <  int n_lipidsY = 10;
42 <  int n_lipids = n_lipidsX * n_lipidsY;
40 >  int n_lipids = 50;
41 >  double water_ratio = 25.0; // water to lipid ratio
42 >  int n_h2o_target = (int)( n_lipids * water_ratio + 0.5 );
43    
44    std::cerr << "n_lipids = " << n_lipids << "\n";
45  
# Line 67 | Line 71 | int main(int argc,char* argv[]){
71    lipidAtoms = entry_plug->atoms;
72    lipidNAtoms = entry_plug->n_atoms;
73  
70 –  int group_nAtoms = n_lipids * lipidNAtoms;
71 –  Atom** group_atoms = new Atom*[group_nAtoms];
72 –  DirectionalAtom* dAtom;
73 –  DirectionalAtom* dAtomNew;
74 –  
75 –  double rotMat[3][3];
76 –  
77 –  rotMat[0][0] = 1.0;
78 –  rotMat[0][1] = 0.0;
79 –  rotMat[0][2] = 0.0;
80 –  
81 –  rotMat[1][0] = 0.0;
82 –  rotMat[1][1] = 1.0;
83 –  rotMat[1][2] = 0.0;
84 –  
85 –  rotMat[2][0] = 0.0;
86 –  rotMat[2][1] = 0.0;
87 –  rotMat[2][2] = 1.0;
74  
75 <  int index =0;
90 <  for(i=0; i<n_lipids; i++ ){
91 <    for(j=0; j<lipidNAtoms; j++){
92 <      
93 <      if( lipidAtoms[j]->isDirectional() ){
94 <        dAtom = (DirectionalAtom *)lipidAtoms[j];
95 <        
96 <        dAtomNew = new DirectionalAtom();
97 <        dAtomNew->setSUx( dAtom->getSUx() );
98 <        dAtomNew->setSUx( dAtom->getSUx() );
99 <        dAtomNew->setSUx( dAtom->getSUx() );
100 <      
101 <        dAtomNew->setA( rotMat );
102 <        
103 <        group_atoms[index] = dAtomNew;
104 <      }
105 <      else{
106 <        
107 <        group_atoms[index] = new GeneralAtom();
108 <      }
109 <      
110 <      group_atoms[index]->setType( lipidAtoms[j]->getType() );
111 <      
112 <      index++;
113 <    }
114 <  }
75 >  // find the width, height, and length of the molecule
76  
116 –  index = 0;
117 –  for(i=0; i<n_lipidsX; i++){
118 –    for(j=0; j<n_lipidsY; j++){
119 –      for(l=0; l<lipidNAtoms; l++){
120 –        
121 –        group_atoms[index]->setX( lipidAtoms[l]->getX() +
122 –                                  i*lipid_spaceing );
123 –        
124 –        group_atoms[index]->setY( lipidAtoms[l]->getY() +
125 –                                  j*lipid_spaceing );
126 –        
127 –        group_atoms[index]->setZ( lipidAtoms[l]->getZ() );
128 –        
129 –        index++;
130 –      }
131 –    }
132 –  }
133 –
77    double min_x, min_y, min_z;
78    double max_x, max_y, max_z;
79    double test_x, test_y, test_z;
80  
81 <  max_x = min_x = group_atoms[0]->getX();
82 <  max_y = min_y = group_atoms[0]->getY();
83 <  max_z = min_z = group_atoms[0]->getZ();
81 >  max_x = min_x = lipidAtoms[0]->getX();
82 >  max_y = min_y = lipidAtoms[0]->getY();
83 >  max_z = min_z = lipidAtoms[0]->getZ();
84    
85 <  for(i=0; i<group_nAtoms; i++){
85 >  for(i=0; i<lipidNAtoms; i++){
86  
87 <    test_x = group_atoms[i]->getX();
88 <    test_y = group_atoms[i]->getY();
89 <    test_z = group_atoms[i]->getZ();
87 >    test_x = lipidAtoms[i]->getX();
88 >    test_y = lipidAtoms[i]->getY();
89 >    test_z = lipidAtoms[i]->getZ();
90  
91      if( test_x < min_x ) min_x = test_x;
92      if( test_y < min_y ) min_y = test_y;
# Line 154 | Line 97 | int main(int argc,char* argv[]){
97      if( test_z > max_z ) max_z = test_z;
98    }
99  
100 <  double box_x = max_x - min_x + 2*water_shell;
101 <  double box_y = max_y - min_y + 2*water_shell;
102 <  double box_z = max_z - min_z + 2*water_shell;
100 >  double ml2 = pow((max_x - min_x), 2 ) + pow((max_y - min_y), 2 )
101 >    + pow((max_x - min_x), 2 );
102 >  double max_length = sqrt( ml2 );
103  
161 –  int n_cellX = (int)(box_x / water_cell + 0.5 );
162 –  int n_cellY = (int)(box_y / water_cell + 0.5 );
163 –  int n_cellZ = (int)(box_z / water_cell + 0.5 );
104    
105 +  // from this information, create the test box
106 +
107 +
108 +  double box_x;
109 +  double box_y;
110 +  double box_z;
111 +
112 +  box_x = box_y = box_z = max_length + water_cell * 4.0; // pad with 4 cells
113 +
114 +  int n_cellX = (int)(box_x / water_cell + 1.0 );
115 +  int n_cellY = (int)(box_y / water_cell + 1.0 );
116 +  int n_cellZ = (int)(box_z / water_cell + 1.0 );
117 +  
118    box_x = water_cell * n_cellX;
119    box_y = water_cell * n_cellY;
120    box_z = water_cell * n_cellZ;
# Line 172 | Line 125 | int main(int argc,char* argv[]){
125    double *waterY = new double[n_water];
126    double *waterZ = new double[n_water];
127    
128 +  
129 +  // find the center of the test lipid, and make it the center of our
130 +  // soon to be created water box.
131 +
132 +
133    double cx, cy, cz;
134  
135    cx = 0.0;
136    cy = 0.0;
137    cz = 0.0;
138 <  for(i=0; i<group_nAtoms; i++){
139 <    cx += group_atoms[i]->getX();
140 <    cy += group_atoms[i]->getY();
141 <    cz += group_atoms[i]->getZ();
138 >  for(i=0; i<lipidNAtoms; i++){
139 >    cx += lipidAtoms[i]->getX();
140 >    cy += lipidAtoms[i]->getY();
141 >    cz += lipidAtoms[i]->getZ();
142    }
143 <  cx /= group_nAtoms;
144 <  cy /= group_nAtoms;
145 <  cz /= group_nAtoms;
143 >  cx /= lipidNAtoms;
144 >  cy /= lipidNAtoms;
145 >  cz /= lipidNAtoms;
146    
147    double x0 = cx - ( box_x * 0.5 );
148    double y0 = cy - ( box_y * 0.5 );
149    double z0 = cz - ( box_z * 0.5 );
150 +
151 +
152 +  // create an fcc lattice in the water box.
153 +
154    
155    index = 0;
156    for( i=0; i < n_cellX; i++ ){
# Line 218 | Line 180 | int main(int argc,char* argv[]){
180      }
181    }
182  
183 +
184 +  // calculate the number of water's displaced by our molecule.
185 +
186    int *isActive = new int[n_water];
187    for(i=0; i<n_water; i++) isActive[i] = 1;
188    
189 <  int n_active = n_water;
189 >  int n_deleted = 0;
190    double dx, dy, dz;
191    double dx2, dy2, dz2, dSqr;
192    double rCutSqr = water_padding * water_padding;
193    
194    for(i=0; ( (i<n_water) && isActive[i] ); i++){
195 <    for(j=0; ( (j<group_nAtoms) && isActive[i] ); j++){
195 >    for(j=0; ( (j<lipidNAtoms) && isActive[i] ); j++){
196  
197 <      dx = waterX[i] - group_atoms[j]->getX();
198 <      dy = waterY[i] - group_atoms[j]->getY();
199 <      dz = waterZ[i] - group_atoms[j]->getZ();
197 >      dx = waterX[i] - lipidAtoms[j]->getX();
198 >      dy = waterY[i] - lipidAtoms[j]->getY();
199 >      dz = waterZ[i] - lipidAtoms[j]->getZ();
200  
201 +      map( &dx, &dy, &dz, cx, cy, cz, box_x, box_y, box_z );
202 +
203        dx2 = dx * dx;
204        dy2 = dy * dy;
205        dz2 = dz * dz;
# Line 240 | Line 207 | int main(int argc,char* argv[]){
207        dSqr = dx2 + dy2 + dz2;
208        if( dSqr < rCutSqr ){
209          isActive[i] = 0;
210 <        n_active--;
210 >        n_deleted++;
211        }
212      }
213    }
214  
215 <  std::cerr << "final n_waters = " << n_active << "\n";
215 >  n_h2o_target += n_deleted * n_lipids;
216  
217 +  // find a box size that best suits the number of waters we need.
218 +
219 +  int done = 0;
220 +  
221 +  if( n_waters < n_h2o_target ){
222 +    
223 +    int n_generated = n_cellX;
224 +    int n_test, nx, ny, nz;
225 +    nx = n_cellX;
226 +    ny = n_cellY;
227 +    nz = n_cellZ;
228 +    
229 +    while( n_test < n_h2o_target ){
230 +      
231 +      nz++;
232 +      n_test = 4 * nx * ny * nz;
233 +    }
234 +    
235 +    int n_diff, goodX, goodY, goodZ;
236 +    
237 +    n_diff = ntest - n_h2o_target;
238 +    goodX = nx;
239 +    goodY = ny;
240 +    goodZ = nz;
241 +    
242 +    int test_diff;
243 +    int n_limit = n_z;
244 +    n_z = n_cellZ;
245 +    
246 +    for( i=n_generated; i<=n_limit; i++ ){
247 +      for( j=i; j<=n_limit; j++ ){
248 +        for( k=j; k<=n_limit; k++ ){
249 +          
250 +          n_test = 4 * i * j * k;
251 +          
252 +          if( n_test > n_h2o_target ){
253 +            
254 +            test_diff = n_test - n_h2o_target;
255 +            
256 +            if( test_diff < n_diff ){
257 +              
258 +              n_diff = test_diff;
259 +              goodX = nx;
260 +              goodY = ny;
261 +              goodZ = nz;
262 +            }
263 +          }
264 +        }
265 +      }
266 +    }
267 +
268 +    n_cellX = goodX;
269 +    n_cellY = goodY;
270 +    n_cellZ = goodZ;
271 +  }
272 +
273 +  // we now have the best box size for the simulation.
274 +    
275 +    
276 +          
277 +
278 +
279 +
280 +
281 +
282 +
283 +
284    int new_nAtoms = group_nAtoms + n_active;
285    Atom** new_atoms = new Atom*[new_nAtoms];
286    
# Line 354 | Line 388 | int main(int argc,char* argv[]){
388    
389    return 0;
390   }
391 +
392 +
393 + void map( x, y, z, centerX, centerY, centerZ, boxX, boxY, boxZ )
394 +     double *x, *y, *z;
395 +     double centerX, centerY, centerZ;
396 +     double boxX, boxY, boxZ;
397 + {
398 +
399 +  *x -= centerX;
400 +  *y -= centerY;
401 +  *z -= centerZ;
402 +
403 +  if(*x < 0) *x -= boxX * (double)( (int)( (*x / boxX)  - 0.5 ) );
404 +  else *x -= boxX * (double)( (int)( (*x / boxX ) + 0.5));
405 +
406 +  if(*y < 0) *y -= boxY * (double)( (int)( (*y / boxY)  - 0.5 ) );
407 +  else *y -= boxY * (double)( (int)( (*y / boxY ) + 0.5));
408 +  
409 +  if(*z < 0) *z -= boxZ * (double)( (int)( (*z / boxZ)  - 0.5 ) );
410 +  else *z -= boxZ * (double)( (int)( (*z / boxZ ) + 0.5));
411 + }

Diff Legend

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