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

Comparing trunk/OOPSE/utils/bilayerSys.cpp (file contents):
Revision 498 by mmeineke, Mon Apr 14 21:22:54 2003 UTC vs.
Revision 536 by mmeineke, Fri May 16 14:28:27 2003 UTC

# Line 1 | Line 1
1 + #include <iostream>
2 +
3   #include <cstdlib>
4   #include <cstring>
5   #include <cmath>
6  
7   #include "simError.h"
6 – #include "parse_me.h"
7 – #include "MakeStamps.hpp"
8 – #include "Globals.hpp"
8   #include "SimInfo.hpp"
9   #include "ReadWrite.hpp"
10  
11 <
11 > #include "MoLocator.hpp"
12   #include "sysBuild.hpp"
13   #include "bilayerSys.hpp"
14  
15  
17 – // this routine is defined in BASS_interface.cpp
18 – extern void set_interface_stamps( MakeStamps* ms, Globals* g );
16  
17 < void buildRandomBilayer( sysBuildInfo info );
17 > void map( double &x, double &y, double &z,
18 >          double boxX, double boxY, double boxZ );
19  
20 + int buildRandomBilayer( void );
21  
22 < void buildBilayer( sysBuildInfo info ){
22 > void getRandomRot( double rot[3][3] );
23 >
24 > int buildBilayer( int isRandom ){
25    
26 <  if( info.isRandom ){
27 <    buildRandomBilayer( info );
26 >  if( isRandom ){
27 >    return buildRandomBilayer();
28    }
29    else{
29 –
30      sprintf( painCave.errMsg,
31               "Cannot currently create a non-random bilayer.\n" );
32      painCave.isFatal = 1;
33      simError();
34 +    return 0;
35    }
36   }
37  
38  
39 + int buildRandomBilayer( void ){
40  
41 < void buildRandomBilayer( sysBuildInfo info ){
41 >  typedef struct{
42 >    double rot[3][3];
43 >    double pos[3];
44 >  } coord;
45  
46 <  MakeStamps* the_stamps;
47 <  Globals* the_globals;
48 <  SimInfo* simnfo;
49 <  bwMolLinked bwInfo;
46 >
47 >  const double waterRho = 0.0334; // number density per cubic angstrom
48 >  const double waterVol = 4.0 / water_rho; // volume occupied by 4 waters
49 >  const double waterCell = 4.929; // fcc unit cell length
50 >
51 >  const double water_padding = 2.5;
52 >  const double lipid_spaceing = 2.5;
53 >
54 >
55 >  int i,j,k;
56 >  int nAtoms, atomIndex, molIndex, molID;
57 >  int* molSeq;
58 >  int* molMap;
59 >  int* molStart;
60 >  int* cardDeck;
61 >  int deckSize;
62 >  int rSite, rCard;
63 >  double cell;
64 >  int nCells, nSites, siteIndex;
65    
66 <  // init the bwInfo
66 >  coord *siteArray;
67 >  coord testSite;
68 >
69 >  MoleculeStamp* lipidStamp;
70 >  MoleculeStamp* waterStamp;  
71 >  MoLocator *lipidLocate;
72 >  MoLocator *waterLocate
73 >  int foundLipid, foundWater;
74 >  int nLipids, lipiNatoms, nWaters;
75 >  double testBox, maxLength;
76    
77 <  bwinfo.components = NULL;
77 >  srand48( RAND_SEED );
78 >
79 >
80 >  // set the the lipidStamp
81 >
82 >  foundLipid = 0;
83 >  for(i=0; i<bsInfo.nComponents; i++){
84 >    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.lipidName ) ){
85 >      
86 >      foundlipid = 1;
87 >      lipidStamp = bsInfo.compStamps[i];
88 >      nLipids = bsInfo.componentsNmol[i];
89 >    }
90 >    if( !strcmp( bsInfo.compStamps[i]->getID(), bsInfo.waterName ) ){
91 >      
92 >      foundWater = 1;
93 >      waterStamp = bsInfo.compStamps[i];
94 >      nWaters = bsInfo.componentsNmol[i];
95 >    }
96 >  }
97 >  if( !foundLipid ){
98 >    sprintf(painCave.errMsg,
99 >            "Could not find lipid \"%s\" in the bass file.\n",
100 >            bsInfo.lipidName );
101 >    painCave.isFatal = 1;
102 >    simError();
103 >  }
104 >  if( !foundWater ){
105 >    sprintf(painCave.errMsg,
106 >            "Could not find water \"%s\" in the bass file.\n",
107 >            bsInfo.waterName );
108 >    painCave.isFatal = 1;
109 >    simError();
110 >  }
111    
112 <  bwInfo.havePressure = 0;
51 <  bwInfo.haveTauBarrostat = 0;
52 <  bwInfo.haveTauTemp = 0;
53 <  bwInfo.haveQmass = 0;
112 >  //create the temp Molocator and atom Arrays
113    
114 +  lipidLocate = new MoLocator( lipidStamp );
115 +  lipidNatoms = lipidStamp->getNAtoms();
116 +  maxLength = lipidLocate->getMaxLength();
117 +
118 +  waterLocate = new MoLocator( waterStamp );
119    
120 <  // create parser and read the Bas file
57 <  
58 <  simnfo = new SimInfo();
59 <  the_stamps = new MakeStamps();
60 <  the_globals = new Globals();
61 <  set_interface_stamps( the_stamps, the_globals );
120 >  nAtoms = nLipids * lipidNatoms;
121  
122 <  yacc_BASS( info.in_name );  
122 >  Atom::createArrays( nAtoms );
123 >  atoms = new Atom*[nAtoms];
124  
125 +  // create the test box for initial water displacement
126  
127 +  testBox = maxLength + waterCell * 4.0; // pad with 4 cells
128 +  int nCells = (int)( testBox / waterCell + 1.0 );
129 +  int testWaters = 4 * nCells * nCells * nCells;
130 +  
131 +  double* waterX = new double[testWaters];
132 +  double* waterX = new double[testWaters];
133 +  double* waterX = new double[testWaters];
134 +  
135 +  double x0 = 0.0 - ( testBox * 0.5 );
136 +  double y0 = 0.0 - ( testBox * 0.5 );
137 +  double z0 = 0.0 - ( testBox * 0.5 );
138  
67 –  // set the easy ones first
68 –  bwInfo.targetTemp = the_globals->getTargetTemp();
69 –  bwInfo.dt = the_globals->getDt();
70 –  bwInfo.runTime = the_globals->getRunTime();
139  
140 <  // get the ones we know are there, yet still may need some work.
73 <  bwInfo.nComponents = the_globals->getNComponents();
74 <  strcpy( bwInfo.forceField, the_globals->getForceField() );
140 >  // create an fcc lattice in the water box.
141  
142 <  // get the ensemble:
143 <  strcpy( bwInfo.ensemble, the_globals->getEnsemble() );
144 <  if( !strcasecmp( bwInfo.ensemble, "NPT" ) ) {
145 <    
146 <    if (the_globals->haveTargetPressure()){
147 <      bwInfo.targetPressure = the_globals->getTargetPressure();
148 <      bwInfo.havePressure = 1;
142 >  int ndx = 0;
143 >  for( i=0; i < nCells; i++ ){
144 >    for( j=0; j < nCells; j++ ){
145 >      for( k=0; k < nCells; k++ ){
146 >        
147 >        waterX[ndx] = i * waterCell + x0;
148 >        waterY[ndx] = j * waterCell + y0;
149 >        waterZ[ndx] = k * waterCell + z0;
150 >        ndx++;
151 >        
152 >        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
153 >        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
154 >        waterZ[ndx] = k * waterCell + z0;
155 >        ndx++;
156 >        
157 >        waterX[ndx] = i * waterCell + x0;
158 >        waterY[ndx] = j * waterCell + 0.5 * waterCell + y0;
159 >        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
160 >        ndx++;
161 >        
162 >        waterX[ndx] = i * waterCell + 0.5 * waterCell + x0;
163 >        waterY[ndx] = j * waterCell + y0;
164 >        waterZ[ndx] = k * waterCell + 0.5 * waterCell + z0;
165 >        ndx++;
166 >      }
167      }
168 <    else {
85 <      sprintf( painCave.errMsg,
86 <               "buildBilayer error: If you use the constant pressure\n"
87 <               "    ensemble, you must set targetPressure.\n"
88 <               "    This was found in the BASS file.\n");
89 <      painCave.isFatal = 1;
90 <      simError();
91 <    }
168 >  }
169  
170 <    if (the_globals->haveTauThermostat()){
94 <      bwInfo.tauThermostat = the_globals->getTauThermostat();
95 <      bwInfo.haveTauThermostat = 1;;
96 <    }
97 <    else if (the_globals->haveQmass()){
98 <      bwinfo.Qmass = the_globals->getQmass();
99 <      bwInfo.haveQmass = 1;
100 <    }
101 <    else {
102 <      sprintf( painCave.errMsg,
103 <               "buildBilayer error: If you use one of the constant temperature\n"
104 <               "    ensembles, you must set either tauThermostat or qMass.\n"
105 <               "    Neither of these was found in the BASS file.\n");
106 <      painCave.isFatal = 1;
107 <      simError();
108 <    }
170 >  // calculate the number of water's displaced by our lipid.
171  
172 <    if (the_globals->haveTauBarostat()){
173 <      bwInfo.tauBarostat = the_globals->getTauBarostat();
174 <      bwInfo.haveTauBarostat = 1;
113 <    }                                    
114 <    else {
115 <      sprintf( painCave.errMsg,
116 <               "SimSetup error: If you use the constant pressure\n"
117 <               "    ensemble, you must set tauBarostat.\n"
118 <               "    This was found in the BASS file.\n");
119 <      painCave.isFatal = 1;
120 <      simError();
121 <    }
172 >  testSite.rot[0][0] = 1.0;
173 >  testSite.rot[0][1] = 0.0;
174 >  testSite.rot[0][2] = 0.0;
175  
176 <  } else if ( !strcasecmp( ensemble, "NVT") ) {
177 <    the_extendedsystem = new ExtendedSystem( simnfo );
178 <    the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
176 >  testSite.rot[1][0] = 0.0;
177 >  testSite.rot[1][1] = 1.0;
178 >  testSite.rot[1][2] = 0.0;
179  
180 <    if (the_globals->haveTauThermostat())
181 <      the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
182 <    else if (the_globals->haveQmass())
183 <      the_extendedsystem->setQmass(the_globals->getQmass());
184 <    else {
185 <      sprintf( painCave.errMsg,
186 <               "SimSetup error: If you use one of the constant temperature\n"
187 <               "    ensembles, you must set either tauThermostat or qMass.\n"
188 <               "    Neither of these was found in the BASS file.\n");
189 <      painCave.isFatal = 1;
190 <      simError();
180 >  testSite.rot[2][0] = 0.0;
181 >  testSite.rot[2][1] = 0.0;
182 >  testSite.rot[2][2] = 1.0;
183 >  
184 >  testSite.pos[0] = 0.0;
185 >  testSite.pos[1] = 0.0;
186 >  testSite.pos[2] = 0.0;
187 >
188 >  lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0 );
189 >
190 >  int *isActive = new int[testWaters];
191 >  for(i=0; i<testWaters; i++) isActive[i] = 1;
192 >  
193 >  int n_deleted = 0;
194 >  double dx, dy, dz;
195 >  double dx2, dy2, dz2, dSqr;
196 >  double rCutSqr = water_padding * water_padding;
197 >  
198 >  for(i=0; ( (i<testWaters) && isActive[i] ); i++){
199 >    for(j=0; ( (j<lipidNatoms) && isActive[i] ); j++){
200 >
201 >      dx = waterX[i] - atoms[j]->getX();
202 >      dy = waterY[i] - atoms[j]->getY();
203 >      dz = waterZ[i] - atoms[j]->getZ();
204 >
205 >      map( dx, dy, dz, testBox, testBox, testBox );
206 >
207 >      dx2 = dx * dx;
208 >      dy2 = dy * dy;
209 >      dz2 = dz * dz;
210 >      
211 >      dSqr = dx2 + dy2 + dz2;
212 >      if( dSqr < rCutSqr ){
213 >        isActive[i] = 0;
214 >        n_deleted++;
215 >      }
216      }
217 +  }
218 +  
219 +  int targetWaters = nWaters + n_deleted * nLipids;
220  
221 <  } else if ( !strcasecmp( ensemble, "NVE") ) {
222 <  } else {
223 <    sprintf( painCave.errMsg,
224 <             "SimSetup Warning. Unrecognized Ensemble -> %s, "
225 <             "reverting to NVE for this simulation.\n",
226 <             ensemble );
227 <    painCave.isFatal = 0;
228 <    simError();
229 <    strcpy( ensemble, "NVE" );
221 >  // find the best box size for the sim
222 >
223 >  int testTot;
224 >  int done = 0;
225 >  ndx = 0;
226 >  while( !done ){
227 >    
228 >    ndx++;
229 >    testTot = 4 * ndx * ndx * ndx;
230 >    
231 >    if( testTot >= targetWaters ) done = 1;
232 >  }
233 >
234 >  nCells = ndx;
235 >    
236 >    
237 >  // create the new water box to the new specifications
238 >  
239 >  int newWaters = nCells * nCells * nCells * 4;
240 >  
241 >  delete[] waterX;
242 >  delete[] waterY;
243 >  delete[] waterZ;
244 >  
245 >  waterX = new double[newWater];
246 >  waterY = new double[newWater];
247 >  waterZ = new double[newWater];
248 >  
249 >  double box_x = waterCell * nCells;
250 >  double box_y = waterCell * nCells;
251 >  double box_z = waterCell * nCells;
252 >    
253 >  // create an fcc lattice in the water box.
254 >  
255 >  ndx = 0;
256 >  for( i=0; i < nCells; i++ ){
257 >    for( j=0; j < nCells; j++ ){
258 >      for( k=0; k < nCells; k++ ){
259 >        
260 >        waterX[ndx] = i * waterCell;
261 >        waterY[ndx] = j * waterCell;
262 >        waterZ[ndx] = k * waterCell;
263 >        ndx++;
264 >        
265 >        waterX[ndx] = i * waterCell + 0.5 * waterCell;
266 >        waterY[ndx] = j * waterCell + 0.5 * waterCell;
267 >        waterZ[ndx] = k * waterCell;
268 >        ndx++;
269 >        
270 >        waterX[ndx] = i * waterCell;
271 >        waterY[ndx] = j * waterCell + 0.5 * waterCell;
272 >        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
273 >        ndx++;
274 >        
275 >        waterX[ndx] = i * waterCell + 0.5 * waterCell;
276 >        waterY[ndx] = j * waterCell;
277 >        waterZ[ndx] = k * waterCell + 0.5 * waterCell;
278 >        ndx++;
279 >      }
280 >    }
281    }  
150 –  strcpy( simnfo->ensemble, ensemble );
282  
283  
284  
285 <  delete stamps;
286 <  delete globals;
285 >
286 >
287 >
288 >
289 >
290 >
291 >
292 >
293 >
294 >
295 >
296 >  // create the real MoLocator and Atom arrays
297 >  
298 >  nAtoms = 0;
299 >  molIndex = 0;
300 >  locate = new MoLocator*[bsInfo.nComponents];
301 >  molSeq = new int[bsInfo.totNmol];
302 >  molStart = new int[bsInfo.totNmol];  
303 >  for(i=0; i<bsInfo.nComponents; i++){
304 >    locate[i] = new MoLocator( bsInfo.compStamps[i] );
305 >    for(j=0; j<bsInfo.componentsNmol[i]; j++){
306 >      molSeq[molIndex] = i;
307 >      molStart[molIndex] = nAtoms;
308 >      molIndex++;
309 >      nAtoms += bsInfo.compStamps[i]->getNAtoms();
310 >    }
311 >  }
312 >  
313 >  Atom::createArrays( nAtoms );
314 >  atoms = new Atom*[nAtoms];
315 >
316 >
317 >  // find the width, height, and length of the molecule
318 >
319 >  
320 >
321 >
322   }
323 +
324 +
325 +
326 + int Old_buildRandomBilayer( void ){
327 +
328 +  int i,j,k;
329 +  int nAtoms, atomIndex, molIndex, molID;
330 +  int* molSeq;
331 +  int* molMap;
332 +  int* molStart;
333 +  int* cardDeck;
334 +  int deckSize;
335 +  int rSite, rCard;
336 +  double cell;
337 +  int nCells, nSites, siteIndex;
338 +  double rot[3][3];
339 +  double pos[3];
340 +
341 +  Atom** atoms;
342 +  SimInfo* simnfo;
343 +  DumpWriter* writer;
344 +  MoLocator** locate;
345 +
346 +  // initialize functions and variables
347 +  
348 +  srand48( RAND_SEED );
349 +  molSeq = NULL;
350 +  molStart = NULL;
351 +  molMap = NULL;
352 +  cardDeck = NULL;
353 +  atoms = NULL;
354 +  locate = NULL;
355 +  simnfo = NULL;
356 +  writer = NULL;
357 +
358 +  // calculate the number of cells in the fcc box
359 +
360 +  nCells = 0;
361 +  nSites = 0;
362 +  while( nSites < bsInfo.totNmol ){
363 +    nCells++;
364 +    nSites = 4.0 * pow( (double)nCells, 3.0 );
365 +  }
366 +
367 +
368 +  // create the molMap and cardDeck arrays
369 +  
370 +  molMap = new int[nSites];
371 +  cardDeck = new int[nSites];
372 +  
373 +  for(i=0; i<nSites; i++){
374 +    molMap[i] = -1;
375 +    cardDeck[i] = i;
376 +  }
377 +
378 +  // randomly place the molecules on the sites
379 +  
380 +  deckSize = nSites;
381 +  for(i=0; i<bsInfo.totNmol; i++){
382 +    rCard = (int)( deckSize * drand48() );
383 +    rSite = cardDeck[rCard];
384 +    molMap[rSite] = i;
385 +
386 +    // book keep the card deck;
387 +    
388 +    deckSize--;
389 +    cardDeck[rCard] = cardDeck[deckSize];
390 +  }
391 +  
392 +  
393 +  // create the MoLocator and Atom arrays
394 +  
395 +  nAtoms = 0;
396 +  molIndex = 0;
397 +  locate = new MoLocator*[bsInfo.nComponents];
398 +  molSeq = new int[bsInfo.totNmol];
399 +  molStart = new int[bsInfo.totNmol];  
400 +  for(i=0; i<bsInfo.nComponents; i++){
401 +    locate[i] = new MoLocator( bsInfo.compStamps[i] );
402 +    for(j=0; j<bsInfo.componentsNmol[i]; j++){
403 +      molSeq[molIndex] = i;
404 +      molStart[molIndex] = nAtoms;
405 +      molIndex++;
406 +      nAtoms += bsInfo.compStamps[i]->getNAtoms();
407 +    }
408 +  }
409 +  
410 +  Atom::createArrays( nAtoms );
411 +  atoms = new Atom*[nAtoms];
412 +  
413 +  
414 +  // place the molecules at each FCC site
415 +  
416 +  cell = 5.0;
417 +  for(i=0; i<bsInfo.nComponents; i++){
418 +    if(cell < locate[i]->getMaxLength() ) cell = locate[i]->getMaxLength();
419 +  }
420 +  cell *= 1.2; // add a little buffer
421 +
422 +  cell *= M_SQRT2;
423 +
424 +  siteIndex = 0;
425 +  for(i=0; i<nCells; i++){
426 +    for(j=0; j<nCells; j++){
427 +      for(k=0; k<nCells; k++){
428 +        
429 +        if( molMap[siteIndex] >= 0 ){
430 +          pos[0] = i * cell;
431 +          pos[1] = j * cell;
432 +          pos[2] = k * cell;
433 +          
434 +          getRandomRot( rot );
435 +          molID = molSeq[molMap[siteIndex]];
436 +          atomIndex = molStart[ molMap[siteIndex] ];
437 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
438 +        }
439 +        siteIndex++;
440 +
441 +        if( molMap[siteIndex] >= 0 ){
442 +          pos[0] = i * cell + (0.5 * cell);
443 +          pos[1] = j * cell;
444 +          pos[2] = k * cell + (0.5 * cell);
445 +
446 +          getRandomRot( rot );
447 +          molID = molSeq[molMap[siteIndex]];
448 +          atomIndex = molStart[ molMap[siteIndex] ];
449 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
450 +        }
451 +        siteIndex++;
452 +
453 +        if( molMap[siteIndex] >= 0 ){
454 +          pos[0] = i * cell + (0.5 * cell);
455 +          pos[1] = j * cell + (0.5 * cell);
456 +          pos[2] = k * cell;
457 +          
458 +          getRandomRot( rot );
459 +          molID = molSeq[molMap[siteIndex]];
460 +          atomIndex = molStart[ molMap[siteIndex] ];
461 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
462 +        }
463 +        siteIndex++;
464 +
465 +        if( molMap[siteIndex] >= 0 ){
466 +          pos[0] = i * cell;
467 +          pos[1] = j * cell + (0.5 * cell);
468 +          pos[2] = k * cell + (0.5 * cell);
469 +
470 +          getRandomRot( rot );
471 +          molID = molSeq[molMap[siteIndex]];
472 +          atomIndex = molStart[ molMap[siteIndex] ];
473 +          locate[molID]->placeMol( pos, rot, atoms, atomIndex );
474 +        }
475 +        siteIndex++;
476 +      }
477 +    }
478 +  }
479 +
480 +  // set up the SimInfo object
481 +  
482 +  bsInfo.boxX = nCells * cell;
483 +  bsInfo.boxY = nCells * cell;
484 +  bsInfo.boxZ = nCells * cell;
485 +  
486 +  simnfo = new SimInfo();
487 +  simnfo->n_atoms = nAtoms;
488 +  simnfo->box_x = bsInfo.boxX;
489 +  simnfo->box_y = bsInfo.boxY;
490 +  simnfo->box_z = bsInfo.boxZ;
491 +  
492 +  sprintf( simnfo->sampleName, "%s.dump", bsInfo.outPrefix );
493 +  sprintf( simnfo->finalName, "%s.init", bsInfo.outPrefix );
494 +
495 +  simnfo->atoms = atoms;
496 +  
497 +  // set up the writer and write out
498 +  
499 +  writer = new DumpWriter( simnfo );
500 +  writer->writeFinal();
501 +    
502 +  // clean up the memory
503 +  
504 +  if( molMap != NULL )   delete[] molMap;
505 +  if( cardDeck != NULL ) delete[] cardDeck;
506 +  if( locate != NULL ){
507 +    for(i=0; i<bsInfo.nComponents; i++){
508 +      delete locate[i];
509 +    }
510 +    delete[] locate;
511 +  }
512 +  if( atoms != NULL ){
513 +    for(i=0; i<nAtoms; i++){
514 +      delete atoms[i];
515 +    }
516 +    Atom::destroyArrays();
517 +    delete[] atoms;
518 +  }
519 +  if( molSeq != NULL ) delete[] molSeq;
520 +  if( simnfo != NULL ) delete simnfo;
521 +  if( writer != NULL ) delete writer;
522 +
523 +  return 1;
524 + }
525 +
526 +
527 + void getRandomRot( double rot[3][3] ){
528 +
529 +  double theta, phi, psi;
530 +  double cosTheta;
531 +
532 +  // select random phi, psi, and cosTheta
533 +
534 +  phi = 2.0 * M_PI * drand48();
535 +  psi = 2.0 * M_PI * drand48();
536 +  cosTheta = (2.0 * drand48()) - 1.0; // sample cos -1 to 1
537 +
538 +  theta = acos( cosTheta );
539 +
540 +  rot[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
541 +  rot[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
542 +  rot[0][2] = sin(theta) * sin(psi);
543 +  
544 +  rot[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
545 +  rot[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
546 +  rot[1][2] = sin(theta) * cos(psi);
547 +
548 +  rot[2][0] = sin(phi) * sin(theta);
549 +  rot[2][1] = -cos(phi) * sin(theta);
550 +  rot[2][2] = cos(theta);  
551 + }
552 +                        
553 +
554 +
555 + void map( double &x, double &y, double &z,
556 +          double boxX, double boxY, double boxZ ){
557 +  
558 +  if(x < 0) x -= boxX * (double)( (int)( (x / boxX)  - 0.5 ) );
559 +  else x -= boxX * (double)( (int)( (x / boxX ) + 0.5));
560 +
561 +  if(y < 0) y -= boxY * (double)( (int)( (y / boxY)  - 0.5 ) );
562 +  else y -= boxY * (double)( (int)( (y / boxY ) + 0.5));
563 +  
564 +  if(z < 0) z -= boxZ * (double)( (int)( (z / boxZ)  - 0.5 ) );
565 +  else z -= boxZ * (double)( (int)( (z / boxZ ) + 0.5));
566 + }

Diff Legend

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