UHBDgrd.hpp

Go to the documentation of this file.

 
#ifndef UHBDGRD_HPP_
#define UHBDGRD_HPP_
 
#include <stdlib.h>
#include <string.h>
#include <sstream>
#include <cstring>
#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
#include "PDBstructure.hpp"
 
template<typename T1>
class UHBD
{
 
public:
    UHBD (char *, char *, bool, float);
    UHBD (char *, float, bool, float);
    UHBD (char *);
    ~UHBD ();
 
    //void loadpdb(char *);
    //void Exclusion();
    //void Extended_Surface;
    void Generate_Skin(float, int);
    void print_slice_grid(int);
    //void Skin_Fill();
 
    std::ifstream file;
    float scale; 
    int grdflag; 
    int one;
    float spacing; 
    char * Title; 
    
    float dummy1;
    float dummy2;
    float dummy3;
    float dummy4;
    float dummy5;
    float dummy6;
    float dummy7;
    int idummy1;
    int idummy2;
    int idummy3;
    int idummy4;
    int idummy5;
    //Dummy variables - end
    
    struct dimensions
    {
        int i; 
        int j; 
        int k; 
    }dim; 
    struct origin
    {
        float x; 
        float y; 
        float z; 
    }o; 
 
    T1 *** Grid; 
    bool *** Skin; 
    //bool *** ExGrid;
    //bool *** ExtGrid;
 
    float probes; 
 
private:
    PDB * structure; 
    int iform; 
    int exclusion_flag; 
    
    void write_uhbd_grid();
    void Test_Format_Grid(char *);
    void Read_Header(char *);
    void Read_Data();
    int excluded_volume(bool ***, double, double, double, double);
    void make_exclusion_grid();
 
    template <typename type>
    type *** Allocate();
 
    template <typename type>
    void Deallocate(type ***);
 
    template <typename type>
    void Initialize(type ***);
 
};

template <typename T1> template <typename type>
type *** UHBD<T1>::Allocate()
{
    type *** grd;
    grd = new type ** [dim.i];
    for (int x=0; x<dim.i; x++){
        grd[x] = new type * [dim.j];
        for (int y=0; y<dim.j; y++)
            grd[x][y] = new type [dim.k];
    }
    return grd;
}


template <typename T1> template <typename type>
void UHBD<T1>::Deallocate(type *** grd)
{
    for (int x=0; x<dim.i; x++){
        for (int y=0; y<dim.j; y++)
            delete[] grd[x][y];
        delete[] grd[x];
   }
}


template <typename T1> template <typename type>
void UHBD<T1>::Initialize(type *** grd)
{
    for (int x=0; x<dim.i; x++){
        for (int y=0; y<dim.j; y++){
            for (int z=0; z<dim.k; z++){
                grd[x][y][z] = 0;
            }
        }
    }
}


template <typename T1>
UHBD<T1>::UHBD(char * pdbfilename, float probe_size, bool pqr_fl, float spcng)
{
    exclusion_flag = 1;
    std::cout << "No UHBD file! \n !!! Assuming an exclusion grid object !!!" << std::endl;
 
    // load pdb file and create a PDB object
    std::cout << "Loading the PDB file ..." << std::endl;
    std::ifstream pdbfile (pdbfilename);
    structure = new PDB(pqr_fl);
    structure->Read_PDB_File(pdbfile);
 
    //Grid, Skin and iform is not initialized here!
    Title = new char [72];
    strcat(Title, "Exclusion grid created by UHBD2dtgrid");
 
    scale = 1.0;
    grdflag = 5;
    one = 0;
    spacing = spcng;
 
    dummy1 = 0.0;
    dummy2 = probe_size; // SDA7 requirement for dummy2 to hold probep!
    dummy3 = 0.0;
    dummy4 = 0.0;
    dummy5 = 0.0;
    dummy6 = 0.0;
    dummy7 = 0.0;
    idummy1 = 0;
    idummy2 = 0;
    idummy3 = 0;
    idummy4 = 0;
    idummy5 = 0;
 
    probes = probe_size;
    make_exclusion_grid(); // call the function directly to create the excluded volume grid.
}
 


