DTgrd.hpp

Go to the documentation of this file.

#ifndef DTGRD_HPP_
#define DTGRD_HPP_

#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string.h>
#include <sstream>
#include <iomanip>
#include <math.h>
#include <sys/stat.h>
#include <vector>
#include "UHBDgrd.hpp"

template<typename T1>
class DT_FILE
{
public:
    DT_FILE (int);
    DT_FILE (int, bool);
    DT_FILE (int, int, int, int);
    ~DT_FILE ();

private:
    int grid_type_flag; 
    bool pqr_flag; 

    int LengthX; 
    int LengthY; 
    int LengthZ; 

    UHBD<T1> * uhbd; 
    int uhbd_alloc;

    std::vector<T1> Value;   
    std::vector<int> xcoord; 
    std::vector<int> ycoord; 
    std::vector<int> zcoord; 
    std::vector<int> acc;    
    std::vector<int> acc2D;  
    std::vector<int> acc3D;  
    std::vector<int> proj1D; 
    std::vector<int> proj2D; 

    int check_neighbours(bool ***, bool ***, int, int, int);
    bool *** exterior_grid(bool ***);
    int write_dtgrdfilebody_ascii(FILE *);
    int write_dtgrdfilebody_binary(FILE *);

public:
    void UHBD_upload(char *, char *, float, float);
    void exclusion(char *, float, float);
    void exclusion(char *);
    int Generate_DTGrid(bool ***, T1 ***);

    int print_layer(bool ***, int, int);
    int write_dtgrdfile_ascii(FILE *);
    int write_dtgrdfile_binary(FILE *);
    int write_dtgrdfile_ascii_nouhbd(FILE *);
    int write_dtgrdfile_binary_nouhbd(FILE *);
};


template <typename T1>
DT_FILE<T1>::DT_FILE(int flag){
    grid_type_flag = flag;
    uhbd_alloc = 0;
    pqr_flag = false;
}


template <typename T1>
DT_FILE<T1>::DT_FILE(int flag, bool pqr_fl){
    grid_type_flag = flag;
    uhbd_alloc = 0;
    pqr_flag = pqr_fl;
}


template <typename T1>
DT_FILE<T1>::DT_FILE(int flag, int i, int j, int k){
    grid_type_flag = flag;
    LengthX = i;
    LengthY = j;
    LengthZ = k;
    uhbd_alloc = 0;
}


template <typename T1>
void DT_FILE<T1>::UHBD_upload(char * PDBFILE, char * UHBDFILE, float probe_size, float skin_size){
    uhbd = new UHBD<T1>(PDBFILE, UHBDFILE, pqr_flag, probe_size);
    // 5 given below is the number of layers inside the vdw surface to be added to the molecular skin
    uhbd->Generate_Skin(skin_size, 5);
    LengthX = uhbd->dim.i;
    LengthY = uhbd->dim.j;
    LengthZ = uhbd->dim.k;
    uhbd_alloc = 1;

    Generate_DTGrid(uhbd->Skin,uhbd->Grid);
}


template <typename T1>
void DT_FILE<T1>::exclusion(char * PDBFILE, float probe_size, float spacing){
    uhbd = new UHBD<T1>(PDBFILE, probe_size, pqr_flag, spacing);

    LengthX = uhbd->dim.i;
    LengthY = uhbd->dim.j;
    LengthZ = uhbd->dim.k;
    uhbd_alloc = 1;

    Generate_DTGrid(uhbd->Skin,uhbd->Grid);
}


template <typename T1>
void DT_FILE<T1>::exclusion(char * UHBDFILE){
    uhbd = new UHBD<T1>(UHBDFILE);

    LengthX = uhbd->dim.i;
    LengthY = uhbd->dim.j;
    LengthZ = uhbd->dim.k;
    uhbd_alloc = 1;

    Generate_DTGrid(uhbd->Grid,uhbd->Grid);
}


