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 33 by mmeineke, Wed Jul 17 20:36:25 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[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    }
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;
209 >      
210 >      dSqr = dx2 + dy2 + dz2;
211 >      if( dSqr < rCutSqr ){
212 >        isActive[i] = 0;
213 >        n_deleted++;
214 >      }
215 >    }
216 >  }
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_water < n_h2o_target ){
225 >    
226 >    int n_generated = n_cellX;
227 >    int n_test, nx, ny, nz;
228 >    nx = n_cellX;
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++;
237 >      n_test = 4 * nx * ny * nz;
238 >    }
239 >    
240 >    int n_diff, goodX, goodY, goodZ;
241 >    
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 = nz;
249 >    nz = n_cellZ;
250 >    
251 >    for( i=n_generated; i<=n_limit; i++ ){
252 >      for( j=i; j<=n_limit; j++ ){
253 >        for( k=j; k<=n_limit; k++ ){
254 >          
255 >          n_test = 4 * i * j * k;
256 >          
257 >          if( n_test > n_h2o_target ){
258 >            
259 >            test_diff = n_test - n_h2o_target;
260 >            
261 >            if( test_diff < n_diff ){
262 >              
263 >              n_diff = test_diff;
264 >              goodX = nx;
265 >              goodY = ny;
266 >              goodZ = nz;
267 >            }
268 >          }
269 >        }
270 >      }
271 >    }
272 >
273 >    n_cellX = goodX;
274 >    n_cellY = goodY;
275 >    n_cellZ = goodZ;
276 >  }
277 >
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 +      map( dx, dy, dz, box_x, box_y, box_z );
472 +
473        dx2 = dx * dx;
474        dy2 = dy * dy;
475        dz2 = dz * dz;
# Line 247 | Line 484 | int main(int argc,char* argv[]){
484  
485    std::cerr << "final n_waters = " << n_active << "\n";
486  
487 <  int new_nAtoms = group_nAtoms + n_active;
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  
286 –
287 –
288 –  
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 354 | Line 591 | int main(int argc,char* argv[]){
591    
592    return 0;
593   }
594 +
595 +
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 +  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 +
609 +
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 +  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