template <typename T1>
UHBD<T1>::UHBD(char * pdbfilename, char * Input, bool pqr_fl, float probe_size)
{
    exclusion_flag = 0;
    //Allocate memory for the Title member of the class
    //Always 72 char!!!
    Title = new char [72];
    std::cout << "Opening the UHBD file..." << std::endl;
    //file.open(Input);
    Test_Format_Grid(Input);
    if (iform == 0)
        file.open(Input, std::ifstream::in | std::ifstream::binary);
    else
        file.open(Input, std::ifstream::in);
    Read_Header(Input);
    Grid = Allocate <T1> ();
    Initialize <T1> (Grid);
    Read_Data();
    // Following statement???
    if (file)
        file.close();
    // load pdb file and create a PDB object
    std::cout << "Loading the PDB file ..." << std::endl;
    std::ifstream pdbfile (pdbfilename);
    structure = new PDB(pqr_fl);
    structure->Read_PDB_File(pdbfile);
 
    probes = probe_size;
}


template <typename T1>
UHBD<T1>::UHBD(char * Input)
{
    exclusion_flag = 1;
    //Allocate memory for the Title member of the class
    //Always 72 char!!!
    Title = new char [72];
    std::cout << "Opening the UHBD file..." << std::endl;
    Test_Format_Grid(Input);
    if (iform == 0)
        file.open(Input, std::ifstream::in | std::ifstream::binary);
    else
        file.open(Input, std::ifstream::in);
    Read_Header(Input);
    Grid = Allocate <T1> ();
    Initialize <T1> (Grid);
    Read_Data();
    if (file)
        file.close();
}

template <typename T1>
UHBD<T1>::~UHBD()
{
    delete[] Title;
    Deallocate <bool> (Skin);
    delete[] Skin;
    delete structure;
    if (exclusion_flag == 0){
        Deallocate <T1> (Grid);
    }
    delete[] Grid;
}


template <typename T1>
void UHBD<T1>::print_slice_grid(int k){
    std::string slice_n, out_name;
    slice_n = static_cast<std::ostringstream*>( &(std::ostringstream() << k) )->str();
    out_name = "exclusion" + slice_n + "slice.txt";
 
    std::ofstream excl_slice;
    excl_slice.open (out_name.c_str(), std::ios::out);
 
    for (int i=0; i< dim.i; i++){
        for (int j=0; j< dim.j; j++){
            if (Skin[i][j][k] == true )
                excl_slice << "*";
            else
                excl_slice << " ";
        }
        excl_slice << "\n";
    }
    excl_slice.close();
}