template <typename T1>
int DT_FILE<T1>::check_neighbours(bool *** Exterior, bool *** Skin_Grid, int x, int y, int z){

    if (x >= LengthX or y >= LengthY or z >= LengthZ )
        return 0;

    if (x < 0 or y < 0 or z < 0 )
        return 0;

    if (Skin_Grid[x][y][z] == false){
        if (Exterior[x][y][z] == false){
            Exterior[x][y][z] = true;
            check_neighbours(Exterior, Skin_Grid, x+1, y, z);
            check_neighbours(Exterior, Skin_Grid, x-1, y, z);
            check_neighbours(Exterior, Skin_Grid, x, y+1, z);
            check_neighbours(Exterior, Skin_Grid, x, y-1, z);
            check_neighbours(Exterior, Skin_Grid, x, y, z+1);
            check_neighbours(Exterior, Skin_Grid, x, y, z-1);
        }
    }
    return 0;
}


template <typename T1>
int DT_FILE<T1>::print_layer(bool *** Exterior, int dimension, int layer){

    if (dimension == 0){
        for (int i=0;i<LengthX;i++){
            for (int j=0;j<LengthY;j++){
                if (Exterior[i][j][layer] == true)
                    std::cout << "X";
                else
                    std::cout << " ";
            }
            std::cout << std::endl;
        }
    }
    else if (dimension == 1){
        for (int j=0;j<LengthY;j++){
            for (int k=0;k<LengthZ;k++){
                if (Exterior[layer][j][k] == true)
                    std::cout << "X";
                else
                    std::cout << " ";
            }
            std::cout << std::endl;
        }
    }
    else if (dimension == 2){
        for (int i=0;i<LengthX;i++){
            for (int k=0;k<LengthZ;k++){
                if (Exterior[i][layer][k] == true)
                    std::cout << "X";
                else
                    std::cout << " ";
            }
            std::cout << std::endl;
        }
    }
    else {
        std::cout << "Check the dimension flag!!" << std::endl;
    }
    return 0;
}


