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 `
#	Content
1	#include <iostream>
2	#include <cstdlib>
3	#include <cstring>
4	#include <cstdio>
5	#include <cmath>
6
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	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
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	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
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	// find the width, height, and length of the molecule
76
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 = 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<lipidNAtoms; i++){
86
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;
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	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
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;
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
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<lipidNAtoms; i++){
139	cx += lipidAtoms[i]->getX();
140	cy += lipidAtoms[i]->getY();
141	cz += lipidAtoms[i]->getZ();
142	}
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++ ){
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
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_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<lipidNAtoms) && isActive[i] ); j++){
196
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;
206
207	dSqr = dx2 + dy2 + dz2;
208	if( dSqr < rCutSqr ){
209	isActive[i] = 0;
210	n_deleted++;
211	}
212	}
213	}
214
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
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
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	}