trunk/makeLipid/makeRandomLipid.cpp

#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>

#include "SimSetup.hpp"
#include "SimInfo.hpp"
#include "Atom.hpp"
#include "Integrator.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"


void map( double *x, double *y, double *z, double centerX, double centerY, 
          double centerZ, double boxX, double boxY, double boxZ );

char* program_name;
using namespace std;

int main(int argc,char* argv[]){
  
  int i, j, k, l;
  unsigned int n_atoms, eo, xo;
  char* in_name;
  SimSetup* startMe;
  SimInfo* entry_plug;
  Thermo* tStats;
  
  int lipidNAtoms;
  Atom** lipidAtoms;
  
  int tot_Natoms;
  Atom** totAtoms;
  
  const double water_rho = 0.0334; // number density per cubic angstrom
  const double water_vol = 4.0 / water_rho; // volume occupied by 4 waters
  const double water_cell = 4.929; // fcc unit cell length

  int n_lipids = 50;
  double water_ratio = 25.0; // water to lipid ratio
  int n_h2o_target = (int)( n_lipids * water_ratio + 0.5 );
  
  std::cerr << "n_lipids = " << n_lipids << "\n";

  double water_shell = 10.0;
  double water_padding = 2.5;
  double lipid_spaceing = 4.0;
  
  srand48( 1337 ); // initialize the random number generator.
 
  program_name = argv[0]; /*save the program name in case we need it*/
  if( argc < 3 ){
    cerr<< "Error, input and output bass files are needed to run.\n"
        << program_name << " <input.bass> <output.bass>\n"; 
    
    exit(8);
  }

  in_name = argv[1];
  char* out_name = argv[2];
  entry_plug = new SimInfo;
  
  startMe = new SimSetup;
  startMe->setSimInfo( entry_plug );
  startMe->parseFile( in_name );
  startMe->createSim();

  delete startMe;

  lipidAtoms = entry_plug->atoms;
  lipidNAtoms = entry_plug->n_atoms;


  // find the width, height, and length of the molecule

  double min_x, min_y, min_z;
  double max_x, max_y, max_z;
  double test_x, test_y, test_z;

  max_x = min_x = lipidAtoms[0]->getX();
  max_y = min_y = lipidAtoms[0]->getY();
  max_z = min_z = lipidAtoms[0]->getZ();
  
  for(i=0; i<lipidNAtoms; i++){

    test_x = lipidAtoms[i]->getX();
    test_y = lipidAtoms[i]->getY();
    test_z = lipidAtoms[i]->getZ();

    if( test_x < min_x ) min_x = test_x;
    if( test_y < min_y ) min_y = test_y;
    if( test_z < min_z ) min_z = test_z;

    if( test_x > max_x ) max_x = test_x;
    if( test_y > max_y ) max_y = test_y;
    if( test_z > max_z ) max_z = test_z;
  }

  double ml2 = pow((max_x - min_x), 2 ) + pow((max_y - min_y), 2 )
    + pow((max_x - min_x), 2 );
  double max_length = sqrt( ml2 );

  
  // from this information, create the test box


  double box_x;
  double box_y;
  double box_z;

  box_x = box_y = box_z = max_length + water_cell * 4.0; // pad with 4 cells

  int n_cellX = (int)(box_x / water_cell + 1.0 );
  int n_cellY = (int)(box_y / water_cell + 1.0 );
  int n_cellZ = (int)(box_z / water_cell + 1.0 );
  
  box_x = water_cell * n_cellX;
  box_y = water_cell * n_cellY;
  box_z = water_cell * n_cellZ;

  int n_water = n_cellX * n_cellY * n_cellZ * 4;
  
  double *waterX = new double[n_water];
  double *waterY = new double[n_water];
  double *waterZ = new double[n_water];
  
  
  // find the center of the test lipid, and make it the center of our
  // soon to be created water box.


  double cx, cy, cz;

  cx = 0.0;
  cy = 0.0;
  cz = 0.0;
  for(i=0; i<lipidNAtoms; i++){
    cx += lipidAtoms[i]->getX();
    cy += lipidAtoms[i]->getY();
    cz += lipidAtoms[i]->getZ();
  }
  cx /= lipidNAtoms;
  cy /= lipidNAtoms;
  cz /= lipidNAtoms;
   
  double x0 = cx - ( box_x * 0.5 );
  double y0 = cy - ( box_y * 0.5 );
  double z0 = cz - ( box_z * 0.5 );


  // create an fcc lattice in the water box.

  
  index = 0;
  for( i=0; i < n_cellX; i++ ){
    for( j=0; j < n_cellY; j++ ){
      for( k=0; k < n_cellZ; k++ ){
        
        waterX[index] = i * water_cell + x0; 
        waterY[index] = j * water_cell + y0;
        waterZ[index] = k * water_cell + z0;
        index++;
        
        waterX[index] = i * water_cell + 0.5 * water_cell + x0;
        waterY[index] = j * water_cell + 0.5 * water_cell + y0;
        waterZ[index] = k * water_cell + z0;
        index++;
        
        waterX[index] = i * water_cell + x0;
        waterY[index] = j * water_cell + 0.5 * water_cell + y0;
        waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
        index++;
        
        waterX[index] = i * water_cell + 0.5 * water_cell + x0;
        waterY[index] = j * water_cell + y0;
        waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
        index++;
      }
    }
  }


  // calculate the number of water's displaced by our molecule.

  int *isActive = new int[n_water];
  for(i=0; i<n_water; i++) isActive[i] = 1;
  
  int n_deleted = 0;
  double dx, dy, dz;
  double dx2, dy2, dz2, dSqr;
  double rCutSqr = water_padding * water_padding;
  
  for(i=0; ( (i<n_water) && isActive[i] ); i++){
    for(j=0; ( (j<lipidNAtoms) && isActive[i] ); j++){

      dx = waterX[i] - lipidAtoms[j]->getX();
      dy = waterY[i] - lipidAtoms[j]->getY();
      dz = waterZ[i] - lipidAtoms[j]->getZ();

      map( &dx, &dy, &dz, cx, cy, cz, box_x, box_y, box_z );

      dx2 = dx * dx;
      dy2 = dy * dy;
      dz2 = dz * dz;
      
      dSqr = dx2 + dy2 + dz2;
      if( dSqr < rCutSqr ){
        isActive[i] = 0;
        n_deleted++;
      }
    }
  }

  n_h2o_target += n_deleted * n_lipids;

  // find a box size that best suits the number of waters we need.

  int done = 0;
  
  if( n_waters < n_h2o_target ){
    
    int n_generated = n_cellX;
    int n_test, nx, ny, nz;
    nx = n_cellX;
    ny = n_cellY;
    nz = n_cellZ;
    
    while( n_test < n_h2o_target ){
      
      nz++;
      n_test = 4 * nx * ny * nz;
    }
    
    int n_diff, goodX, goodY, goodZ;
    
    n_diff = ntest - n_h2o_target;
    goodX = nx;
    goodY = ny;
    goodZ = nz;
    
    int test_diff;
    int n_limit = n_z;
    n_z = n_cellZ;
    
    for( i=n_generated; i<=n_limit; i++ ){
      for( j=i; j<=n_limit; j++ ){
        for( k=j; k<=n_limit; k++ ){
          
          n_test = 4 * i * j * k;
          
          if( n_test > n_h2o_target ){
            
            test_diff = n_test - n_h2o_target;
            
            if( test_diff < n_diff ){
              
              n_diff = test_diff;
              goodX = nx;
              goodY = ny;
              goodZ = nz;
            }
          }
        }
      }
    }

    n_cellX = goodX;
    n_cellY = goodY;
    n_cellZ = goodZ;
  }

  // we now have the best box size for the simulation.
    
    
  int new_nAtoms = group_nAtoms + n_active;
  Atom** new_atoms = new Atom*[new_nAtoms];
  
  index = 0;
  for(i=0; i<group_nAtoms; i++ ){
    
    if( group_atoms[i]->isDirectional() ){
      dAtom = (DirectionalAtom *)group_atoms[i];
      
      dAtomNew = new DirectionalAtom();
      dAtomNew->setSUx( dAtom->getSUx() );
      dAtomNew->setSUx( dAtom->getSUx() );
      dAtomNew->setSUx( dAtom->getSUx() );
      
      dAtomNew->setA( rotMat );
      
      new_atoms[index] = dAtomNew;
    }
    else{
      
      new_atoms[index] = new GeneralAtom();
    }

    new_atoms[index]->setType( group_atoms[i]->getType() );    
    
    new_atoms[index]->setX( group_atoms[i]->getX() );    
    new_atoms[index]->setY( group_atoms[i]->getY() );
    new_atoms[index]->setZ( group_atoms[i]->getZ() );

    new_atoms[index]->set_vx( 0.0 );
    new_atoms[index]->set_vy( 0.0 );
    new_atoms[index]->set_vz( 0.0 );
    
    index++;
  }


  for(i=0; i<n_water; i++){
    if(isActive[i]){
      
      new_atoms[index] = new DirectionalAtom();
      new_atoms[index]->setType( "SSD" );

      new_atoms[index]->setX( waterX[i] );    
      new_atoms[index]->setY( waterY[i] );
      new_atoms[index]->setZ( waterZ[i] );
      
      new_atoms[index]->set_vx( 0.0 );
      new_atoms[index]->set_vy( 0.0 );
      new_atoms[index]->set_vz( 0.0 );
      
      dAtom = (DirectionalAtom *) new_atoms[index];

      dAtom->setSUx( 0.0 );
      dAtom->setSUy( 0.0 );
      dAtom->setSUz( 1.0 );
      
      dAtom->setA( rotMat );
      
      index++;
    }
  }

  entry_plug->n_atoms = new_nAtoms;
  entry_plug->atoms = new_atoms;
  
  entry_plug->box_x = box_x;
  entry_plug->box_y = box_y;
  entry_plug->box_z = box_z;

  DumpWriter* xyz_out = new DumpWriter( entry_plug );
  xyz_out->writeFinal();
  delete xyz_out;

  FILE* out_file;
  
  out_file = fopen( out_name, "w" );
  
  fprintf(out_file,
          "#include \"water.mdl\"\n"
          "#include \"lipid.mdl\"\n"
          "\n"
          "nComponents = 2;\n"
          "component{\n"
          "  type = \"theLipid\";\n"
          "  nMol = %d;\n"
          "}\n"
          "\n"
          "component{\n"
          "  type = \"SSD_water\";\n"
          "  nMol = %d;\n"
          "}\n"
          "\n"
          "initialConfig = \"%s\";\n"
          "\n"
          "boxX = %lf;\n"
          "boxY = %lf;\n"
          "boxZ = %lf;\n",
          n_lipids, n_active, entry_plug->finalName,
          box_x, box_y, box_z );
 
  fclose( out_file );
  
  return 0;
}