template <typename T1>
bool *** DT_FILE<T1>::exterior_grid(bool *** Skin_Grid){
    bool *** Exterior;
    Exterior = new bool ** [LengthX];
    for (int i=0;i<LengthX;i++){
        Exterior[i] = new bool * [LengthY];
        for (int j=0;j<LengthY;j++){
            Exterior[i][j] = new bool [LengthZ];
        }
    }

    for (int i=0;i<LengthX;i++){
        for (int j=0;j<LengthY;j++){
            for (int k=0;k<LengthZ;k++){
                Exterior[i][j][k] = false;
            }
        }
    }

    for (int i=0;i<LengthX;i++){
        for (int j=0;j<LengthY;j++){
            if (Skin_Grid[i][j][0] == false){
                if (Exterior[i][j][0] == false){
                    Exterior[i][j][0] = true;
                    check_neighbours(Exterior, Skin_Grid, i+1, j, 0);
                    check_neighbours(Exterior, Skin_Grid, i-1, j, 0);
                    check_neighbours(Exterior, Skin_Grid, i, j+1, 0);
                    check_neighbours(Exterior, Skin_Grid, i, j-1, 0);
                    check_neighbours(Exterior, Skin_Grid, i, j,   1);
                }
            }
            if (Skin_Grid[i][j][LengthZ-1] == false){
                if (Exterior[i][j][LengthZ-1] == false){
                    Exterior[i][j][LengthZ-1] = true;
                    check_neighbours(Exterior, Skin_Grid, i+1, j, LengthZ-1);
                    check_neighbours(Exterior, Skin_Grid, i-1, j, LengthZ-1);
                    check_neighbours(Exterior, Skin_Grid, i, j+1, LengthZ-1);
                    check_neighbours(Exterior, Skin_Grid, i, j-1, LengthZ-1);
                    check_neighbours(Exterior, Skin_Grid, i, j,   LengthZ-2);
                }
            }
        }
    }
    for (int j=0;j<LengthY;j++){
        for (int k=0;k<LengthZ;k++){
            if (Skin_Grid[0][j][k] == false){
                if (Exterior[0][j][k] == false){
                    Exterior[0][j][k] = true;
                    check_neighbours(Exterior, Skin_Grid, 0, j, k+1);
                    check_neighbours(Exterior, Skin_Grid, 0, j, k-1);
                    check_neighbours(Exterior, Skin_Grid, 0, j+1, k);
                    check_neighbours(Exterior, Skin_Grid, 0, j-1, k);
                    check_neighbours(Exterior, Skin_Grid, 1, j,   k);
                }
            }
            if (Skin_Grid[LengthX-1][j][k] == false){
                if (Exterior[LengthX-1][j][k] == false){
                    Exterior[LengthX-1][j][k] = true;
                    check_neighbours(Exterior, Skin_Grid, LengthX-1, j, k+1);
                    check_neighbours(Exterior, Skin_Grid, LengthX-1, j, k-1);
                    check_neighbours(Exterior, Skin_Grid, LengthX-1, j+1, k);
                    check_neighbours(Exterior, Skin_Grid, LengthX-1, j-1, k);
                    check_neighbours(Exterior, Skin_Grid, LengthX-2, j,   k);
                }
            }
        }
    }

    for (int i=0;i<LengthX;i++){
        for (int k=0;k<LengthZ;k++){
            if (Skin_Grid[i][0][k] == false){
                if (Exterior[i][0][k] == false){
                    Exterior[i][0][k] = true;
                    check_neighbours(Exterior, Skin_Grid, i, 0, k+1);
                    check_neighbours(Exterior, Skin_Grid, i, 0, k-1);
                    check_neighbours(Exterior, Skin_Grid, i+1, 0, k);
                    check_neighbours(Exterior, Skin_Grid, i-1, 0, k);
                    check_neighbours(Exterior, Skin_Grid, i, 1, k);
                }
            }
            if (Skin_Grid[i][LengthY-1][k] == false){
                if (Exterior[i][LengthY-1][k] == false){
                    Exterior[i][LengthY-1][k] = true;
                    check_neighbours(Exterior, Skin_Grid, i, LengthY-1, k+1);
                    check_neighbours(Exterior, Skin_Grid, i, LengthY-1, k-1);
                    check_neighbours(Exterior, Skin_Grid, i+1, LengthY-1, k);
                    check_neighbours(Exterior, Skin_Grid, i-1, LengthY-1, k);
                    check_neighbours(Exterior, Skin_Grid, i, LengthY-2,   k);
                }
            }
        }
    }
    return Exterior;
}