template <typename T1>
int UHBD<T1>::excluded_volume(bool *** IEV, double oex, double oey, double oez, double probep){
    int klim, nb_true_excl;
    int dummy_lim;
    //int nxex, nyex, nzex;
    int c3d_1, c3d_2, c3d_3;
    int i, j, k;
    int ix, jx, kk, jj, ii;
    double rexmax, hdexcl, rexcl, rexcl2;
    double xv1, xv2, xv3, xa1, xa2, xa3;
    double zk, ztemp, yj, ytemp, xi, xtemp, dist2;
    double hexclusion;
    //Atom::Coordinates center;
 
    //center = structure->Center_of_Geometry();
    nb_true_excl = 0;
 
    hexclusion = (double) spacing;
    rexmax = (double) structure->max_vdw + probep;
    hdexcl = hexclusion / 2.0;
    klim = (int) ((rexmax + hexclusion * 1.5) / hexclusion);
 
#pragma omp parallel shared(hdexcl, hexclusion, klim, nb_true_excl, probep, IEV) \
    private(i, j, k, c3d_1, c3d_2, c3d_3, ix, jx, kk, jj, ii, xv1, xv2, xv3, xa1, xa2, xa3, \
            rexcl, rexcl2, zk, ztemp, yj, ytemp, xi, xtemp, dist2, dummy_lim)
{
    #pragma omp for schedule(dynamic) nowait reduction(+:nb_true_excl)
    for (ix = 0; ix < structure->Number_of_Residues; ix++){
        for (jx = 0; jx < structure->Residues[ix].Number_of_Atoms; jx++){
            xv1 = (double) (structure->Residues[ix].Atoms[jx].Coord.X - structure->CoG.X);
            xv2 = (double) (structure->Residues[ix].Atoms[jx].Coord.Y - structure->CoG.Y);
            xv3 = (double) (structure->Residues[ix].Atoms[jx].Coord.Z - structure->CoG.Z);
 
            xa1 = xv1 + hdexcl;
            xa2 = xv2 + hdexcl;
            xa3 = xv3 + hdexcl;
 
            rexcl = (double) structure->Residues[ix].Atoms[jx].radii + probep;
            rexcl2 = rexcl * rexcl;
 
            c3d_1 = (int) ((xa1 - oex) / hexclusion);
            c3d_2 = (int) ((xa2 - oey) / hexclusion);
            c3d_3 = (int) ((xa3 - oez) / hexclusion);
 
            for (kk = -klim; kk <= klim; kk++){
                k = kk + c3d_3;
 
                if (k <= 0 or k > dim.k){
                    //std::cout << "WARNING: k out of range!" << std::endl;
                    continue;
                }
 
                zk = (double) k * hexclusion + oez;
                ztemp = zk - xv3;
 
                if (ztemp > rexcl)
                    continue;
 
                for (jj = -klim; jj <= klim; jj++){
                    j = jj + c3d_2;
 
                    if (j <= 0 or j > dim.j){
                        //std::cout << "WARNING: j out of range!" << std::endl;
                        continue;
                    }
 
                    yj = (double) j * hexclusion + oey;
                    ytemp = yj - xv2;
 
                    if (ytemp > rexcl)
                        continue;
 
                    for (ii = -klim; ii <= klim; ii++){
                        i = ii + c3d_1;
 
                        if (i <= 0 or i > dim.i){
                            //std::cout << "WARNING: i out of range!" << std::endl;
                            continue;
                        }
 
                        if (IEV[i-1][j-1][k-1] == true)
                            continue;
 
                        xi = (double) i * hexclusion + oex;
                        xtemp = xi - xv1;
                        dist2 = xtemp * xtemp + ytemp * ytemp + ztemp * ztemp;
 
                        if (dist2 > rexcl2)
                            continue;
 
#pragma omp critical
                        IEV[i-1][j-1][k-1] = true;
 
                        // When run parallel, the line below would probably not give any reliable numbers.
                        // Not essential for the program though
                        nb_true_excl += 1;
                    }
                }
            }
        }
    }
 
    // If there is any dummy atom in the input pdb file:
    if (structure->Dummy.Number_of_Atoms > 0){
        for (jx = 0; jx < structure->Dummy.Number_of_Atoms; jx++){
            xv1 = (double) (structure->Dummy.Atoms[jx].Coord.X - structure->CoG.X);
            xv2 = (double) (structure->Dummy.Atoms[jx].Coord.Y - structure->CoG.Y);
            xv3 = (double) (structure->Dummy.Atoms[jx].Coord.Z - structure->CoG.Z);
 
            xa1 = xv1 + hdexcl;
            xa2 = xv2 + hdexcl;
            xa3 = xv3 + hdexcl;
 
            rexcl = (double) structure->Dummy.Atoms[jx].radii + probep;
            rexcl2 = rexcl * rexcl;
 
            c3d_1 = (int) ((xa1 - oex) / hexclusion);
            c3d_2 = (int) ((xa2 - oey) / hexclusion);
            c3d_3 = (int) ((xa3 - oez) / hexclusion);
 
            dummy_lim = (int) ((rexcl + hexclusion * 1.5) / hexclusion);
 
            #pragma omp for schedule(dynamic) nowait
            for (kk = -dummy_lim; kk <= dummy_lim; kk++){
                k = kk + c3d_3;
 
                if (k <= 0 or k > dim.k)
                    continue;
 
                zk = (double) k * hexclusion + oez;
                ztemp = zk - xv3;
 
                if (ztemp > rexcl)
                    continue;
 
                for (jj = -dummy_lim; jj <= dummy_lim; jj++){
                    j = jj + c3d_2;
 
                    if (j <= 0 or j > dim.j)
                        continue;
 
                    yj = (double) j * hexclusion + oey;
                    ytemp = yj - xv2;
 
                    if (ytemp > rexcl)
                        continue;
 
                    for (ii = -dummy_lim; ii <= dummy_lim; ii++){
                        i = ii + c3d_1;
 
                        if (i <= 0 or i > dim.i)
                            continue;
 
                        if (IEV[i-1][j-1][k-1] == true)
                            continue;
 
                        xi = (double) i * hexclusion + oex;
                        xtemp = xi - xv1;
                        dist2 = xtemp * xtemp + ytemp * ytemp + ztemp * ztemp;
 
                        if (dist2 > rexcl2)
                            continue;
 
#pragma omp critical
                        IEV[i-1][j-1][k-1] = true;
                        nb_true_excl += 1;
                    }
                }
            }
        }
    }
}
    return nb_true_excl;
}


