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 = 5;
  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];
  
  rotMat[0][0] = 1.0;
  rotMat[0][1] = 0.0;
  rotMat[0][2] = 0.0;
  
  rotMat[1][0] = 0.0;
  rotMat[1][1] = 1.0;
  rotMat[1][2] = 0.0;
  
  rotMat[2][0] = 0.0;
  rotMat[2][1] = 0.0;
  rotMat[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() );
      
        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() );
      
      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;
}
Revision:	17
Committed:	Tue Jul 9 18:45:04 2002 UTC (21 years, 11 months ago) by mmeineke
File size:	8232 byte(s)
Log Message:	This commit was generated by cvs2svn to compensate for changes in r16, which included commits to RCS files with non-trunk default branches.
#	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			int n_lipidsY = 5;
38			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
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;
88
89			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			}
115
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
134			double min_x, min_y, min_z;
135			double max_x, max_y, max_z;
136			double test_x, test_y, test_z;
137
138			max_x = min_x = group_atoms[0]->getX();
139			max_y = min_y = group_atoms[0]->getY();
140			max_z = min_z = group_atoms[0]->getZ();
141
142			for(i=0; i<group_nAtoms; i++){
143
144			test_x = group_atoms[i]->getX();
145			test_y = group_atoms[i]->getY();
146			test_z = group_atoms[i]->getZ();
147
148			if( test_x < min_x ) min_x = test_x;
149			if( test_y < min_y ) min_y = test_y;
150			if( test_z < min_z ) min_z = test_z;
151
152			if( test_x > max_x ) max_x = test_x;
153			if( test_y > max_y ) max_y = test_y;
154			if( test_z > max_z ) max_z = test_z;
155			}
156
157			double box_x = max_x - min_x + 2*water_shell;
158			double box_y = max_y - min_y + 2*water_shell;
159			double box_z = max_z - min_z + 2*water_shell;
160
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 );
164
165			box_x = water_cell * n_cellX;
166			box_y = water_cell * n_cellY;
167			box_z = water_cell * n_cellZ;
168
169			int n_water = n_cellX * n_cellY * n_cellZ * 4;
170
171			double *waterX = new double[n_water];
172			double *waterY = new double[n_water];
173			double *waterZ = new double[n_water];
174
175			double cx, cy, cz;
176
177			cx = 0.0;
178			cy = 0.0;
179			cz = 0.0;
180			for(i=0; i<group_nAtoms; i++){
181			cx += group_atoms[i]->getX();
182			cy += group_atoms[i]->getY();
183			cz += group_atoms[i]->getZ();
184			}
185			cx /= group_nAtoms;
186			cy /= group_nAtoms;
187			cz /= group_nAtoms;
188
189			double x0 = cx - ( box_x * 0.5 );
190			double y0 = cy - ( box_y * 0.5 );
191			double z0 = cz - ( box_z * 0.5 );
192
193			index = 0;
194			for( i=0; i < n_cellX; i++ ){
195			for( j=0; j < n_cellY; j++ ){
196			for( k=0; k < n_cellZ; k++ ){
197
198			waterX[index] = i * water_cell + x0;
199			waterY[index] = j * water_cell + y0;
200			waterZ[index] = k * water_cell + z0;
201			index++;
202
203			waterX[index] = i * water_cell + 0.5 * water_cell + x0;
204			waterY[index] = j * water_cell + 0.5 * water_cell + y0;
205			waterZ[index] = k * water_cell + z0;
206			index++;
207
208			waterX[index] = i * water_cell + x0;
209			waterY[index] = j * water_cell + 0.5 * water_cell + y0;
210			waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
211			index++;
212
213			waterX[index] = i * water_cell + 0.5 * water_cell + x0;
214			waterY[index] = j * water_cell + y0;
215			waterZ[index] = k * water_cell + 0.5 * water_cell + z0;
216			index++;
217			}
218			}
219			}
220
221			int *isActive = new int[n_water];
222			for(i=0; i<n_water; i++) isActive[i] = 1;
223
224			int n_active = n_water;
225			double dx, dy, dz;
226			double dx2, dy2, dz2, dSqr;
227			double rCutSqr = water_padding * water_padding;
228
229			for(i=0; ( (i<n_water) && isActive[i] ); i++){
230			for(j=0; ( (j<group_nAtoms) && isActive[i] ); j++){
231
232			dx = waterX[i] - group_atoms[j]->getX();
233			dy = waterY[i] - group_atoms[j]->getY();
234			dz = waterZ[i] - group_atoms[j]->getZ();
235
236			dx2 = dx * dx;
237			dy2 = dy * dy;
238			dz2 = dz * dz;
239
240			dSqr = dx2 + dy2 + dz2;
241			if( dSqr < rCutSqr ){
242			isActive[i] = 0;
243			n_active--;
244			}
245			}
246			}
247
248			std::cerr << "final n_waters = " << n_active << "\n";
249
250			int new_nAtoms = group_nAtoms + n_active;
251			Atom** new_atoms = new Atom*[new_nAtoms];
252
253			index = 0;
254			for(i=0; i<group_nAtoms; i++ ){
255
256			if( group_atoms[i]->isDirectional() ){
257			dAtom = (DirectionalAtom *)group_atoms[i];
258
259			dAtomNew = new DirectionalAtom();
260			dAtomNew->setSUx( dAtom->getSUx() );
261			dAtomNew->setSUx( dAtom->getSUx() );
262			dAtomNew->setSUx( dAtom->getSUx() );
263
264			dAtomNew->setA( rotMat );
265
266			new_atoms[index] = dAtomNew;
267			}
268			else{
269
270			new_atoms[index] = new GeneralAtom();
271			}
272
273			new_atoms[index]->setType( group_atoms[i]->getType() );
274
275			new_atoms[index]->setX( group_atoms[i]->getX() );
276			new_atoms[index]->setY( group_atoms[i]->getY() );
277			new_atoms[index]->setZ( group_atoms[i]->getZ() );
278
279			new_atoms[index]->set_vx( 0.0 );
280			new_atoms[index]->set_vy( 0.0 );
281			new_atoms[index]->set_vz( 0.0 );
282
283			index++;
284			}
285
286
287
288
289			for(i=0; i<n_water; i++){
290			if(isActive[i]){
291
292			new_atoms[index] = new DirectionalAtom();
293			new_atoms[index]->setType( "SSD" );
294
295			new_atoms[index]->setX( waterX[i] );
296			new_atoms[index]->setY( waterY[i] );
297			new_atoms[index]->setZ( waterZ[i] );
298
299			new_atoms[index]->set_vx( 0.0 );
300			new_atoms[index]->set_vy( 0.0 );
301			new_atoms[index]->set_vz( 0.0 );
302
303			dAtom = (DirectionalAtom *) new_atoms[index];
304
305			dAtom->setSUx( 0.0 );
306			dAtom->setSUy( 0.0 );
307			dAtom->setSUz( 1.0 );
308
309			dAtom->setA( rotMat );
310
311			index++;
312			}
313			}
314
315			entry_plug->n_atoms = new_nAtoms;
316			entry_plug->atoms = new_atoms;
317
318			entry_plug->box_x = box_x;
319			entry_plug->box_y = box_y;
320			entry_plug->box_z = box_z;
321
322			DumpWriter* xyz_out = new DumpWriter( entry_plug );
323			xyz_out->writeFinal();
324			delete xyz_out;
325
326			FILE* out_file;
327
328			out_file = fopen( out_name, "w" );
329
330			fprintf(out_file,
331			"#include \"water.mdl\"\n"
332			"#include \"lipid.mdl\"\n"
333			"\n"
334			"nComponents = 2;\n"
335			"component{\n"
336			" type = \"theLipid\";\n"
337			" nMol = %d;\n"
338			"}\n"
339			"\n"
340			"component{\n"
341			" type = \"SSD_water\";\n"
342			" nMol = %d;\n"
343			"}\n"
344			"\n"
345			"initialConfig = \"%s\";\n"
346			"\n"
347			"boxX = %lf;\n"
348			"boxY = %lf;\n"
349			"boxZ = %lf;\n",
350			n_lipids, n_active, entry_plug->finalName,
351			box_x, box_y, box_z );
352
353			fclose( out_file );
354
355			return 0;
356			}