template <typename T1>
int DT_FILE<T1>::Generate_DTGrid(bool *** Skin_Grid, T1 *** UHBD_Grid){

    std::cout << "Generating the DT-Grid structure ..." << std::endl;
    int componentY, componentZ;
    int prev_componentY, prev_componentZ;

    prev_componentY = -1;

    for (int i=0;i<LengthX;i++){
        prev_componentZ = -1;
        componentY = 0;
        for (int j=0; j<LengthY; j++){
            componentZ = 0;
            for (int k=0; k<LengthZ; k++){
                if (grid_type_flag == 0){
                    if (Skin_Grid[i][j][k] == true){
                        if (k==0){
                            componentZ += 1;
                            zcoord.push_back(k);
                            acc.push_back(0);
                            proj1D.push_back(0);
                            proj1D.push_back(((int) zcoord.size())-1);
                        }
                        else{
                            if (Skin_Grid[i][j][k-1] == false){
                                componentZ += 1;
                                zcoord.push_back(k);
                                acc.push_back(0);
                                if (componentZ == 1){
                                    proj1D.push_back(0);
                                    proj1D.push_back(((int) zcoord.size())-1);
                                }
                            }
                        }
                    }
                }
                else{
                    if (Skin_Grid[i][j][k] == true){
                        Value.push_back(UHBD_Grid[i][j][k]);

                        if (k==0){
                            componentZ += 1;
                            zcoord.push_back(k);
                            acc.push_back(((int) Value.size())-1);
                            proj1D.push_back(((int) Value.size())-1);
                            proj1D.push_back(((int) zcoord.size())-1);
                        }
                        else{
                            if (Skin_Grid[i][j][k-1] == false){
                                componentZ += 1;
                                zcoord.push_back(k);
                                acc.push_back(((int) Value.size())-1);
                                if (componentZ == 1){
                                    proj1D.push_back(((int) Value.size())-1);
                                    proj1D.push_back(((int) zcoord.size())-1);
                                }
                            }
                        }
                    }
                }

                if (k == (LengthZ-1)){
                    if (Skin_Grid[i][j][k] == true){
                        zcoord.push_back(k);
                    }
                    else{
                        if (Skin_Grid[i][j][k-1] == true){
                            zcoord.push_back((k-1));
                        }
                    }
                }
                else{
                    if (k != 0){
                        if (Skin_Grid[i][j][k-1] == true){
                            if (Skin_Grid[i][j][k] == false){
                                zcoord.push_back((k-1));
                            }
                        }
                    }
                }
            }

            if (j==0){
                if (componentZ > 0){
                    ycoord.push_back(j);
                    acc2D.push_back(((int) proj1D.size())-2);
                    componentY += 1;
                    proj2D.push_back(((int) proj1D.size())-2);
                    proj2D.push_back(((int) ycoord.size())-1);

                }
            }
            else{
                if (prev_componentZ == 0){
                    if (componentZ > 0){
                        ycoord.push_back(j);
                        acc2D.push_back(((int) proj1D.size())-2);
                        componentY += 1;
                        if (componentY == 1){
                            proj2D.push_back(((int) proj1D.size())-2);
                            proj2D.push_back(((int) ycoord.size())-1);
                        }
                    }
                }
            }
            if (j == (LengthY-1)){
                if (componentZ > 0){
                    ycoord.push_back(j);
                }
                else{
                    if (prev_componentZ > 0){
                        ycoord.push_back(j-1);
                    }
                }
            }
            else{
                if (prev_componentZ > 0){
                    if (componentZ == 0){
                        ycoord.push_back(j-1);
                    }
                }
            }
            prev_componentZ = componentZ;
        }
        if (i == 0){
            if (componentY > 0){
                xcoord.push_back(i);
                acc3D.push_back(((int) proj2D.size())-2);
            }
        }
        else{
            if (prev_componentY == 0){
                if (componentY > 0){
                    xcoord.push_back(i);
                    acc3D.push_back(((int) proj2D.size())-2);
                }
            }
        }
        if (i == (LengthX-1)){
            if (componentY > 0){
                xcoord.push_back(i);
            }
            else{
                if (prev_componentY > 0){
                    xcoord.push_back(i-1);
                }
            }
        }
        else{
            if (prev_componentY > 0){
                if (componentY == 0){
                    xcoord.push_back(i-1);
                }
            }
        }
        prev_componentY = componentY;
    }
    std::cout << "DT-Grid structure is generated successfully!" << std::endl;
    return 0;
}


template <typename T1>
int DT_FILE<T1>::write_dtgrdfile_ascii(FILE *output_file)
{
    //Write the uhbd file header first
    std::cout << "Writing the DT-Grid (ASCII) to the output file ..." << std::endl;
    //std::cout << uhbd->Title << std::endl;
    fprintf(output_file, "%72s\n", uhbd->Title);
    fprintf(output_file, "%12.5E%12.5E%7d%7d%7d%7d%7d\n", uhbd->scale, \
    uhbd->dummy1, 3, uhbd->idummy1, uhbd->idummy2, uhbd->one, uhbd->idummy3);
    fprintf(output_file, "%7d%7d%7d%12.5E%12.5E%12.5E%12.5E\n", uhbd->dim.i, uhbd->dim.j, \
    uhbd->dim.k, uhbd->spacing, uhbd->o.x, uhbd->o.y, uhbd->o.z);
    fprintf(output_file, "%12.5E%12.5E%12.5E%12.5E\n", uhbd->dummy2, uhbd->dummy3, uhbd->dummy4, uhbd->dummy5);
    fprintf(output_file, "%12.5E%12.5E%7d%7d\n", uhbd->dummy6, uhbd->dummy7, uhbd->idummy4, uhbd->idummy5);

    write_dtgrdfilebody_ascii(output_file); //call to write the body of the grid file
    return 0;
}