template <typename T1>
void UHBD<T1>::make_exclusion_grid()
{
    // temporary variables
    int xmin, xmax, ymin, ymax, zmin, zmax;
    int klim, nb_true_excl;
    int nxex, nyex, nzex;
    int c3d_1, c3d_2, c3d_3;
    int i, j, k;
    double rexmax, hdexcl, rexcl, rexcl2;
    double xv1, xv2, xv3, xa1, xa2, xa3;
    double oex, oey, oez;
    double zk, ztemp, yj, ytemp, xi, xtemp, dist2;
    double hexclusion, probep;
    //Atom::Coordinates center;
    xmin = 0;
    ymin = 0;
    zmin = 0;
    xmax = 0;
    ymax = 0;
    zmax = 0;
    nb_true_excl = 0;
 
    hexclusion = (double) spacing;
    probep = (double) probes;
    rexmax = (double) structure->max_vdw + probep;
    hdexcl = hexclusion / 2.0;
    klim = (int) ((rexmax + hexclusion * 1.5) / hexclusion);
    //klim = (int) ceil((rexmax + hexclusion * 1.5d) / hexclusion);
    //++klim;
 
    //std::cout << "klim: " << klim << std::endl;
    std::cout << "Generate exclusion grid with hexclusion: " << hexclusion << " probep: " << probep << std::endl;
 
    //center = structure->Center_of_Geometry();
 
    // find the extents of the input molecule in terms of its coordinates.
    for (int ix = 0; ix < structure->Number_of_Residues; ix++){
        //std::cout << structure->Residues[ix].Number_of_Atoms << " " << ix << std::endl;
        for (int jx = 0; jx < structure->Residues[ix].Number_of_Atoms; jx++){
            xv1 = (double) (structure->Residues[ix].Atoms[jx].Coord.X - structure->CoG.X);
            xv2 = (double) (structure->Residues[ix].Atoms[jx].Coord.Y - structure->CoG.Y);
            xv3 = (double) (structure->Residues[ix].Atoms[jx].Coord.Z - structure->CoG.Z);
            c3d_1 = (int) (xv1  / hexclusion);
            c3d_2 = (int) (xv2  / hexclusion);
            c3d_3 = (int) (xv3  / hexclusion);
            if (c3d_1 < xmin) xmin = c3d_1;
            if (c3d_1 > xmax) xmax = c3d_1;
            if (c3d_2 < ymin) ymin = c3d_2;
            if (c3d_2 > ymax) ymax = c3d_2;
            if (c3d_3 < zmin) zmin = c3d_3;
            if (c3d_3 > zmax) zmax = c3d_3;
        }
    }
    
    // estimate the dimensions of the excluded volume.
    xmin = xmin - klim;
    xmax = xmax + klim;
    ymin = ymin - klim;
    ymax = ymax + klim;
    zmin = zmin - klim;
    zmax = zmax + klim;
    std::cout << "Min & max grid points x: " << xmin << " " << xmax << std::endl;
    std::cout << "Min & max grid points y: " << ymin << " " << ymax << std::endl;
    std::cout << "Min & max grid points z: " << zmin << " " << zmax << std::endl;
 
    nxex = xmax - xmin + 2;
    nyex = ymax - ymin + 2;
    nzex = zmax - zmin + 2;
 
    oex = hexclusion * ((double) xmin - 0.5) - hexclusion;
    oey = hexclusion * ((double) ymin - 0.5) - hexclusion;
    oez = hexclusion * ((double) zmin - 0.5) - hexclusion;
 
    dim.i = nxex;
    dim.j = nyex;
    dim.k = nzex;
 
    // calculate the origin of the grid
    o.x = (float)((double) structure->CoG.X + oex);
    o.y = (float)((double) structure->CoG.Y + oey);
    o.z = (float)((double) structure->CoG.Z + oez);
 
    std::cout << "Exclusion grid information ..." << std::endl;
    std::cout << "Grid dimensions     : " << dim.i << "  " << dim.j << "  " << dim.k << std::endl;
    std::cout << "Grid spacing        : " << spacing << std::endl;
    std::cout << "Grid origin         : " << o.x << "  " << o.y << "  " << o.z << std::endl;
    std::cout << std::endl;
 
    // molecular skin (interaction field) initialize
    Skin = Allocate <bool> ();
    Initialize(Skin);
 
    //construct the excluded volume
    nb_true_excl = excluded_volume(Skin, oex, oey, oez, (double) probes);
 
    std::cout << "The number of points in exclusion grid assigned value 1: " << nb_true_excl << std::endl;
    std::cout << "The number of atoms in the PDB file: " << structure->Number_of_Atoms << std::endl;
    
    /*
    print_slice_grid(45);
    
    std::cout << "Only for testing purposes: Write the excluded volume to a UHBD file " << std::endl;
    write_uhbd_grid();
    */
}


