trunk/makeLipid/makeLipid.cpp

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

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


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_lipidsX = 5;
  int n_lipidsY = 10;
  int n_lipids = n_lipidsX * n_lipidsY;
  
  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;

  int group_nAtoms = n_lipids * lipidNAtoms;
  Atom** group_atoms = new Atom*[group_nAtoms];
  DirectionalAtom* dAtom;
  DirectionalAtom* dAtomNew;
  
  double rotMat[3][3];
  double unitRotMat[3][3];
  
  unitRotMat[0][0] = 1.0;
  unitRotMat[0][1] = 0.0;
  unitRotMat[0][2] = 0.0;
  
  unitRotMat[1][0] = 0.0;
  unitRotMat[1][1] = 1.0;
  unitRotMat[1][2] = 0.0;
  
  unitRotMat[2][0] = 0.0;
  unitRotMat[2][1] = 0.0;
  unitRotMat[2][2] = 1.0;

  int index =0;
  for(i=0; i<n_lipids; i++ ){
    for(j=0; j<lipidNAtoms; j++){
      
      if( lipidAtoms[j]->isDirectional() ){
        dAtom = (DirectionalAtom *)lipidAtoms[j];
        
        dAtomNew = new DirectionalAtom();
        dAtomNew->setSUx( dAtom->getSUx() );
        dAtomNew->setSUx( dAtom->getSUx() );
        dAtomNew->setSUx( dAtom->getSUx() );
        
        dAtom->getA( rotMat );
        dAtomNew->setA( rotMat );
        
        group_atoms[index] = dAtomNew;
      }
      else{
        
        group_atoms[index] = new GeneralAtom();
      }
      
      group_atoms[index]->setType( lipidAtoms[j]->getType() );
      
      index++;
    }
  }

  index = 0;
  for(i=0; i<n_lipidsX; i++){
    for(j=0; j<n_lipidsY; j++){
      for(l=0; l<lipidNAtoms; l++){
        
        group_atoms[index]->setX( lipidAtoms[l]->getX() + 
                                  i*lipid_spaceing );
        
        group_atoms[index]->setY( lipidAtoms[l]->getY() + 
                                  j*lipid_spaceing );
        
        group_atoms[index]->setZ( lipidAtoms[l]->getZ() );
        
        index++;
      }
    }
  }

  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 = group_atoms[0]->getX();
  max_y = min_y = group_atoms[0]->getY();
  max_z = min_z = group_atoms[0]->getZ();
  
  for(i=0; i<group_nAtoms; i++){

    test_x = group_atoms[i]->getX();
    test_y = group_atoms[i]->getY();
    test_z = group_atoms[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 box_x = max_x - min_x + 2*water_shell;
  double box_y = max_y - min_y + 2*water_shell;
  double box_z = max_z - min_z + 2*water_shell;

  int n_cellX = (int)(box_x / water_cell + 0.5 );
  int n_cellY = (int)(box_y / water_cell + 0.5 );
  int n_cellZ = (int)(box_z / water_cell + 0.5 );
  
  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];
  
  double cx, cy, cz;

  cx = 0.0;
  cy = 0.0;
  cz = 0.0;
  for(i=0; i<group_nAtoms; i++){
    cx += group_atoms[i]->getX();
    cy += group_atoms[i]->getY();
    cz += group_atoms[i]->getZ();
  }
  cx /= group_nAtoms;
  cy /= group_nAtoms;
  cz /= group_nAtoms;
   
  double x0 = cx - ( box_x * 0.5 );
  double y0 = cy - ( box_y * 0.5 );
  double z0 = cz - ( box_z * 0.5 );
  
  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++;
      }
    }
  }

  int *isActive = new int[n_water];
  for(i=0; i<n_water; i++) isActive[i] = 1;
  
  int n_active = n_water;
  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<group_nAtoms) && isActive[i] ); j++){

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

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

  std::cerr << "final n_waters = " << n_active << "\n";

  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() );
      
      dAtom->getA(rotMat);
      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( unitRotMat );
      
      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;
}
Revision:	32
Committed:	Wed Jul 17 20:33:56 2002 UTC (21 years, 11 months ago) by mmeineke
File size:	8346 byte(s)
Log Message:	fixed makeRandom to set up a correct simulation
#	User	Rev	Content
1	mmeineke	16	#include <iostream>
2			#include <cstdlib>
3			#include <cstring>
4			#include <cstdio>
5
6			#include "SimSetup.hpp"
7			#include "SimInfo.hpp"
8			#include "Atom.hpp"
9			#include "Integrator.hpp"
10			#include "Thermo.hpp"
11			#include "ReadWrite.hpp"
12
13
14			char* program_name;
15			using namespace std;
16
17			int main(int argc,char* argv[]){
18
19			int i, j, k, l;
20			unsigned int n_atoms, eo, xo;
21			char* in_name;
22			SimSetup* startMe;
23			SimInfo* entry_plug;
24			Thermo* tStats;
25
26			int lipidNAtoms;
27			Atom** lipidAtoms;
28
29			int tot_Natoms;
30			Atom** totAtoms;
31
32			const double water_rho = 0.0334; // number density per cubic angstrom
33			const double water_vol = 4.0 / water_rho; // volume occupied by 4 waters
34			const double water_cell = 4.929; // fcc unit cell length
35
36			int n_lipidsX = 5;
37	mmeineke	27	int n_lipidsY = 10;
38	mmeineke	16	int n_lipids = n_lipidsX * n_lipidsY;
39
40			std::cerr << "n_lipids = " << n_lipids << "\n";
41
42			double water_shell = 10.0;
43			double water_padding = 2.5;
44			double lipid_spaceing = 4.0;
45
46			srand48( 1337 ); // initialize the random number generator.
47
48			program_name = argv[0]; /save the program name in case we need it/
49			if( argc < 3 ){
50			cerr<< "Error, input and output bass files are needed to run.\n"
51			<< program_name << " <input.bass> <output.bass>\n";
52
53			exit(8);
54			}
55
56			in_name = argv[1];
57			char* out_name = argv[2];
58			entry_plug = new SimInfo;
59
60			startMe = new SimSetup;
61			startMe->setSimInfo( entry_plug );
62			startMe->parseFile( in_name );
63			startMe->createSim();
64
65			delete startMe;
66
67			lipidAtoms = entry_plug->atoms;
68			lipidNAtoms = entry_plug->n_atoms;
69
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	mmeineke	32	double unitRotMat[3][3];
77	mmeineke	16
78	mmeineke	32	unitRotMat[0][0] = 1.0;
79			unitRotMat[0][1] = 0.0;
80			unitRotMat[0][2] = 0.0;
81	mmeineke	16
82	mmeineke	32	unitRotMat[1][0] = 0.0;
83			unitRotMat[1][1] = 1.0;
84			unitRotMat[1][2] = 0.0;
85	mmeineke	16
86	mmeineke	32	unitRotMat[2][0] = 0.0;
87			unitRotMat[2][1] = 0.0;
88			unitRotMat[2][2] = 1.0;
89	mmeineke	16
90			int index =0;
91			for(i=0; i<n_lipids; i++ ){
92			for(j=0; j<lipidNAtoms; j++){
93
94			if( lipidAtoms[j]->isDirectional() ){
95			dAtom = (DirectionalAtom *)lipidAtoms[j];
96
97			dAtomNew = new DirectionalAtom();
98			dAtomNew->setSUx( dAtom->getSUx() );
99			dAtomNew->setSUx( dAtom->getSUx() );
100			dAtomNew->setSUx( dAtom->getSUx() );
101	mmeineke	32
102			dAtom->getA( rotMat );
103	mmeineke	16	dAtomNew->setA( rotMat );
104
105			group_atoms[index] = dAtomNew;
106			}
107			else{
108
109			group_atoms[index] = new GeneralAtom();
110			}
111
112			group_atoms[index]->setType( lipidAtoms[j]->getType() );
113
114			index++;
115			}
116			}
117
118			index = 0;
119			for(i=0; i<n_lipidsX; i++){
120			for(j=0; j<n_lipidsY; j++){
121			for(l=0; l<lipidNAtoms; l++){
122
123			group_atoms[index]->setX( lipidAtoms[l]->getX() +
124			i*lipid_spaceing );
125
126			group_atoms[index]->setY( lipidAtoms[l]->getY() +
127			j*lipid_spaceing );
128
129			group_atoms[index]->setZ( lipidAtoms[l]->getZ() );
130
131			index++;
132			}
133			}
134			}
135
136			double min_x, min_y, min_z;
137			double max_x, max_y, max_z;
138			double test_x, test_y, test_z;
139
140			max_x = min_x = group_atoms[0]->getX();
141			max_y = min_y = group_atoms[0]->getY();
142			max_z = min_z = group_atoms[0]->getZ();
143
144			for(i=0; i<group_nAtoms; i++){
145
146			test_x = group_atoms[i]->getX();
147			test_y = group_atoms[i]->getY();
148			test_z = group_atoms[i]->getZ();
149
150			if( test_x < min_x ) min_x = test_x;
151			if( test_y < min_y ) min_y = test_y;
152			if( test_z < min_z ) min_z = test_z;
153
154			if( test_x > max_x ) max_x = test_x;
155			if( test_y > max_y ) max_y = test_y;
156			if( test_z > max_z ) max_z = test_z;
157			}
158
159			double box_x = max_x - min_x + 2*water_shell;
160			double box_y = max_y - min_y + 2*water_shell;
161			double box_z = max_z - min_z + 2*water_shell;
162
163			int n_cellX = (int)(box_x / water_cell + 0.5 );
164			int n_cellY = (int)(box_y / water_cell + 0.5 );
165			int n_cellZ = (int)(box_z / water_cell + 0.5 );
166
167			box_x = water_cell * n_cellX;
168			box_y = water_cell * n_cellY;
169			box_z = water_cell * n_cellZ;
170
171			int n_water = n_cellX * n_cellY * n_cellZ * 4;
172
173			double *waterX = new double[n_water];
174			double *waterY = new double[n_water];
175			double *waterZ = new double[n_water];
176
177			double cx, cy, cz;
178
179			cx = 0.0;
180			cy = 0.0;
181			cz = 0.0;
182			for(i=0; i<group_nAtoms; i++){
183			cx += group_atoms[i]->getX();
184			cy += group_atoms[i]->getY();
185			cz += group_atoms[i]->getZ();
186			}
187			cx /= group_nAtoms;
188			cy /= group_nAtoms;
189			cz /= group_nAtoms;
190
191			double x0 = cx - ( box_x * 0.5 );
192			double y0 = cy - ( box_y * 0.5 );
193			double z0 = cz - ( box_z * 0.5 );
194
195			index = 0;
196			for( i=0; i < n_cellX; i++ ){
197			for( j=0; j < n_cellY; j++ ){
198			for( k=0; k < n_cellZ; k++ ){
199
200			waterX[index] = i * water_cell + x0;
201			waterY[index] = j * water_cell + y0;
202			waterZ[index] = k * water_cell + z0;
203			index++;
204
205			waterX[index] = i * water_cell + 0.5 * water_cell + x0;
206			waterY[index] = j * water_cell + 0.5 * water_cell + y0;
207			waterZ[index] = k * water_cell + z0;
208			index++;
209
210			waterX[index] = i * water_cell + x0;
211			waterY[index] = j * water_cell + 0.5 * water_cell + y0;
212			waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
213			index++;
214
215			waterX[index] = i * water_cell + 0.5 * water_cell + x0;
216			waterY[index] = j * water_cell + y0;
217			waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
218			index++;
219			}
220			}
221			}
222
223			int *isActive = new int[n_water];
224			for(i=0; i<n_water; i++) isActive[i] = 1;
225
226			int n_active = n_water;
227			double dx, dy, dz;
228			double dx2, dy2, dz2, dSqr;
229			double rCutSqr = water_padding * water_padding;
230
231			for(i=0; ( (i<n_water) && isActive[i] ); i++){
232			for(j=0; ( (j<group_nAtoms) && isActive[i] ); j++){
233
234			dx = waterX[i] - group_atoms[j]->getX();
235			dy = waterY[i] - group_atoms[j]->getY();
236			dz = waterZ[i] - group_atoms[j]->getZ();
237
238			dx2 = dx * dx;
239			dy2 = dy * dy;
240			dz2 = dz * dz;
241
242			dSqr = dx2 + dy2 + dz2;
243			if( dSqr < rCutSqr ){
244			isActive[i] = 0;
245			n_active--;
246			}
247			}
248			}
249
250			std::cerr << "final n_waters = " << n_active << "\n";
251
252			int new_nAtoms = group_nAtoms + n_active;
253			Atom** new_atoms = new Atom*[new_nAtoms];
254
255			index = 0;
256			for(i=0; i<group_nAtoms; i++ ){
257
258			if( group_atoms[i]->isDirectional() ){
259			dAtom = (DirectionalAtom *)group_atoms[i];
260
261			dAtomNew = new DirectionalAtom();
262			dAtomNew->setSUx( dAtom->getSUx() );
263			dAtomNew->setSUx( dAtom->getSUx() );
264			dAtomNew->setSUx( dAtom->getSUx() );
265
266	mmeineke	32	dAtom->getA(rotMat);
267	mmeineke	16	dAtomNew->setA( rotMat );
268
269			new_atoms[index] = dAtomNew;
270			}
271			else{
272
273			new_atoms[index] = new GeneralAtom();
274			}
275
276			new_atoms[index]->setType( group_atoms[i]->getType() );
277
278			new_atoms[index]->setX( group_atoms[i]->getX() );
279			new_atoms[index]->setY( group_atoms[i]->getY() );
280			new_atoms[index]->setZ( group_atoms[i]->getZ() );
281
282			new_atoms[index]->set_vx( 0.0 );
283			new_atoms[index]->set_vy( 0.0 );
284			new_atoms[index]->set_vz( 0.0 );
285
286			index++;
287			}
288
289
290
291
292			for(i=0; i<n_water; i++){
293			if(isActive[i]){
294
295			new_atoms[index] = new DirectionalAtom();
296			new_atoms[index]->setType( "SSD" );
297
298			new_atoms[index]->setX( waterX[i] );
299			new_atoms[index]->setY( waterY[i] );
300			new_atoms[index]->setZ( waterZ[i] );
301
302			new_atoms[index]->set_vx( 0.0 );
303			new_atoms[index]->set_vy( 0.0 );
304			new_atoms[index]->set_vz( 0.0 );
305
306			dAtom = (DirectionalAtom *) new_atoms[index];
307
308			dAtom->setSUx( 0.0 );
309			dAtom->setSUy( 0.0 );
310			dAtom->setSUz( 1.0 );
311
312	mmeineke	32	dAtom->setA( unitRotMat );
313	mmeineke	16
314			index++;
315			}
316			}
317
318			entry_plug->n_atoms = new_nAtoms;
319			entry_plug->atoms = new_atoms;
320
321			entry_plug->box_x = box_x;
322			entry_plug->box_y = box_y;
323			entry_plug->box_z = box_z;
324
325			DumpWriter* xyz_out = new DumpWriter( entry_plug );
326			xyz_out->writeFinal();
327			delete xyz_out;
328
329			FILE* out_file;
330
331			out_file = fopen( out_name, "w" );
332
333			fprintf(out_file,
334			"#include \"water.mdl\"\n"
335			"#include \"lipid.mdl\"\n"
336			"\n"
337			"nComponents = 2;\n"
338			"component{\n"
339			" type = \"theLipid\";\n"
340			" nMol = %d;\n"
341			"}\n"
342			"\n"
343			"component{\n"
344			" type = \"SSD_water\";\n"
345			" nMol = %d;\n"
346			"}\n"
347			"\n"
348			"initialConfig = \"%s\";\n"
349			"\n"
350			"boxX = %lf;\n"
351			"boxY = %lf;\n"
352			"boxZ = %lf;\n",
353			n_lipids, n_active, entry_plug->finalName,
354			box_x, box_y, box_z );
355
356			fclose( out_file );
357
358			return 0;
359			}