template <typename T1>
int DT_FILE<T1>::write_dtgrdfile_binary(FILE *output_file)
{
    int binary_dt_type = 2;
    int flag_binary = 160;
    //WRITE THE UHBD FILE HEADER FIRST
    std::cout << "Writing the DT-Grid (Binary) to the output file ..." << std::endl;
    fwrite(&flag_binary, sizeof(int),1,output_file);
    fwrite(uhbd->Title,sizeof(char), 72, output_file);
    fwrite(&(uhbd->scale), sizeof(float),1,output_file);
    fwrite(&(uhbd->dummy1), sizeof(float),1,output_file);
    fwrite(&binary_dt_type, sizeof(int),1,output_file);
    fwrite(&(uhbd->idummy1), sizeof(int),1,output_file);
    fwrite(&(uhbd->idummy2), sizeof(int),1,output_file);
    fwrite(&(uhbd->one), sizeof(int),1,output_file);
    fwrite(&(uhbd->idummy3), sizeof(int),1,output_file);
    fwrite(&(uhbd->dim.i), sizeof(int),1,output_file);
    fwrite(&(uhbd->dim.j), sizeof(int),1,output_file);
    fwrite(&(uhbd->dim.k), sizeof(int),1,output_file);
    fwrite(&(uhbd->spacing), sizeof(float),1,output_file);
    fwrite(&(uhbd->o.x), sizeof(float),1,output_file);
    fwrite(&(uhbd->o.y), sizeof(float),1,output_file);
    fwrite(&(uhbd->o.z), sizeof(float),1,output_file);
    fwrite(&(uhbd->dummy2), sizeof(float),1,output_file);
    fwrite(&(uhbd->dummy3), sizeof(float),1,output_file);
    fwrite(&(uhbd->dummy4), sizeof(float),1,output_file);
    fwrite(&(uhbd->dummy5), sizeof(float),1,output_file);
    fwrite(&(uhbd->dummy6), sizeof(float),1,output_file);
    fwrite(&(uhbd->dummy7), sizeof(float),1,output_file);
    fwrite(&(uhbd->idummy4), sizeof(int),1,output_file);
    fwrite(&(uhbd->idummy5), sizeof(int),1,output_file);

    write_dtgrdfilebody_binary(output_file);
    return 0;
}


template <typename T1>
int DT_FILE<T1>::write_dtgrdfile_ascii_nouhbd(FILE *output_file)
{
    char output_file_name [] = "sphere.grd"; // the input file name is not used! why? to be checked. 
    char Title [] = "Sphere                                                                  ";

    float scale = 1.0, spacing = 0.1;
    float dummy1 = 0.0, dummy2 = 0.0, dummy3 = 0.0, dummy4 = 0.0, dummy5 = 0.0, dummy6 = 0.0, dummy7 = 0.0;
    float o_x = 0.0, o_y = 0.0, o_z = 0.0;
    int idummy1 = 0, idummy2 = 0, idummy3 = 0, idummy4 = 0, idummy5 = 0;
    int one = 1, dimi = LengthX, dimj = LengthY, dimk = LengthZ;

    std::cout << "Writing the DT-Grid (ASCII) to the output file ..." << std::endl;
    fprintf(output_file, "%72s\n", Title);
    fprintf(output_file, "%12.5E%12.5E%7d%7d%7d%7d%7d\n", scale, dummy1, 3, idummy1, idummy2, one, idummy3);
    fprintf(output_file, "%7d%7d%7d%12.5E%12.5E%12.5E%12.5E\n", dimi, dimj, dimk, spacing, o_x, o_y, o_z);
    fprintf(output_file, "%12.5E%12.5E%12.5E%12.5E\n", dummy2, dummy3, dummy4, dummy5);
    fprintf(output_file, "%12.5E%12.5E%7d%7d\n", dummy6, dummy7, idummy4, idummy5);

    write_dtgrdfilebody_ascii(output_file);
    return 0;
}