template <typename T1>
void UHBD<T1>::write_uhbd_grid()
{
    FILE* output_file;
    output_file = fopen("test_excluded_volume.bin.grd","wb"); // name of the output grid file
 
    //Write the header of the UHBD file
    int binary_uhbd_type = 4;
    int binary_file_vmd_flag = 160;
    std::cout << "Writing the UHBD (binary) to the output file ..." << std::endl;
    fwrite(&binary_file_vmd_flag, sizeof(int),1,output_file);
    fwrite(Title, sizeof(char), 72, output_file);
    fwrite(&scale, sizeof(float),1,output_file);
    fwrite(&dummy1, sizeof(float),1,output_file);
    fwrite(&binary_uhbd_type, sizeof(int),1,output_file);
    fwrite(&idummy1, sizeof(int),1,output_file);
    fwrite(&idummy2, sizeof(int),1,output_file);
    fwrite(&one, sizeof(int),1,output_file);
    fwrite(&idummy3, sizeof(int),1,output_file);
    fwrite(&dim.i, sizeof(int),1,output_file);
    fwrite(&dim.j, sizeof(int),1,output_file);
    fwrite(&dim.k, sizeof(int),1,output_file);
    fwrite(&spacing, sizeof(float),1,output_file);
    fwrite(&o.x, sizeof(float),1,output_file);
    fwrite(&o.y, sizeof(float),1,output_file);
    fwrite(&o.z, sizeof(float),1,output_file);
    fwrite(&dummy2, sizeof(float),1,output_file);
    fwrite(&dummy3, sizeof(float),1,output_file);
    fwrite(&dummy4, sizeof(float),1,output_file);
    fwrite(&dummy5, sizeof(float),1,output_file);
    fwrite(&dummy6, sizeof(float),1,output_file);
    fwrite(&dummy7, sizeof(float),1,output_file);
    fwrite(&idummy4, sizeof(int),1,output_file);
    fwrite(&idummy5, sizeof(int),1,output_file);
 
    //Write the body of the UHBD file
    float f_value;
    int i_value;
 
    i_value = sizeof(char) * 72 + sizeof(int) * 11 + sizeof(float) * 11;
    fwrite(&i_value, sizeof(int),1,output_file);
 
    for (int k=0; k < dim.k ; k++){
        i_value = sizeof(int) * 3;
        fwrite(&i_value, sizeof(int),1,output_file);
        i_value = k+1;
        fwrite(&i_value, sizeof(int),1,output_file);
        fwrite(&(dim.i), sizeof(int),1,output_file);
        fwrite(&(dim.j), sizeof(int),1,output_file);
        i_value = sizeof(int) * 3;
        fwrite(&i_value, sizeof(int),1,output_file);
        i_value = sizeof(float) * dim.i * dim.j;
        fwrite(&i_value, sizeof(int),1,output_file);
        for (int j=0; j < dim.j ; j++){
            for (int i=0; i < dim.i ; i++){
                f_value = (float) Skin[i][j][k];
                fwrite(&f_value, sizeof(float),1,output_file);
            }
        }
        i_value = sizeof(float) * dim.i * dim.j;
        fwrite(&i_value, sizeof(int),1,output_file);
    }
    fclose(output_file);
}