void map( x, y, z, centerX, centerY, centerZ, boxX, boxY, boxZ )
     double *x, *y, *z;
     double centerX, centerY, centerZ;
     double boxX, boxY, boxZ;
{ 

  *x -= centerX;
  *y -= centerY;
  *z -= centerZ;

  if(*x < 0) *x -= boxX * (double)( (int)( (*x / boxX)  - 0.5 ) );
  else *x -= boxX * (double)( (int)( (*x / boxX ) + 0.5));

  if(*y < 0) *y -= boxY * (double)( (int)( (*y / boxY)  - 0.5 ) );
  else *y -= boxY * (double)( (int)( (*y / boxY ) + 0.5));
  
  if(*z < 0) *z -= boxZ * (double)( (int)( (*z / boxZ)  - 0.5 ) );
  else *z -= boxZ * (double)( (int)( (*z / boxZ ) + 0.5));
}
Revision:	29
Committed:	Tue Jul 16 01:22:04 2002 UTC (21 years, 11 months ago) by mmeineke
File size:	9271 byte(s)
Log Message:	adding the routines to the proogram `
#	User	Rev	Content
1	mmeineke	28	#include <iostream>
2			#include <cstdlib>
3			#include <cstring>
4			#include <cstdio>
5	mmeineke	29	#include <cmath>
6	mmeineke	28
7			#include "SimSetup.hpp"
8			#include "SimInfo.hpp"
9			#include "Atom.hpp"
10			#include "Integrator.hpp"
11			#include "Thermo.hpp"
12			#include "ReadWrite.hpp"
13
14
15	mmeineke	29	void map( double x, double y, double *z, double centerX, double centerY,
16			double centerZ, double boxX, double boxY, double boxZ );
17
18	mmeineke	28	char* program_name;
19			using namespace std;
20
21			int main(int argc,char* argv[]){
22
23			int i, j, k, l;
24			unsigned int n_atoms, eo, xo;
25			char* in_name;
26			SimSetup* startMe;
27			SimInfo* entry_plug;
28			Thermo* tStats;
29
30			int lipidNAtoms;
31			Atom** lipidAtoms;
32
33			int tot_Natoms;
34			Atom** totAtoms;
35
36			const double water_rho = 0.0334; // number density per cubic angstrom
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	mmeineke	29	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	mmeineke	28
44			std::cerr << "n_lipids = " << n_lipids << "\n";
45
46			double water_shell = 10.0;
47			double water_padding = 2.5;
48			double lipid_spaceing = 4.0;
49
50			srand48( 1337 ); // initialize the random number generator.
51
52			program_name = argv[0]; /save the program name in case we need it/
53			if( argc < 3 ){
54			cerr<< "Error, input and output bass files are needed to run.\n"
55			<< program_name << " <input.bass> <output.bass>\n";
56
57			exit(8);
58			}
59
60			in_name = argv[1];
61			char* out_name = argv[2];
62			entry_plug = new SimInfo;
63
64			startMe = new SimSetup;
65			startMe->setSimInfo( entry_plug );
66			startMe->parseFile( in_name );
67			startMe->createSim();
68
69			delete startMe;
70
71			lipidAtoms = entry_plug->atoms;
72			lipidNAtoms = entry_plug->n_atoms;
73
74
75	mmeineke	29	// find the width, height, and length of the molecule
76	mmeineke	28
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	mmeineke	29	max_x = min_x = lipidAtoms[0]->getX();
82			max_y = min_y = lipidAtoms[0]->getY();
83			max_z = min_z = lipidAtoms[0]->getZ();
84	mmeineke	28
85	mmeineke	29	for(i=0; i<lipidNAtoms; i++){
86	mmeineke	28
87	mmeineke	29	test_x = lipidAtoms[i]->getX();
88			test_y = lipidAtoms[i]->getY();
89			test_z = lipidAtoms[i]->getZ();
90	mmeineke	28
91			if( test_x < min_x ) min_x = test_x;
92			if( test_y < min_y ) min_y = test_y;
93			if( test_z < min_z ) min_z = test_z;
94
95			if( test_x > max_x ) max_x = test_x;
96			if( test_y > max_y ) max_y = test_y;
97			if( test_z > max_z ) max_z = test_z;
98			}
99
100	mmeineke	29	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	mmeineke	28
104
105	mmeineke	29	// 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	mmeineke	28	box_x = water_cell * n_cellX;
119			box_y = water_cell * n_cellY;
120			box_z = water_cell * n_cellZ;
121
122			int n_water = n_cellX * n_cellY * n_cellZ * 4;
123
124			double *waterX = new double[n_water];
125			double *waterY = new double[n_water];
126			double *waterZ = new double[n_water];
127
128	mmeineke	29
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	mmeineke	28	double cx, cy, cz;
134
135			cx = 0.0;
136			cy = 0.0;
137			cz = 0.0;
138	mmeineke	29	for(i=0; i<lipidNAtoms; i++){
139			cx += lipidAtoms[i]->getX();
140			cy += lipidAtoms[i]->getY();
141			cz += lipidAtoms[i]->getZ();
142	mmeineke	28	}
143	mmeineke	29	cx /= lipidNAtoms;
144			cy /= lipidNAtoms;
145			cz /= lipidNAtoms;
146	mmeineke	28
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	mmeineke	29
151
152			// create an fcc lattice in the water box.
153
154	mmeineke	28
155			index = 0;
156			for( i=0; i < n_cellX; i++ ){
157			for( j=0; j < n_cellY; j++ ){
158			for( k=0; k < n_cellZ; k++ ){
159
160			waterX[index] = i * water_cell + x0;
161			waterY[index] = j * water_cell + y0;
162			waterZ[index] = k * water_cell + z0;
163			index++;
164
165			waterX[index] = i * water_cell + 0.5 * water_cell + x0;
166			waterY[index] = j * water_cell + 0.5 * water_cell + y0;
167			waterZ[index] = k * water_cell + z0;
168			index++;
169
170			waterX[index] = i * water_cell + x0;
171			waterY[index] = j * water_cell + 0.5 * water_cell + y0;
172			waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
173			index++;
174
175			waterX[index] = i * water_cell + 0.5 * water_cell + x0;
176			waterY[index] = j * water_cell + y0;
177			waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
178			index++;
179			}
180			}
181			}
182
183	mmeineke	29
184			// calculate the number of water's displaced by our molecule.
185
186	mmeineke	28	int *isActive = new int[n_water];
187			for(i=0; i<n_water; i++) isActive[i] = 1;
188
189	mmeineke	29	int n_deleted = 0;
190	mmeineke	28	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	mmeineke	29	for(j=0; ( (j<lipidNAtoms) && isActive[i] ); j++){
196	mmeineke	28
197	mmeineke	29	dx = waterX[i] - lipidAtoms[j]->getX();
198			dy = waterY[i] - lipidAtoms[j]->getY();
199			dz = waterZ[i] - lipidAtoms[j]->getZ();
200	mmeineke	28
201	mmeineke	29	map( &dx, &dy, &dz, cx, cy, cz, box_x, box_y, box_z );
202
203	mmeineke	28	dx2 = dx * dx;
204			dy2 = dy * dy;
205			dz2 = dz * dz;
206
207			dSqr = dx2 + dy2 + dz2;
208			if( dSqr < rCutSqr ){
209			isActive[i] = 0;
210	mmeineke	29	n_deleted++;
211	mmeineke	28	}
212			}
213			}
214
215	mmeineke	29	n_h2o_target += n_deleted * n_lipids;
216	mmeineke	28
217	mmeineke	29	// 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	mmeineke	28	int new_nAtoms = group_nAtoms + n_active;
285			Atom** new_atoms = new Atom*[new_nAtoms];
286
287			index = 0;
288			for(i=0; i<group_nAtoms; i++ ){
289
290			if( group_atoms[i]->isDirectional() ){
291			dAtom = (DirectionalAtom *)group_atoms[i];
292
293			dAtomNew = new DirectionalAtom();
294			dAtomNew->setSUx( dAtom->getSUx() );
295			dAtomNew->setSUx( dAtom->getSUx() );
296			dAtomNew->setSUx( dAtom->getSUx() );
297
298			dAtomNew->setA( rotMat );
299
300			new_atoms[index] = dAtomNew;
301			}
302			else{
303
304			new_atoms[index] = new GeneralAtom();
305			}
306
307			new_atoms[index]->setType( group_atoms[i]->getType() );
308
309			new_atoms[index]->setX( group_atoms[i]->getX() );
310			new_atoms[index]->setY( group_atoms[i]->getY() );
311			new_atoms[index]->setZ( group_atoms[i]->getZ() );
312
313			new_atoms[index]->set_vx( 0.0 );
314			new_atoms[index]->set_vy( 0.0 );
315			new_atoms[index]->set_vz( 0.0 );
316
317			index++;
318			}
319
320
321
322
323			for(i=0; i<n_water; i++){
324			if(isActive[i]){
325
326			new_atoms[index] = new DirectionalAtom();
327			new_atoms[index]->setType( "SSD" );
328
329			new_atoms[index]->setX( waterX[i] );
330			new_atoms[index]->setY( waterY[i] );
331			new_atoms[index]->setZ( waterZ[i] );
332
333			new_atoms[index]->set_vx( 0.0 );
334			new_atoms[index]->set_vy( 0.0 );
335			new_atoms[index]->set_vz( 0.0 );
336
337			dAtom = (DirectionalAtom *) new_atoms[index];
338
339			dAtom->setSUx( 0.0 );
340			dAtom->setSUy( 0.0 );
341			dAtom->setSUz( 1.0 );
342
343			dAtom->setA( rotMat );
344
345			index++;
346			}
347			}
348
349			entry_plug->n_atoms = new_nAtoms;
350			entry_plug->atoms = new_atoms;
351
352			entry_plug->box_x = box_x;
353			entry_plug->box_y = box_y;
354			entry_plug->box_z = box_z;
355
356			DumpWriter* xyz_out = new DumpWriter( entry_plug );
357			xyz_out->writeFinal();
358			delete xyz_out;
359
360			FILE* out_file;
361
362			out_file = fopen( out_name, "w" );
363
364			fprintf(out_file,
365			"#include \"water.mdl\"\n"
366			"#include \"lipid.mdl\"\n"
367			"\n"
368			"nComponents = 2;\n"
369			"component{\n"
370			" type = \"theLipid\";\n"
371			" nMol = %d;\n"
372			"}\n"
373			"\n"
374			"component{\n"
375			" type = \"SSD_water\";\n"
376			" nMol = %d;\n"
377			"}\n"
378			"\n"
379			"initialConfig = \"%s\";\n"
380			"\n"
381			"boxX = %lf;\n"
382			"boxY = %lf;\n"
383			"boxZ = %lf;\n",
384			n_lipids, n_active, entry_plug->finalName,
385			box_x, box_y, box_z );
386
387			fclose( out_file );
388
389			return 0;
390			}
391	mmeineke	29
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			}