template <typename T1>
int DT_FILE<T1>::write_dtgrdfile_binary_nouhbd(FILE *output_file)
{
    char Title [] = "Sphere                                                                  ";

    float scale = 1.0, spacing = 0.1;
    float dummy1 = 0.0, dummy2 = 0.0, dummy3 = 0.0, dummy4 = 0.0, dummy5 = 0.0, dummy6 = 0.0, dummy7 = 0.0;
    float o_x = 0.0, o_y = 0.0, o_z = 0.0;
    int idummy1 = 0, idummy2 = 0, idummy3 = 0, idummy4 = 0, idummy5 = 0;
    int one = 1, dimi = LengthX, dimj = LengthY, dimk = LengthZ;

    int binary_dt_type = 2;
    int flag_binary = 160;
    std::cout << "Writing the DT-Grid to the output file ..." << std::endl;
    fwrite(&flag_binary, 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_dt_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(&(dimi), sizeof(int),1,output_file);
    fwrite(&(dimj), sizeof(int),1,output_file);
    fwrite(&(dimk), 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_dtgrdfilebody_binary(output_file);
    return 0;
}


template <typename T1>
int DT_FILE<T1>::write_dtgrdfilebody_ascii(FILE *output_file)
{
    //THE DT-GRID HEADER PART
    fprintf(output_file, "%8d%12d\n", grid_type_flag, (int) Value.size());
    fprintf(output_file, "%10d%10d%10d\n", (int) xcoord.size(), (int) ycoord.size(), (int) zcoord.size());
    fprintf(output_file, "%10d%10d%10d\n", (int) acc.size(), (int) acc2D.size(), (int) acc3D.size());
    fprintf(output_file, "%10d%10d\n", ((int) proj1D.size())+2, ((int) proj2D.size())+2);

    // DT-GRID FILE BODY TEXT
    int counter;
    typename std::vector<T1>::iterator itf;
    std::vector<int>::iterator it;

    if (grid_type_flag == 1){
        counter= 0;
        for (itf = Value.begin(); itf != Value.end(); itf++){
            fprintf(output_file, "\t%E", (float) *itf);
            counter++;
            if (counter == 6){
                fprintf(output_file, "\n");
                counter = 0;
            }
        }
        fprintf(output_file, "\n");
    }
    else{
        //no value to be written!
    }

    counter = 0;
    for (it = zcoord.begin(); it != zcoord.end(); it++){
        fprintf(output_file, "\t%d", *it);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }

    counter = 0;
    fprintf(output_file, "\n");
    for (it = acc.begin(); it != acc.end(); it++){
        fprintf(output_file, "\t%d", *it);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }

    counter = 0;
    fprintf(output_file, "\n");
    for (it = proj1D.begin(); it != proj1D.end(); it++){
        if (fmod(counter,2) == 0)
            fprintf(output_file, "\t%d", *it);
        else
            fprintf(output_file, "\t%d", (*it)/2);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }
    
    it--;
    it--;
    fprintf(output_file, "\t%d", *it);
    it++;
    fprintf(output_file, "\t%d", (*it)/2+1);
    
    //fprintf(output_file, "\t%d", (int) Value.size());
    //fprintf(output_file, "\t%d", ((int) zcoord.size())/2 - 1);


    counter = 0;
    fprintf(output_file, "\n");
    for (it = ycoord.begin(); it != ycoord.end(); it++){
        fprintf(output_file, "\t%d", *it);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }

    counter = 0;
    fprintf(output_file, "\n");
    for (it = acc2D.begin(); it != acc2D.end(); it++){
        fprintf(output_file, "\t%d", (*it)/2);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }

    counter = 0;
    fprintf(output_file, "\n");
    for (it = proj2D.begin(); it != proj2D.end(); it++){
        fprintf(output_file, "\t%d", (*it)/2);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }
    
    it--;
    it--;
    fprintf(output_file, "\t%d", (*it)/2);
    it++;
    fprintf(output_file, "\t%d", (*it)/2+1);
    
    //fprintf(output_file, "\t%d", ((int) proj1D.size())/2);
    //fprintf(output_file, "\t%d", ((int)ycoord.size())/2 - 1);

    counter = 0;
    fprintf(output_file, "\n");
    for (it = xcoord.begin(); it != xcoord.end(); it++){
        fprintf(output_file, "\t%d", *it);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }

    counter = 0;
    fprintf(output_file, "\n");
    for (it = acc3D.begin(); it != acc3D.end(); it++){
        fprintf(output_file, "\t%d", (*it)/2);
        counter++;
        if (counter == 12){
            fprintf(output_file, "\n");
            counter = 0;
        }
    }
    fprintf(output_file, "\n");
    return 0;
}


template <typename T1>
int DT_FILE<T1>::write_dtgrdfilebody_binary(FILE *output_file)
{
    //THE DT-GRID HEADER PART
    int dt_header_buffer [10] = {grid_type_flag, (int) Value.size(), (int) xcoord.size(), (int) ycoord.size(), \
            (int) zcoord.size(), (int) acc.size(), (int) acc2D.size(), (int) acc3D.size(), \
            ((int) proj1D.size())+2, ((int) proj2D.size())+2};
    fwrite(dt_header_buffer, sizeof(int), 10, output_file);

    // DT-GRID FILE BODY TEXT
    typename std::vector<T1>::iterator itf;
    std::vector<int>::iterator it;

    int counter;
    float itf_val;
    int it_val;

    if (grid_type_flag == 1){
        itf_val = 0;
        for (itf = Value.begin(); itf != Value.end(); itf++){
            itf_val = (float) *itf;
            fwrite(&itf_val, sizeof(float),1, output_file);
        }
    }
    else{
        //nothing to write for excluded volume grids!
    }

    it_val = 0;
    for (it = zcoord.begin(); it != zcoord.end(); it++){
        it_val = *it;
        fwrite(&it_val, sizeof(int),1, output_file);
    }

    it_val = 0;
    for (it = acc.begin(); it != acc.end(); it++){
        it_val = *it;
        fwrite(&it_val, sizeof(int),1, output_file);
    }

    counter = 0;
    it_val = 0;
    for (it = proj1D.begin(); it != proj1D.end(); it++){
        if (fmod(counter,2) == 0){
            it_val = *it;
            fwrite(&it_val, sizeof(int),1, output_file);
        }
        else{
            it_val = *it / 2;
            fwrite(&it_val, sizeof(int),1, output_file);
        }
        counter++;
    }
    
    /*
     * it--;
     * it--;
     * it_val = *it;
     * it++;
     * it_val = (*it)/2+1;
    */
    
    it_val = (int) Value.size();
    fwrite(&it_val, sizeof(int),1, output_file);
    
    it_val = ((int) zcoord.size() / 2) - 1;
    fwrite(&it_val, sizeof(int),1, output_file);

    it_val = 0;
    for (it = ycoord.begin(); it != ycoord.end(); it++){
        it_val = *it;
        fwrite(&it_val, sizeof(int),1, output_file);
    }

    it_val = 0;
    for (it = acc2D.begin(); it != acc2D.end(); it++){
        it_val = *it / 2;
        fwrite(&it_val, sizeof(int),1, output_file);
    }

    it_val = 0;
    for (it = proj2D.begin(); it != proj2D.end(); it++){
        it_val = *it / 2;
        fwrite(&it_val, sizeof(int),1, output_file);
    }
    /*
     * it--;
     * it--;
     * it_val = *it / 2;
     * it++;
     * it_val = *it / 2 + 1;
    */
    
    it_val = (int) proj1D.size() / 2;
    fwrite(&it_val, sizeof(int),1, output_file);
    
    it_val = ((int) ycoord.size() / 2) - 1;
    fwrite(&it_val, sizeof(int),1, output_file);

    it_val = 0;
    for (it = xcoord.begin(); it != xcoord.end(); it++){
        it_val = *it;
        fwrite(&it_val, sizeof(int),1, output_file);
    }

    it_val = 0;
    for (it = acc3D.begin(); it != acc3D.end(); it++){
        it_val = *it / 2;
        fwrite(&it_val, sizeof(int),1, output_file);
    }

    return 0;
}

#endif /* DTGRD_HPP_ */