template <typename T1>
void UHBD<T1>::Test_Format_Grid(char * Input)
{
    std::ifstream file_bt;
    file_bt.open(Input, std::ifstream::in | std::ifstream::binary);
 
    int flag;
    file_bt.read((char *) &flag, sizeof(int));
 
    //Read the first 120 characters
    char * buffer = new char [120];
    file_bt.read(buffer, 120);
 
    char * binary_test;
    binary_test = (char*) memchr (buffer, '\n', 120);
 
    if (flag == 160){
        iform = 0; //Binary file
    }
    else{
        if (binary_test == NULL) { // NO return statement
            binary_test = (char*) memchr (buffer, '\0', 120);
            if (binary_test == NULL) {
                iform = 1; //ASCII file
            }
            else{
                iform = 0; //Binary file
            }
        }
        else{
            iform = 1; //ASCII file
        }
    }
    file_bt.close();
    delete[] buffer;
    //std::cout << "TEST filetype: " << iform << std::endl;
}


template <typename T1>
void UHBD<T1>::Read_Header(char * Input)
{
    // if binary
    if (iform == 0){
        //file.open(Input, std::ifstream::in | std::ifstream::binary);
        if (file) {
            /* to be deleted, check first!
            //get length of file:
            file.seekg (0, file.end);
            int length = file.tellg();
            file.seekg (0, file.beg);
 
            char * buffer = new char [length];
 
            std::cout << "Reading " << length << " characters... ";
            // read data as a block:
            file.read (buffer,length);
 
            if (file)
                std::cout << "UHBD file read successfully." << std::endl;
            else
                std::cout << "error: only " << file.gcount() << " could be read" << std::endl;
            Title = new char [72];
            int buffer_i;
            float buffer_f;
 
            char * buffer_c;
            buffer_c = new char [72];
            */
            
            //binary uhbd reading from position 4.
            //Fortran puts 4 bytes in the beginning.
            file.seekg (4, file.beg);
            file.read(Title, 72);
            file.read((char *) &scale, sizeof(float));
            file.read((char *) &dummy1, sizeof(float));
            file.read((char *) &grdflag, sizeof(int));
            file.read((char *) &idummy1, sizeof(int));
            file.read((char *) &idummy2, sizeof(int));
            file.read((char *) &one, sizeof(int));
            file.read((char *) &idummy3, sizeof(int));
            file.read((char *) &dim.i, sizeof(int));
            file.read((char *) &dim.j, sizeof(int));
            file.read((char *) &dim.k, sizeof(int));
            file.read((char *) &spacing, sizeof(float));
            file.read((char *) &o.x, sizeof(float));
            file.read((char *) &o.y, sizeof(float));
            file.read((char *) &o.z, sizeof(float));
            file.read((char *) &dummy2, sizeof(float));
            file.read((char *) &dummy3, sizeof(float));
            file.read((char *) &dummy4, sizeof(float));
            file.read((char *) &dummy5, sizeof(float));
            file.read((char *) &dummy6, sizeof(float));
            file.read((char *) &dummy7, sizeof(float));
            file.read((char *) &idummy4, sizeof(int));
            file.read((char *) &idummy5, sizeof(int));
 
            //file.close();
            //delete[] buffer;
        }
        else{
            std::cout << "ERROR: Cannot open the UHBD file!" << std::endl;
            exit(0);
        }
    }
    // if ascii
    else{
        std::string title_str, line;
        std::string scale_str, dum1_str, grdflag_str, idum1_str, idum2_str, one_str, idum3_str;
        std::string i_str, j_str, k_str;
        std::string spacing_str, ox_str, oy_str, oz_str;
        std::string dum2_str, dum3_str, dum4_str, dum5_str;
        std::string dum6_str, dum7_str, idum4_str, idum5_str;
        //file.open(Input, std::ifstream::in);
        if (file) {
            getline(file,line);
            title_str = line.substr(0,72);
            getline(file,line);
            scale_str = line.substr(0,12);
            dum1_str = line.substr(12,12);
            grdflag_str = line.substr(24,7);
            idum1_str = line.substr(31,7);
            idum2_str = line.substr(38,7);
            one_str = line.substr(45,7);
            idum3_str = line.substr(52,7);
            getline(file,line);
            i_str = line.substr(0,7);
            j_str = line.substr(7,7);
            k_str = line.substr(14,7);
            spacing_str = line.substr(21,12);
            ox_str = line.substr(33,12);
            oy_str = line.substr(45,12);
            oz_str = line.substr(57,12);
            getline(file,line);
            dum2_str = line.substr(0,12);
            dum3_str = line.substr(12,12);
            dum4_str = line.substr(24,12);
            dum5_str = line.substr(36,12);
            getline(file,line);
            dum6_str = line.substr(0,12);
            dum7_str = line.substr(12,12);
            idum4_str = line.substr(24,7);
            idum5_str = line.substr(31,7);
 
            //Title = title_str.c_str();
            std::strncpy (Title, title_str.c_str(), 72);
 
            std::stringstream (scale_str) >> scale;
            std::stringstream (dum1_str) >> dummy1;
            std::stringstream (grdflag_str) >> grdflag;
            std::stringstream (idum1_str) >> idummy1;
            std::stringstream (idum2_str) >> idummy2;
            std::stringstream (one_str) >> one;
            std::stringstream (idum3_str) >> idummy3;
            std::stringstream (i_str) >> dim.i;
            std::stringstream (j_str) >> dim.j;
            std::stringstream (k_str) >> dim.k;
            std::stringstream (spacing_str) >> spacing;
            std::stringstream (ox_str) >> o.x;
            std::stringstream (oy_str) >> o.y;
            std::stringstream (oz_str) >> o.z;
            std::stringstream (dum2_str) >> dummy2;
            std::stringstream (dum3_str) >> dummy3;
            std::stringstream (dum4_str) >> dummy4;
            std::stringstream (dum5_str) >> dummy5;
            std::stringstream (dum6_str) >> dummy6;
            std::stringstream (dum7_str) >> dummy7;
            std::stringstream (idum4_str) >> idummy4;
            std::stringstream (idum5_str) >> idummy5;
        }
        else{
            std::cout << "ERROR: Cannot open the UHBD file" << std::endl;
            exit(0);
        }
    }
}


