SDA (SDA flex)  7.2
Simulation of Diffusional Association
Loading...
Searching...
No Matches
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
23template<typename T1>
24class UHBD
25{
26
27public:
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;
60 //Dummy variables - end
61
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
82private:
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);
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
106template <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
123template <typename T1> template <typename type>
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
137template <typename T1> template <typename type>
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
159template <typename T1>
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
213template <typename T1>
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
250template <typename T1>
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
272template <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
292template <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
326template <typename T1>
327int 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.0;
346 klim = (int) ((rexmax + hexclusion * 1.5) / 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.5) / 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
509template <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.0;
536 klim = (int) ((rexmax + hexclusion * 1.5) / 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.5) - hexclusion;
580 oey = hexclusion * ((double) ymin - 0.5) - hexclusion;
581 oez = hexclusion * ((double) zmin - 0.5) - 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
623template <typename T1>
625{
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 }
686}
687
689
692template <typename T1>
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
736template <typename T1>
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;
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
878template <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>
949void 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
970template <typename T1>
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;
989
990 probsk = (double) probes;
991 nb_true_excl = excluded_volume(IEX, oex, oey, oez, probsk);
992
993 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
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>();
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 //
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 DTGrid3D.hpp:65
int idummy3
Definition DTGrid3D.hpp:148
int grdflag
Definition DTGrid3D.hpp:144
float dummy5
Definition DTGrid3D.hpp:157
float dummy7
Definition DTGrid3D.hpp:159
float z
Definition DTGrid3D.hpp:122
float dummy3
Definition DTGrid3D.hpp:155
float dummy6
Definition DTGrid3D.hpp:158
char * Title
Definition DTGrid3D.hpp:143
float x
Definition DTGrid3D.hpp:120
int iform
grid file format: binary = 0; ascii = 1
Definition DTGrid3D.hpp:86
struct DTGrid3D::@2::origin o
float y
Definition DTGrid3D.hpp:121
float scale
scaling factor
Definition DTGrid3D.hpp:92
int j
Definition DTGrid3D.hpp:114
float dummy1
Definition DTGrid3D.hpp:153
float spacing
Definition DTGrid3D.hpp:93
int i
Definition DTGrid3D.hpp:113
int idummy2
Definition DTGrid3D.hpp:147
int idummy4
Definition DTGrid3D.hpp:149
float dummy4
Definition DTGrid3D.hpp:156
int k
Definition DTGrid3D.hpp:115
int idummy5
Definition DTGrid3D.hpp:150
struct DTGrid3D::@2::dimensions dim
int idummy1
Definition DTGrid3D.hpp:146
int one
Definition DTGrid3D.hpp:145
float dummy2
Definition DTGrid3D.hpp:154
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
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