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 32 by mmeineke, Wed Jul 17 20:33:56 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,
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 33 | Line 40 | int main(int argc,char* argv[]){
40    const double water_vol = 4.0 / water_rho; // volume occupied by 4 waters
41    const double water_cell = 4.929; // fcc unit cell length
42  
43 <  int n_lipidsX = 5;
44 <  int n_lipidsY = 10;
45 <  int n_lipids = n_lipidsX * n_lipidsY;
43 >  int n_lipids = 50;
44 >  double water_ratio = 25.0; // water to lipid ratio
45 >  int n_h2o_target = (int)( n_lipids * water_ratio + 0.5 );
46    
47    std::cerr << "n_lipids = " << n_lipids << "\n";
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 67 | Line 74 | int main(int argc,char* argv[]){
74    lipidAtoms = entry_plug->atoms;
75    lipidNAtoms = entry_plug->n_atoms;
76  
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;
77  
78 <  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 <  }
78 >  // find the width, height, and length of the molecule
79  
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 –
80    double min_x, min_y, min_z;
81    double max_x, max_y, max_z;
82    double test_x, test_y, test_z;
83  
84 <  max_x = min_x = group_atoms[0]->getX();
85 <  max_y = min_y = group_atoms[0]->getY();
86 <  max_z = min_z = group_atoms[0]->getZ();
84 >  max_x = min_x = lipidAtoms[0]->getX();
85 >  max_y = min_y = lipidAtoms[0]->getY();
86 >  max_z = min_z = lipidAtoms[0]->getZ();
87    
88 <  for(i=0; i<group_nAtoms; i++){
88 >  for(i=0; i<lipidNAtoms; i++){
89  
90 <    test_x = group_atoms[i]->getX();
91 <    test_y = group_atoms[i]->getY();
92 <    test_z = group_atoms[i]->getZ();
90 >    test_x = lipidAtoms[i]->getX();
91 >    test_y = lipidAtoms[i]->getY();
92 >    test_z = lipidAtoms[i]->getZ();
93  
94      if( test_x < min_x ) min_x = test_x;
95      if( test_y < min_y ) min_y = test_y;
# Line 154 | Line 100 | int main(int argc,char* argv[]){
100      if( test_z > max_z ) max_z = test_z;
101    }
102  
103 <  double box_x = max_x - min_x + 2*water_shell;
104 <  double box_y = max_y - min_y + 2*water_shell;
105 <  double box_z = max_z - min_z + 2*water_shell;
103 >  double ml2 = pow((max_x - min_x), 2 ) + pow((max_y - min_y), 2 )
104 >    + pow((max_x - min_x), 2 );
105 >  double max_length = sqrt( ml2 );
106  
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 );
107    
108 +  // from this information, create the test box
109 +
110 +
111 +  double box_x;
112 +  double box_y;
113 +  double box_z;
114 +
115 +  box_x = box_y = box_z = max_length + water_cell * 4.0; // pad with 4 cells
116 +
117 +  int n_cellX = (int)(box_x / water_cell + 1.0 );
118 +  int n_cellY = (int)(box_y / water_cell + 1.0 );
119 +  int n_cellZ = (int)(box_z / water_cell + 1.0 );
120 +  
121    box_x = water_cell * n_cellX;
122    box_y = water_cell * n_cellY;
123    box_z = water_cell * n_cellZ;
# Line 172 | Line 128 | int main(int argc,char* argv[]){
128    double *waterY = new double[n_water];
129    double *waterZ = new double[n_water];
130    
131 +  
132 +  // find the center of the test lipid, and make it the center of our
133 +  // soon to be created water box.
134 +
135 +
136    double cx, cy, cz;
137  
138    cx = 0.0;
139    cy = 0.0;
140    cz = 0.0;
141 <  for(i=0; i<group_nAtoms; i++){
142 <    cx += group_atoms[i]->getX();
143 <    cy += group_atoms[i]->getY();
144 <    cz += group_atoms[i]->getZ();
141 >  for(i=0; i<lipidNAtoms; i++){
142 >    cx += lipidAtoms[i]->getX();
143 >    cy += lipidAtoms[i]->getY();
144 >    cz += lipidAtoms[i]->getZ();
145    }
146 <  cx /= group_nAtoms;
147 <  cy /= group_nAtoms;
148 <  cz /= group_nAtoms;
146 >  cx /= lipidNAtoms;
147 >  cy /= lipidNAtoms;
148 >  cz /= lipidNAtoms;
149    
150    double x0 = cx - ( box_x * 0.5 );
151    double y0 = cy - ( box_y * 0.5 );
152    double z0 = cz - ( box_z * 0.5 );
153 +
154 +
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[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 + 0.5 * water_cell + x0;
174 <        waterY[index] = j * water_cell + 0.5 * water_cell + y0;
175 <        waterZ[index] = k * 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 + x0;
179 <        waterY[index] = j * water_cell + 0.5 * water_cell + y0;
180 <        waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
181 <        index++;
212 <        
213 <        waterX[index] = i * water_cell + 0.5 * water_cell + x0;
214 <        waterY[index] = j * water_cell + y0;
215 <        waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
216 <        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    }
185  
186 +
187 +  // calculate the number of water's displaced by our molecule.
188 +
189    int *isActive = new int[n_water];
190    for(i=0; i<n_water; i++) isActive[i] = 1;
191    
192 <  int n_active = n_water;
192 >  int n_deleted = 0;
193    double dx, dy, dz;
194    double dx2, dy2, dz2, dSqr;
195    double rCutSqr = water_padding * water_padding;
196    
197    for(i=0; ( (i<n_water) && isActive[i] ); i++){
198 <    for(j=0; ( (j<group_nAtoms) && isActive[i] ); j++){
198 >    for(j=0; ( (j<lipidNAtoms) && isActive[i] ); j++){
199  
200 <      dx = waterX[i] - group_atoms[j]->getX();
201 <      dy = waterY[i] - group_atoms[j]->getY();
202 <      dz = waterZ[i] - group_atoms[j]->getZ();
200 >      dx = waterX[i] - lipidAtoms[j]->getX();
201 >      dy = waterY[i] - lipidAtoms[j]->getY();
202 >      dz = waterZ[i] - lipidAtoms[j]->getZ();
203  
204 +      map( dx, dy, dz, box_x, box_y, box_z );
205 +
206        dx2 = dx * dx;
207        dy2 = dy * dy;
208        dz2 = dz * dz;
# Line 240 | Line 210 | int main(int argc,char* argv[]){
210        dSqr = dx2 + dy2 + dz2;
211        if( dSqr < rCutSqr ){
212          isActive[i] = 0;
213 +        n_deleted++;
214 +      }
215 +    }
216 +  }
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_water < n_h2o_target ){
230 +    
231 +    int n_generated = n_cellX;
232 +    int n_test, nx, ny, nz;
233 +    nx = n_cellX;
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++;
242 +      n_test = 4 * nx * ny * nz;
243 +    }
244 +    
245 +    int n_diff, goodX, goodY, goodZ;
246 +    
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 = nz;
254 +    nz = n_cellZ;
255 +    
256 +    for( i=n_generated; i<=n_limit; i++ ){
257 +      for( j=i; j<=n_limit; j++ ){
258 +        for( k=j; k<=n_limit; k++ ){
259 +          
260 +          n_test = 4 * i * j * k;
261 +          
262 +          if( n_test > n_h2o_target ){
263 +            
264 +            test_diff = n_test - n_h2o_target;
265 +            
266 +            if( test_diff < n_diff ){
267 +              
268 +              n_diff = test_diff;
269 +              goodX = nx;
270 +              goodY = ny;
271 +              goodZ = nz;
272 +            }
273 +          }
274 +        }
275 +      }
276 +    }
277 +
278 +    n_cellX = goodX;
279 +    n_cellY = goodY;
280 +    n_cellZ = goodZ;
281 +  }
282 +
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      }
# Line 247 | Line 491 | int main(int argc,char* argv[]){
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  
286 –
287 –
288 –  
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 354 | Line 598 | int main(int argc,char* argv[]){
598    
599    return 0;
600   }
601 +
602 +
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 +  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 +
616 +
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 +  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