template <typename T1>
void UHBD<T1>::Read_Data()
{
    std::cout << "Reading the UHBD grid ..." << std::endl;
    std::cout << "Grid dimensions     : " << dim.i << "  " << dim.j << "  " << dim.k << std::endl;
    std::cout << "Grid spacing        : " << spacing << std::endl;
    std::cout << "Grid origin         : " << o.x << "  " << o.y << "  " << o.z << std::endl;
    std::cout << std::endl;
 
    // if file is binary
    if (iform == 0){
        if (file) {
            int kx, im, jm, ufo1, ufo2, ufo3, ufo4 ;
            float data_buffer;
            for (int k = 0; k<dim.k; k++){
 
                file.read((char *) &ufo1, sizeof(float));
                file.read((char *) &ufo2, sizeof(float));
                file.read((char *) &kx, sizeof(float));
                file.read((char *) &im, sizeof(float));
                file.read((char *) &jm, sizeof(float));
                file.read((char *) &ufo3, sizeof(float));
                file.read((char *) &ufo4, sizeof(float));
                //std::cout << (float) kx << " " << (float) im << " " <<  (float) jm << " " << (float) ufo1 << " " << (float) ufo2 << " " << (float) ufo3 << " " << (float) ufo4 << std::endl;
                //std::cout << kx << " " << im << " " << jm << " " << ufo1 << " " << ufo2 << " " << ufo3 << " " << ufo4 << std::endl;
                for (int j = 0; j<dim.j; j++){
                    for (int i = 0; i<dim.i; i++){
                        file.read((char *) &data_buffer, sizeof(float));
                        //std::cout << data_buffer << "\t" << (int) data_buffer << std::endl;
                        Grid[i][j][k] = data_buffer;
                        //std::cout << i << " " << j << " " << k << " " << Grid[i][j][k] << std::endl;
                    }
                }
            }
            //file.close();
        }
        else{
            std::cout << "ERROR: Cannot open the UHBD file" << std::endl;
            exit(0);
        }
    }
    else{
        if (file) {
            std::string line;
            //int ij_length = dim.i * dim.j;
            //T1 value;
 
            getline(file,line);
            for (int k = 0; k<dim.k; k++){
                for (int j = 0; j<dim.j; j++){
                    for (int i = 0; i<dim.i; i++){
                        file >> Grid[i][j][k];
                        //std::cout << i << " " << j << " " << k << " " << Grid[i][j][k] << std::endl;
                    }
                }
                getline(file,line);
                //std::cout << line << std::endl;
                getline(file,line);
                //std::cout << line << std::endl;
            }
 
            //file.close();
        }
        else{
            std::cout << "ERROR: Cannot open the UHBD file" << std::endl;
            exit(0);
        }
    }
}
 
