SDA (SDA flex)  7.2
Simulation of Diffusional Association
UHBDgrd.hpp
Go to the documentation of this file.
1 
10 #ifndef UHBDGRD_HPP_
11 #define UHBDGRD_HPP_
12 
13 #include <stdlib.h>
14 #include <string.h>
15 #include <sstream>
16 #include <cstring>
17 #include <iostream>
18 #include <fstream>
19 #include <sstream>
20 #include <math.h>
21 #include "PDBstructure.hpp"
22 
23 template<typename T1>
24 class UHBD
25 {
26 
27 public:
28  UHBD (char *, char *, bool, float);
29  UHBD (char *, float, bool, float);
30  UHBD (char *);
31  ~UHBD ();
32 
33  //void loadpdb(char *);
34  //void Exclusion();
35  //void Extended_Surface;
36  void Generate_Skin(float, int);
37  void print_slice_grid(int);
38  //void Skin_Fill();
39 
40  std::ifstream file;
41  float scale;
42  int grdflag;
43  int one;
44  float spacing;
45  char * Title;
46 
48  float dummy1;
49  float dummy2;
50  float dummy3;
51  float dummy4;
52  float dummy5;
53  float dummy6;
54  float dummy7;
55  int idummy1;
56  int idummy2;
57  int idummy3;
58  int idummy4;
59  int idummy5;
60  //Dummy variables - end
61 
62  struct dimensions
63  {
64  int i;
65  int j;
66  int k;
67  }dim;
68  struct origin
69  {
70  float x;
71  float y;
72  float z;
73  }o;
74 
75  T1 *** Grid;
76  bool *** Skin;
77  //bool *** ExGrid;
78  //bool *** ExtGrid;
79 
80  float probes;
81 
82 private:
84  int iform;
86 
87  void write_uhbd_grid();
88  void Test_Format_Grid(char *);
89  void Read_Header(char *);
90  void Read_Data();
91  int excluded_volume(bool ***, double, double, double, double);
92  void make_exclusion_grid();
93 
94  template <typename type>
95  type *** Allocate();
96 
97  template <typename type>
98  void Deallocate(type ***);
99 
100  template <typename type>
101  void Initialize(type ***);
102 
103 };
104 
106 template <typename T1> template <typename type>
108 {
109  type *** grd;
110  grd = new type ** [dim.i];
111  for (int x=0; x<dim.i; x++){
112  grd[x] = new type * [dim.j];
113  for (int y=0; y<dim.j; y++)
114  grd[x][y] = new type [dim.k];
115  }
116  return grd;
117 }
118 
120 
123 template <typename T1> template <typename type>
124 void UHBD<T1>::Deallocate(type *** grd)
125 {
126  for (int x=0; x<dim.i; x++){
127  for (int y=0; y<dim.j; y++)
128  delete[] grd[x][y];
129  delete[] grd[x];
130  }
131 }
132 
134 
137 template <typename T1> template <typename type>
138 void UHBD<T1>::Initialize(type *** grd)
139 {
140  for (int x=0; x<dim.i; x++){
141  for (int y=0; y<dim.j; y++){
142  for (int z=0; z<dim.k; z++){
143  grd[x][y][z] = 0;
144  }
145  }
146  }
147 }
148 
150 
159 template <typename T1>
160 UHBD<T1>::UHBD(char * pdbfilename, float probe_size, bool pqr_fl, float spcng)
161 {
162  exclusion_flag = 1;
163  std::cout << "No UHBD file! \n !!! Assuming an exclusion grid object !!!" << std::endl;
164 
165  // load pdb file and create a PDB object
166  std::cout << "Loading the PDB file ..." << std::endl;
167  std::ifstream pdbfile (pdbfilename);
168  structure = new PDB(pqr_fl);
169  structure->Read_PDB_File(pdbfile);
170 
171  //Grid, Skin and iform is not initialized here!
172  Title = new char [72];
173  strcat(Title, "Exclusion grid created by UHBD2dtgrid");
174 
175  scale = 1.0;
176  grdflag = 5;
177  one = 0;
178  spacing = spcng;
179 
180  dummy1 = 0.0;
181  dummy2 = probe_size; // SDA7 requirement for dummy2 to hold probep!
182  dummy3 = 0.0;
183  dummy4 = 0.0;
184  dummy5 = 0.0;
185  dummy6 = 0.0;
186  dummy7 = 0.0;
187  idummy1 = 0;
188  idummy2 = 0;
189  idummy3 = 0;
190  idummy4 = 0;
191  idummy5 = 0;
192 
193  probes = probe_size;
194  make_exclusion_grid(); // call the function directly to create the excluded volume grid.
195 }
196 
197 
199 
213 template <typename T1>
214 UHBD<T1>::UHBD(char * pdbfilename, char * Input, bool pqr_fl, float probe_size)
215 {
216  exclusion_flag = 0;
217  //Allocate memory for the Title member of the class
218  //Always 72 char!!!
219  Title = new char [72];
220  std::cout << "Opening the UHBD file..." << std::endl;
221  //file.open(Input);
222  Test_Format_Grid(Input);
223  if (iform == 0)
224  file.open(Input, std::ifstream::in | std::ifstream::binary);
225  else
226  file.open(Input, std::ifstream::in);
227  Read_Header(Input);
228  Grid = Allocate <T1> ();
229  Initialize <T1> (Grid);
230  Read_Data();
231  // Following statement???
232  if (file)
233  file.close();
234  // load pdb file and create a PDB object
235  std::cout << "Loading the PDB file ..." << std::endl;
236  std::ifstream pdbfile (pdbfilename);
237  structure = new PDB(pqr_fl);
238  structure->Read_PDB_File(pdbfile);
239 
240  probes = probe_size;
241 }
242 
244 
250 template <typename T1>
251 UHBD<T1>::UHBD(char * Input)
252 {
253  exclusion_flag = 1;
254  //Allocate memory for the Title member of the class
255  //Always 72 char!!!
256  Title = new char [72];
257  std::cout << "Opening the UHBD file..." << std::endl;
258  Test_Format_Grid(Input);
259  if (iform == 0)
260  file.open(Input, std::ifstream::in | std::ifstream::binary);
261  else
262  file.open(Input, std::ifstream::in);
263  Read_Header(Input);
264  Grid = Allocate <T1> ();
265  Initialize <T1> (Grid);
266  Read_Data();
267  if (file)
268  file.close();
269 }
270 
272 template <typename T1>
274 {
275  delete[] Title;
276  Deallocate <bool> (Skin);
277  delete[] Skin;
278  delete structure;
279  if (exclusion_flag == 0){
280  Deallocate <T1> (Grid);
281  }
282  delete[] Grid;
283 }
284 
286 
292 template <typename T1>
294  std::string slice_n, out_name;
295  slice_n = static_cast<std::ostringstream*>( &(std::ostringstream() << k) )->str();
296  out_name = "exclusion" + slice_n + "slice.txt";
297 
298  std::ofstream excl_slice;
299  excl_slice.open (out_name.c_str(), std::ios::out);
300 
301  for (int i=0; i< dim.i; i++){
302  for (int j=0; j< dim.j; j++){
303  if (Skin[i][j][k] == true )
304  excl_slice << "*";
305  else
306  excl_slice << " ";
307  }
308  excl_slice << "\n";
309  }
310  excl_slice.close();
311 }
312 
314 
326 template <typename T1>
327 int UHBD<T1>::excluded_volume(bool *** IEV, double oex, double oey, double oez, double probep){
328  int klim, nb_true_excl;
329  int dummy_lim;
330  //int nxex, nyex, nzex;
331  int c3d_1, c3d_2, c3d_3;
332  int i, j, k;
333  int ix, jx, kk, jj, ii;
334  double rexmax, hdexcl, rexcl, rexcl2;
335  double xv1, xv2, xv3, xa1, xa2, xa3;
336  double zk, ztemp, yj, ytemp, xi, xtemp, dist2;
337  double hexclusion;
338  //Atom::Coordinates center;
339 
340  //center = structure->Center_of_Geometry();
341  nb_true_excl = 0;
342 
343  hexclusion = (double) spacing;
344  rexmax = (double) structure->max_vdw + probep;
345  hdexcl = hexclusion / 2.0d;
346  klim = (int) ((rexmax + hexclusion * 1.5d) / hexclusion);
347 
348 #pragma omp parallel shared(hdexcl, hexclusion, klim, nb_true_excl, probep, IEV) \
349  private(i, j, k, c3d_1, c3d_2, c3d_3, ix, jx, kk, jj, ii, xv1, xv2, xv3, xa1, xa2, xa3, \
350  rexcl, rexcl2, zk, ztemp, yj, ytemp, xi, xtemp, dist2, dummy_lim)
351 {
352  #pragma omp for schedule(dynamic) nowait reduction(+:nb_true_excl)
353  for (ix = 0; ix < structure->Number_of_Residues; ix++){
354  for (jx = 0; jx < structure->Residues[ix].Number_of_Atoms; jx++){
355  xv1 = (double) (structure->Residues[ix].Atoms[jx].Coord.X - structure->CoG.X);
356  xv2 = (double) (structure->Residues[ix].Atoms[jx].Coord.Y - structure->CoG.Y);
357  xv3 = (double) (structure->Residues[ix].Atoms[jx].Coord.Z - structure->CoG.Z);
358 
359  xa1 = xv1 + hdexcl;
360  xa2 = xv2 + hdexcl;
361  xa3 = xv3 + hdexcl;
362 
363  rexcl = (double) structure->Residues[ix].Atoms[jx].radii + probep;
364  rexcl2 = rexcl * rexcl;
365 
366  c3d_1 = (int) ((xa1 - oex) / hexclusion);
367  c3d_2 = (int) ((xa2 - oey) / hexclusion);
368  c3d_3 = (int) ((xa3 - oez) / hexclusion);
369 
370  for (kk = -klim; kk <= klim; kk++){
371  k = kk + c3d_3;
372 
373  if (k <= 0 or k > dim.k){
374  //std::cout << "WARNING: k out of range!" << std::endl;
375  continue;
376  }
377 
378  zk = (double) k * hexclusion + oez;
379  ztemp = zk - xv3;
380 
381  if (ztemp > rexcl)
382  continue;
383 
384  for (jj = -klim; jj <= klim; jj++){
385  j = jj + c3d_2;
386 
387  if (j <= 0 or j > dim.j){
388  //std::cout << "WARNING: j out of range!" << std::endl;
389  continue;
390  }
391 
392  yj = (double) j * hexclusion + oey;
393  ytemp = yj - xv2;
394 
395  if (ytemp > rexcl)
396  continue;
397 
398  for (ii = -klim; ii <= klim; ii++){
399  i = ii + c3d_1;
400 
401  if (i <= 0 or i > dim.i){
402  //std::cout << "WARNING: i out of range!" << std::endl;
403  continue;
404  }
405 
406  if (IEV[i-1][j-1][k-1] == true)
407  continue;
408 
409  xi = (double) i * hexclusion + oex;
410  xtemp = xi - xv1;
411  dist2 = xtemp * xtemp + ytemp * ytemp + ztemp * ztemp;
412 
413  if (dist2 > rexcl2)
414  continue;
415 
416 #pragma omp critical
417  IEV[i-1][j-1][k-1] = true;
418 
419  // When run parallel, the line below would probably not give any reliable numbers.
420  // Not essential for the program though
421  nb_true_excl += 1;
422  }
423  }
424  }
425  }
426  }
427 
428  // If there is any dummy atom in the input pdb file:
429  if (structure->Dummy.Number_of_Atoms > 0){
430  for (jx = 0; jx < structure->Dummy.Number_of_Atoms; jx++){
431  xv1 = (double) (structure->Dummy.Atoms[jx].Coord.X - structure->CoG.X);
432  xv2 = (double) (structure->Dummy.Atoms[jx].Coord.Y - structure->CoG.Y);
433  xv3 = (double) (structure->Dummy.Atoms[jx].Coord.Z - structure->CoG.Z);
434 
435  xa1 = xv1 + hdexcl;
436  xa2 = xv2 + hdexcl;
437  xa3 = xv3 + hdexcl;
438 
439  rexcl = (double) structure->Dummy.Atoms[jx].radii + probep;
440  rexcl2 = rexcl * rexcl;
441 
442  c3d_1 = (int) ((xa1 - oex) / hexclusion);
443  c3d_2 = (int) ((xa2 - oey) / hexclusion);
444  c3d_3 = (int) ((xa3 - oez) / hexclusion);
445 
446  dummy_lim = (int) ((rexcl + hexclusion * 1.5d) / hexclusion);
447 
448  #pragma omp for schedule(dynamic) nowait
449  for (kk = -dummy_lim; kk <= dummy_lim; kk++){
450  k = kk + c3d_3;
451 
452  if (k <= 0 or k > dim.k)
453  continue;
454 
455  zk = (double) k * hexclusion + oez;
456  ztemp = zk - xv3;
457 
458  if (ztemp > rexcl)
459  continue;
460 
461  for (jj = -dummy_lim; jj <= dummy_lim; jj++){
462  j = jj + c3d_2;
463 
464  if (j <= 0 or j > dim.j)
465  continue;
466 
467  yj = (double) j * hexclusion + oey;
468  ytemp = yj - xv2;
469 
470  if (ytemp > rexcl)
471  continue;
472 
473  for (ii = -dummy_lim; ii <= dummy_lim; ii++){
474  i = ii + c3d_1;
475 
476  if (i <= 0 or i > dim.i)
477  continue;
478 
479  if (IEV[i-1][j-1][k-1] == true)
480  continue;
481 
482  xi = (double) i * hexclusion + oex;
483  xtemp = xi - xv1;
484  dist2 = xtemp * xtemp + ytemp * ytemp + ztemp * ztemp;
485 
486  if (dist2 > rexcl2)
487  continue;
488 
489 #pragma omp critical
490  IEV[i-1][j-1][k-1] = true;
491  nb_true_excl += 1;
492  }
493  }
494  }
495  }
496  }
497 }
498  return nb_true_excl;
499 }
500 
502 
509 template <typename T1>
511 {
512  // temporary variables
513  int xmin, xmax, ymin, ymax, zmin, zmax;
514  int klim, nb_true_excl;
515  int nxex, nyex, nzex;
516  int c3d_1, c3d_2, c3d_3;
517  int i, j, k;
518  double rexmax, hdexcl, rexcl, rexcl2;
519  double xv1, xv2, xv3, xa1, xa2, xa3;
520  double oex, oey, oez;
521  double zk, ztemp, yj, ytemp, xi, xtemp, dist2;
522  double hexclusion, probep;
523  //Atom::Coordinates center;
524  xmin = 0;
525  ymin = 0;
526  zmin = 0;
527  xmax = 0;
528  ymax = 0;
529  zmax = 0;
530  nb_true_excl = 0;
531 
532  hexclusion = (double) spacing;
533  probep = (double) probes;
534  rexmax = (double) structure->max_vdw + probep;
535  hdexcl = hexclusion / 2.0d;
536  klim = (int) ((rexmax + hexclusion * 1.5d) / hexclusion);
537  //klim = (int) ceil((rexmax + hexclusion * 1.5d) / hexclusion);
538  //++klim;
539 
540  //std::cout << "klim: " << klim << std::endl;
541  std::cout << "Generate exclusion grid with hexclusion: " << hexclusion << " probep: " << probep << std::endl;
542 
543  //center = structure->Center_of_Geometry();
544 
545  // find the extents of the input molecule in terms of its coordinates.
546  for (int ix = 0; ix < structure->Number_of_Residues; ix++){
547  //std::cout << structure->Residues[ix].Number_of_Atoms << " " << ix << std::endl;
548  for (int jx = 0; jx < structure->Residues[ix].Number_of_Atoms; jx++){
549  xv1 = (double) (structure->Residues[ix].Atoms[jx].Coord.X - structure->CoG.X);
550  xv2 = (double) (structure->Residues[ix].Atoms[jx].Coord.Y - structure->CoG.Y);
551  xv3 = (double) (structure->Residues[ix].Atoms[jx].Coord.Z - structure->CoG.Z);
552  c3d_1 = (int) (xv1 / hexclusion);
553  c3d_2 = (int) (xv2 / hexclusion);
554  c3d_3 = (int) (xv3 / hexclusion);
555  if (c3d_1 < xmin) xmin = c3d_1;
556  if (c3d_1 > xmax) xmax = c3d_1;
557  if (c3d_2 < ymin) ymin = c3d_2;
558  if (c3d_2 > ymax) ymax = c3d_2;
559  if (c3d_3 < zmin) zmin = c3d_3;
560  if (c3d_3 > zmax) zmax = c3d_3;
561  }
562  }
563 
564  // estimate the dimensions of the excluded volume.
565  xmin = xmin - klim;
566  xmax = xmax + klim;
567  ymin = ymin - klim;
568  ymax = ymax + klim;
569  zmin = zmin - klim;
570  zmax = zmax + klim;
571  std::cout << "Min & max grid points x: " << xmin << " " << xmax << std::endl;
572  std::cout << "Min & max grid points y: " << ymin << " " << ymax << std::endl;
573  std::cout << "Min & max grid points z: " << zmin << " " << zmax << std::endl;
574 
575  nxex = xmax - xmin + 2;
576  nyex = ymax - ymin + 2;
577  nzex = zmax - zmin + 2;
578 
579  oex = hexclusion * ((double) xmin - 0.5d) - hexclusion;
580  oey = hexclusion * ((double) ymin - 0.5d) - hexclusion;
581  oez = hexclusion * ((double) zmin - 0.5d) - hexclusion;
582 
583  dim.i = nxex;
584  dim.j = nyex;
585  dim.k = nzex;
586 
587  // calculate the origin of the grid
588  o.x = (float)((double) structure->CoG.X + oex);
589  o.y = (float)((double) structure->CoG.Y + oey);
590  o.z = (float)((double) structure->CoG.Z + oez);
591 
592  std::cout << "Exclusion grid information ..." << std::endl;
593  std::cout << "Grid dimensions : " << dim.i << " " << dim.j << " " << dim.k << std::endl;
594  std::cout << "Grid spacing : " << spacing << std::endl;
595  std::cout << "Grid origin : " << o.x << " " << o.y << " " << o.z << std::endl;
596  std::cout << std::endl;
597 
598  // molecular skin (interaction field) initialize
599  Skin = Allocate <bool> ();
600  Initialize(Skin);
601 
602  //construct the excluded volume
603  nb_true_excl = excluded_volume(Skin, oex, oey, oez, (double) probes);
604 
605  std::cout << "The number of points in exclusion grid assigned value 1: " << nb_true_excl << std::endl;
606  std::cout << "The number of atoms in the PDB file: " << structure->Number_of_Atoms << std::endl;
607 
608  /*
609  print_slice_grid(45);
610 
611  std::cout << "Only for testing purposes: Write the excluded volume to a UHBD file " << std::endl;
612  write_uhbd_grid();
613  */
614 }
616 
623 template <typename T1>
625 {
626  FILE* output_file;
627  output_file = fopen("test_excluded_volume.bin.grd","wb"); // name of the output grid file
628 
629  //Write the header of the UHBD file
630  int binary_uhbd_type = 4;
631  int binary_file_vmd_flag = 160;
632  std::cout << "Writing the UHBD (binary) to the output file ..." << std::endl;
633  fwrite(&binary_file_vmd_flag, sizeof(int),1,output_file);
634  fwrite(Title, sizeof(char), 72, output_file);
635  fwrite(&scale, sizeof(float),1,output_file);
636  fwrite(&dummy1, sizeof(float),1,output_file);
637  fwrite(&binary_uhbd_type, sizeof(int),1,output_file);
638  fwrite(&idummy1, sizeof(int),1,output_file);
639  fwrite(&idummy2, sizeof(int),1,output_file);
640  fwrite(&one, sizeof(int),1,output_file);
641  fwrite(&idummy3, sizeof(int),1,output_file);
642  fwrite(&dim.i, sizeof(int),1,output_file);
643  fwrite(&dim.j, sizeof(int),1,output_file);
644  fwrite(&dim.k, sizeof(int),1,output_file);
645  fwrite(&spacing, sizeof(float),1,output_file);
646  fwrite(&o.x, sizeof(float),1,output_file);
647  fwrite(&o.y, sizeof(float),1,output_file);
648  fwrite(&o.z, sizeof(float),1,output_file);
649  fwrite(&dummy2, sizeof(float),1,output_file);
650  fwrite(&dummy3, sizeof(float),1,output_file);
651  fwrite(&dummy4, sizeof(float),1,output_file);
652  fwrite(&dummy5, sizeof(float),1,output_file);
653  fwrite(&dummy6, sizeof(float),1,output_file);
654  fwrite(&dummy7, sizeof(float),1,output_file);
655  fwrite(&idummy4, sizeof(int),1,output_file);
656  fwrite(&idummy5, sizeof(int),1,output_file);
657 
658  //Write the body of the UHBD file
659  float f_value;
660  int i_value;
661 
662  i_value = sizeof(char) * 72 + sizeof(int) * 11 + sizeof(float) * 11;
663  fwrite(&i_value, sizeof(int),1,output_file);
664 
665  for (int k=0; k < dim.k ; k++){
666  i_value = sizeof(int) * 3;
667  fwrite(&i_value, sizeof(int),1,output_file);
668  i_value = k+1;
669  fwrite(&i_value, sizeof(int),1,output_file);
670  fwrite(&(dim.i), sizeof(int),1,output_file);
671  fwrite(&(dim.j), sizeof(int),1,output_file);
672  i_value = sizeof(int) * 3;
673  fwrite(&i_value, sizeof(int),1,output_file);
674  i_value = sizeof(float) * dim.i * dim.j;
675  fwrite(&i_value, sizeof(int),1,output_file);
676  for (int j=0; j < dim.j ; j++){
677  for (int i=0; i < dim.i ; i++){
678  f_value = (float) Skin[i][j][k];
679  fwrite(&f_value, sizeof(float),1,output_file);
680  }
681  }
682  i_value = sizeof(float) * dim.i * dim.j;
683  fwrite(&i_value, sizeof(int),1,output_file);
684  }
685  fclose(output_file);
686 }
687 
689 
692 template <typename T1>
693 void UHBD<T1>::Test_Format_Grid(char * Input)
694 {
695  std::ifstream file_bt;
696  file_bt.open(Input, std::ifstream::in | std::ifstream::binary);
697 
698  int flag;
699  file_bt.read((char *) &flag, sizeof(int));
700 
701  //Read the first 120 characters
702  char * buffer = new char [120];
703  file_bt.read(buffer, 120);
704 
705  char * binary_test;
706  binary_test = (char*) memchr (buffer, '\n', 120);
707 
708  if (flag == 160){
709  iform = 0; //Binary file
710  }
711  else{
712  if (binary_test == NULL) { // NO return statement
713  binary_test = (char*) memchr (buffer, '\0', 120);
714  if (binary_test == NULL) {
715  iform = 1; //ASCII file
716  }
717  else{
718  iform = 0; //Binary file
719  }
720  }
721  else{
722  iform = 1; //ASCII file
723  }
724  }
725  file_bt.close();
726  delete[] buffer;
727  //std::cout << "TEST filetype: " << iform << std::endl;
728 }
729 
731 
736 template <typename T1>
737 void UHBD<T1>::Read_Header(char * Input)
738 {
739  // if binary
740  if (iform == 0){
741  //file.open(Input, std::ifstream::in | std::ifstream::binary);
742  if (file) {
743  /* to be deleted, check first!
744  //get length of file:
745  file.seekg (0, file.end);
746  int length = file.tellg();
747  file.seekg (0, file.beg);
748 
749  char * buffer = new char [length];
750 
751  std::cout << "Reading " << length << " characters... ";
752  // read data as a block:
753  file.read (buffer,length);
754 
755  if (file)
756  std::cout << "UHBD file read successfully." << std::endl;
757  else
758  std::cout << "error: only " << file.gcount() << " could be read" << std::endl;
759  Title = new char [72];
760  int buffer_i;
761  float buffer_f;
762 
763  char * buffer_c;
764  buffer_c = new char [72];
765  */
766 
767  //binary uhbd reading from position 4.
768  //Fortran puts 4 bytes in the beginning.
769  file.seekg (4, file.beg);
770  file.read(Title, 72);
771  file.read((char *) &scale, sizeof(float));
772  file.read((char *) &dummy1, sizeof(float));
773  file.read((char *) &grdflag, sizeof(int));
774  file.read((char *) &idummy1, sizeof(int));
775  file.read((char *) &idummy2, sizeof(int));
776  file.read((char *) &one, sizeof(int));
777  file.read((char *) &idummy3, sizeof(int));
778  file.read((char *) &dim.i, sizeof(int));
779  file.read((char *) &dim.j, sizeof(int));
780  file.read((char *) &dim.k, sizeof(int));
781  file.read((char *) &spacing, sizeof(float));
782  file.read((char *) &o.x, sizeof(float));
783  file.read((char *) &o.y, sizeof(float));
784  file.read((char *) &o.z, sizeof(float));
785  file.read((char *) &dummy2, sizeof(float));
786  file.read((char *) &dummy3, sizeof(float));
787  file.read((char *) &dummy4, sizeof(float));
788  file.read((char *) &dummy5, sizeof(float));
789  file.read((char *) &dummy6, sizeof(float));
790  file.read((char *) &dummy7, sizeof(float));
791  file.read((char *) &idummy4, sizeof(int));
792  file.read((char *) &idummy5, sizeof(int));
793 
794  //file.close();
795  //delete[] buffer;
796  }
797  else{
798  std::cout << "ERROR: Cannot open the UHBD file!" << std::endl;
799  exit(0);
800  }
801  }
802  // if ascii
803  else{
804  std::string title_str, line;
805  std::string scale_str, dum1_str, grdflag_str, idum1_str, idum2_str, one_str, idum3_str;
806  std::string i_str, j_str, k_str;
807  std::string spacing_str, ox_str, oy_str, oz_str;
808  std::string dum2_str, dum3_str, dum4_str, dum5_str;
809  std::string dum6_str, dum7_str, idum4_str, idum5_str;
810  //file.open(Input, std::ifstream::in);
811  if (file) {
812  getline(file,line);
813  title_str = line.substr(0,72);
814  getline(file,line);
815  scale_str = line.substr(0,12);
816  dum1_str = line.substr(12,12);
817  grdflag_str = line.substr(24,7);
818  idum1_str = line.substr(31,7);
819  idum2_str = line.substr(38,7);
820  one_str = line.substr(45,7);
821  idum3_str = line.substr(52,7);
822  getline(file,line);
823  i_str = line.substr(0,7);
824  j_str = line.substr(7,7);
825  k_str = line.substr(14,7);
826  spacing_str = line.substr(21,12);
827  ox_str = line.substr(33,12);
828  oy_str = line.substr(45,12);
829  oz_str = line.substr(57,12);
830  getline(file,line);
831  dum2_str = line.substr(0,12);
832  dum3_str = line.substr(12,12);
833  dum4_str = line.substr(24,12);
834  dum5_str = line.substr(36,12);
835  getline(file,line);
836  dum6_str = line.substr(0,12);
837  dum7_str = line.substr(12,12);
838  idum4_str = line.substr(24,7);
839  idum5_str = line.substr(31,7);
840 
841  //Title = title_str.c_str();
842  std::strncpy (Title, title_str.c_str(), 72);
843 
844  std::stringstream (scale_str) >> scale;
845  std::stringstream (dum1_str) >> dummy1;
846  std::stringstream (grdflag_str) >> grdflag;
847  std::stringstream (idum1_str) >> idummy1;
848  std::stringstream (idum2_str) >> idummy2;
849  std::stringstream (one_str) >> one;
850  std::stringstream (idum3_str) >> idummy3;
851  std::stringstream (i_str) >> dim.i;
852  std::stringstream (j_str) >> dim.j;
853  std::stringstream (k_str) >> dim.k;
854  std::stringstream (spacing_str) >> spacing;
855  std::stringstream (ox_str) >> o.x;
856  std::stringstream (oy_str) >> o.y;
857  std::stringstream (oz_str) >> o.z;
858  std::stringstream (dum2_str) >> dummy2;
859  std::stringstream (dum3_str) >> dummy3;
860  std::stringstream (dum4_str) >> dummy4;
861  std::stringstream (dum5_str) >> dummy5;
862  std::stringstream (dum6_str) >> dummy6;
863  std::stringstream (dum7_str) >> dummy7;
864  std::stringstream (idum4_str) >> idummy4;
865  std::stringstream (idum5_str) >> idummy5;
866  }
867  else{
868  std::cout << "ERROR: Cannot open the UHBD file" << std::endl;
869  exit(0);
870  }
871  }
872 }
873 
875 
878 template <typename T1>
880 {
881  std::cout << "Reading the UHBD grid ..." << std::endl;
882  std::cout << "Grid dimensions : " << dim.i << " " << dim.j << " " << dim.k << std::endl;
883  std::cout << "Grid spacing : " << spacing << std::endl;
884  std::cout << "Grid origin : " << o.x << " " << o.y << " " << o.z << std::endl;
885  std::cout << std::endl;
886 
887  // if file is binary
888  if (iform == 0){
889  if (file) {
890  int kx, im, jm, ufo1, ufo2, ufo3, ufo4 ;
891  float data_buffer;
892  for (int k = 0; k<dim.k; k++){
893 
894  file.read((char *) &ufo1, sizeof(float));
895  file.read((char *) &ufo2, sizeof(float));
896  file.read((char *) &kx, sizeof(float));
897  file.read((char *) &im, sizeof(float));
898  file.read((char *) &jm, sizeof(float));
899  file.read((char *) &ufo3, sizeof(float));
900  file.read((char *) &ufo4, sizeof(float));
901  //std::cout << (float) kx << " " << (float) im << " " << (float) jm << " " << (float) ufo1 << " " << (float) ufo2 << " " << (float) ufo3 << " " << (float) ufo4 << std::endl;
902  //std::cout << kx << " " << im << " " << jm << " " << ufo1 << " " << ufo2 << " " << ufo3 << " " << ufo4 << std::endl;
903  for (int j = 0; j<dim.j; j++){
904  for (int i = 0; i<dim.i; i++){
905  file.read((char *) &data_buffer, sizeof(float));
906  //std::cout << data_buffer << "\t" << (int) data_buffer << std::endl;
907  Grid[i][j][k] = data_buffer;
908  //std::cout << i << " " << j << " " << k << " " << Grid[i][j][k] << std::endl;
909  }
910  }
911  }
912  //file.close();
913  }
914  else{
915  std::cout << "ERROR: Cannot open the UHBD file" << std::endl;
916  exit(0);
917  }
918  }
919  else{
920  if (file) {
921  std::string line;
922  //int ij_length = dim.i * dim.j;
923  //T1 value;
924 
925  getline(file,line);
926  for (int k = 0; k<dim.k; k++){
927  for (int j = 0; j<dim.j; j++){
928  for (int i = 0; i<dim.i; i++){
929  file >> Grid[i][j][k];
930  //std::cout << i << " " << j << " " << k << " " << Grid[i][j][k] << std::endl;
931  }
932  }
933  getline(file,line);
934  //std::cout << line << std::endl;
935  getline(file,line);
936  //std::cout << line << std::endl;
937  }
938 
939  //file.close();
940  }
941  else{
942  std::cout << "ERROR: Cannot open the UHBD file" << std::endl;
943  exit(0);
944  }
945  }
946 }
947 
948 /*template <typename T1>
949 void UHBD<T1>::loadpdb(char *pdbfilename)
950 {
951  std::cout << "Loading the PDB file ..." << std::endl;
952  std::ifstream pdbfile (pdbfilename);
953  structure = new PDB;
954  structure->Read_PDB_File(pdbfile,MAX);
955 }*/
956 
958 
970 template <typename T1>
971 void UHBD<T1>::Generate_Skin(float skin, int iterate_inside_skin)
972 {
973  std::cout << "Generating the molecule skin points ..." << std::endl;
974 
975  double probsk;
976  double oex, oey, oez;
977  int ih, jh, kh, KK, JJ, II;
978  int nb_true_excl;
979  //Atom::Coordinates center;
980 
981  //center = structure->Center_of_Geometry();
982  oex = (double) o.x - structure->CoG.X;
983  oey = (double) o.y - structure->CoG.Y;
984  oez = (double) o.z - structure->CoG.Z;
985 
986  bool *** IEX;
987  IEX = Allocate<bool>();
988  Initialize<bool>(IEX);
989 
990  probsk = (double) probes;
991  nb_true_excl = excluded_volume(IEX, oex, oey, oez, probsk);
992 
993  bool *** ISK;
994  ISK = Allocate<bool>();
995  Initialize<bool>(ISK);
996 
997  probsk = probes + skin;
998  nb_true_excl = excluded_volume(ISK, oex, oey, oez, probsk);
999 
1000  //SKIN FILL PART
1001  Skin = Allocate <bool> ();
1002  Initialize(Skin);
1003  for (int i = 0; i < dim.i; i ++){
1004  for (int j = 0; j < dim.j; j++){
1005  for (int k = 0; k < dim.k; k++){
1006  if (IEX[i][j][k]== true){
1007  Skin[i][j][k]=false;
1008  }
1009  else
1010  if (ISK[i][j][k]== true)
1011  Skin[i][j][k]=true;
1012  else
1013  Skin[i][j][k]=false;
1014  }
1015  }
1016  }
1017 
1018  Deallocate<bool>(ISK);
1019  delete[] ISK;
1020 
1021  // This part is added to include several layers of data points inside the protein
1022  // number of layers is given by iterate_inside_skin variable
1023  bool *** IIN;
1024  IIN = Allocate<bool>();
1025  Initialize<bool>(IIN);
1026 
1027  if (iterate_inside_skin != 0){
1028  for (int i_skin = 0; i_skin < iterate_inside_skin; i_skin++){
1029  for (int i = 0; i < dim.i; i ++){
1030  for (int j = 0; j < dim.j; j++){
1031  for (int k = 0; k < dim.k; k++){
1032  if (Skin[i][j][k] == true){
1033  if ((i+1 < dim.i) && (IEX[i+1][j][k] == true))
1034  IIN[i+1][j][k] = 1;
1035  if ((i-1 >= 0) && (IEX[i-1][j][k] == true))
1036  IIN[i-1][j][k] = 1;
1037  if ((j+1 < dim.j) && (IEX[i][j+1][k] == true))
1038  IIN[i][j+1][k] = 1;
1039  if ((j-1 >= 0) && (IEX[i][j-1][k] == true))
1040  IIN[i][j-1][k] = 1;
1041  if ((k+1 < dim.k) && (IEX[i][j][k+1] == true))
1042  IIN[i][j][k+1] = 1;
1043  if ((k-1 >= 0) && (IEX[i][j][k-1] == true))
1044  IIN[i][j][k-1] = 1;
1045  }
1046  else
1047  continue;
1048  }
1049  }
1050  }
1051  for (int i = 0; i < dim.i; i ++){
1052  for (int j = 0; j < dim.j; j++){
1053  for (int k = 0; k < dim.k; k++){
1054  if (IIN[i][j][k] == true){
1055  Skin[i][j][k] = true;
1056  IEX[i][j][k] = false;
1057  }
1058  }
1059  }
1060  }
1061  }
1062  }
1063  //
1064  Deallocate<bool>(IIN);
1065  Deallocate<bool>(IEX);
1066  delete[] IIN;
1067  delete[] IEX;
1068 }
1069 
1070 #endif /* UHBDGRD_HPP_ */
Group functions to read-write PDB files and describe the molecular structures inside (Header)
Definition: PDBstructure.hpp:78
Definition: UHBDgrd.hpp:25
float spacing
Definition: UHBDgrd.hpp:44
float dummy6
Definition: UHBDgrd.hpp:53
void Test_Format_Grid(char *)
Determine the type of grid: binary or ascii?
Definition: UHBDgrd.hpp:693
int idummy5
Definition: UHBDgrd.hpp:59
void print_slice_grid(int)
Print a 2D slice of an excluded volume grid.
Definition: UHBDgrd.hpp:293
float dummy3
Definition: UHBDgrd.hpp:50
void Generate_Skin(float, int)
Calculate the molecule skin (molecular interaction field)
Definition: UHBDgrd.hpp:971
char * Title
grid spacing/resolution
Definition: UHBDgrd.hpp:45
type *** Allocate()
create a 3d array and allocate memory for it.
Definition: UHBDgrd.hpp:107
float probes
Boolean for Skin! expanded for other types, though.
Definition: UHBDgrd.hpp:80
int exclusion_flag
Binary : 0, ASCII : 1.
Definition: UHBDgrd.hpp:85
UHBD(char *, char *, bool, float)
UHBD constructor - 2.
Definition: UHBDgrd.hpp:214
int idummy4
Definition: UHBDgrd.hpp:58
int excluded_volume(bool ***, double, double, double, double)
Calculate the excluded volume.
Definition: UHBDgrd.hpp:327
float dummy5
Definition: UHBDgrd.hpp:52
struct UHBD::dimensions dim
void Deallocate(type ***)
deallocate the memory for a 3D array
Definition: UHBDgrd.hpp:124
float dummy2
Definition: UHBDgrd.hpp:49
int grdflag
scaling for every value
Definition: UHBDgrd.hpp:42
bool *** Skin
Store the grid data in a 3d array.
Definition: UHBDgrd.hpp:76
~UHBD()
UHBD destructor.
Definition: UHBDgrd.hpp:273
void make_exclusion_grid()
construct the excluded volume grid
Definition: UHBDgrd.hpp:510
void Read_Header(char *)
Read the header of an input UHBD file.
Definition: UHBDgrd.hpp:737
float scale
Definition: UHBDgrd.hpp:41
float dummy1
Text header kept in the UHBD file.
Definition: UHBDgrd.hpp:48
float dummy4
Definition: UHBDgrd.hpp:51
int iform
structure object for the input molecule
Definition: UHBDgrd.hpp:84
void write_uhbd_grid()
Exclusion : 1, Non-exclusion : 0.
Definition: UHBDgrd.hpp:624
void Initialize(type ***)
initialize values of a given 3D array to zero
Definition: UHBDgrd.hpp:138
int idummy1
Definition: UHBDgrd.hpp:55
int one
flag for the type of grid: UHBD, DT-Grid, etc.
Definition: UHBDgrd.hpp:43
std::ifstream file
Definition: UHBDgrd.hpp:40
float dummy7
Definition: UHBDgrd.hpp:54
T1 *** Grid
grid origin coordinates
Definition: UHBDgrd.hpp:75
int idummy3
Definition: UHBDgrd.hpp:57
struct UHBD::origin o
void Read_Data()
Read the body of the input UHBD file.
Definition: UHBDgrd.hpp:879
PDB * structure
probe size
Definition: UHBDgrd.hpp:83
int idummy2
Definition: UHBDgrd.hpp:56
subroutine make_exclusion_grid(sgrid, tmp_atom_pos, tmp_nat, center, atom_vdw, max_vdw, probep, hexclusion)
Use intermediate module mod_exclusion_grid.
Definition: make_exclusion_grid.f90:54
real(kind=8), parameter k
Definition: mod_gouy_chapman.f90:34
Definition: UHBDgrd.hpp:63
int k
number of grid points in y direction
Definition: UHBDgrd.hpp:66
int i
Definition: UHBDgrd.hpp:64
int j
number of grid points in x direction
Definition: UHBDgrd.hpp:65
dimensions of the grid (number of grid points)
Definition: UHBDgrd.hpp:69
float x
Definition: UHBDgrd.hpp:70
float z
origin y coordinate
Definition: UHBDgrd.hpp:72
float y
origin x coordinate
Definition: UHBDgrd.hpp:71
Imprint/Privacy