/*template <typename T1>
void UHBD<T1>::loadpdb(char *pdbfilename)
{
    std::cout << "Loading the PDB file ..." << std::endl;
    std::ifstream pdbfile (pdbfilename);
    structure = new PDB;
    structure->Read_PDB_File(pdbfile,MAX);
}*/


template <typename T1>
void UHBD<T1>::Generate_Skin(float skin, int iterate_inside_skin)
{
    std::cout << "Generating the molecule skin points ..." << std::endl;
 
    double probsk;
    double oex, oey, oez;
    int ih, jh, kh, KK, JJ, II;
    int nb_true_excl;
    //Atom::Coordinates center;
 
    //center = structure->Center_of_Geometry();
    oex = (double) o.x - structure->CoG.X;
    oey = (double) o.y - structure->CoG.Y;
    oez = (double) o.z - structure->CoG.Z;
 
    bool *** IEX;
    IEX = Allocate<bool>();
    Initialize<bool>(IEX);
 
    probsk = (double) probes;
    nb_true_excl = excluded_volume(IEX, oex, oey, oez, probsk);
 
    bool *** ISK;
    ISK = Allocate<bool>();
    Initialize<bool>(ISK);
 
    probsk = probes + skin;
    nb_true_excl = excluded_volume(ISK, oex, oey, oez, probsk);
 
    //SKIN FILL PART
    Skin = Allocate <bool> ();
    Initialize(Skin);
    for (int i = 0; i < dim.i; i ++){
        for (int j = 0; j < dim.j; j++){
            for (int k = 0; k < dim.k; k++){
                if (IEX[i][j][k]== true){
                    Skin[i][j][k]=false;
                }
                else
                    if (ISK[i][j][k]== true)
                        Skin[i][j][k]=true;
                    else
                        Skin[i][j][k]=false;
            }
        }
    }
 
    Deallocate<bool>(ISK);
    delete[] ISK;
 
    // This part is added to include several layers of data points inside the protein
    // number of layers is given by iterate_inside_skin variable
    bool *** IIN;
    IIN = Allocate<bool>();
    Initialize<bool>(IIN);
 
    if (iterate_inside_skin != 0){
        for (int i_skin = 0; i_skin < iterate_inside_skin; i_skin++){
            for (int i = 0; i < dim.i; i ++){
                for (int j = 0; j < dim.j; j++){
                    for (int k = 0; k < dim.k; k++){
                        if (Skin[i][j][k] == true){
                            if ((i+1 < dim.i) && (IEX[i+1][j][k] == true))
                                IIN[i+1][j][k] = 1;
                            if ((i-1 >= 0) && (IEX[i-1][j][k] == true))
                                IIN[i-1][j][k] = 1;
                            if ((j+1 < dim.j) && (IEX[i][j+1][k] == true))
                                IIN[i][j+1][k] = 1;
                            if ((j-1 >= 0) && (IEX[i][j-1][k] == true))
                                IIN[i][j-1][k] = 1;
                            if ((k+1 < dim.k) && (IEX[i][j][k+1] == true))
                                IIN[i][j][k+1] = 1;
                            if ((k-1 >= 0) && (IEX[i][j][k-1] == true))
                                IIN[i][j][k-1] = 1;
                        }
                        else
                            continue;
                    }
                }
            }
            for (int i = 0; i < dim.i; i ++){
                for (int j = 0; j < dim.j; j++){
                    for (int k = 0; k < dim.k; k++){
                        if (IIN[i][j][k] == true){
                            Skin[i][j][k] = true;
                            IEX[i][j][k] = false;
                        }
                    }
                }
            }
        }
    }
    //
    Deallocate<bool>(IIN);
    Deallocate<bool>(IEX);
    delete[] IIN;
    delete[] IEX;
}
 
#endif /* UHBDGRD_HPP_ */