first commit for HEX

src/INTBuilder.c

@@ -284,10 +284,12 @@ int main(int argc, char** argv){
     pdbR = readPDB(receptor);
     pdbL = readPDB(ligand);
 
+    if(MAXDo || HCMD2){
         /* Compute the angles for the first and fifth CA of the protein
          * This is essential to compute the rotation angles
          */
         getAnglesFromCA1andCA5_sophie(pdbL,&alpha,&beta,&gamma);
+    }
 
     if(complexPDB){    /* If we are building for a given complex, no rotation */
         /* Write the number of the conformation and get the interface */
@@ -298,6 +300,7 @@ int main(int argc, char** argv){
         }
 
     }else if(target_conf > 0){    /* If we are building a specific conformation */
+        if(MAXDo || HCMD2){
             dalpha = dock_res->alpha[target_conf] - alpha;
             dbeta  = dock_res->beta[target_conf]  - beta;
             dgamma = dock_res->gamma[target_conf] - gamma;
@@ -306,6 +309,16 @@ int main(int argc, char** argv){
             newPDB = rotate_sophie(pdbL,pdbR,NULL,dock_res->distCenters[target_conf],
                                 dock_res->theta[target_conf],dock_res->phi[target_conf],
                                 dalpha,dbeta,dgamma,dock_res,target_conf);
+        }else if(HEX){
+            newPDB = rotate_HEX(pdbL,pdbR,NULL,dock_res->trans_X[target_conf],
+                                dock_res->trans_Y[target_conf],dock_res->trans_Z[target_conf],
+                                dock_res->alpha[target_conf],dock_res->beta[target_conf],
+                                dock_res->gamma[target_conf]);
+        }else{
+            newPDB = rotate_global(pdbL,pdbR,NULL,dock_res->distCenters[target_conf],
+                                dock_res->theta[target_conf],dock_res->phi[target_conf],
+                                alpha,beta,gamma);
+        }
 
         /* Write the number of the conformation and get the interface */
         if(!doNotOutputINT){
@@ -325,6 +338,7 @@ int main(int argc, char** argv){
                 continue;
             }
 
+            if(MAXDo || HCMD2){
                 /* compute the rotation angles dalpha, dbeta and dgamma
                  * from the angles alpha and beta of the CA 1 and the
                  * angle gamma computed after rotation of alpha and beta
@@ -340,6 +354,16 @@ int main(int argc, char** argv){
                 newPDB = rotate_sophie(pdbL, pdbR, newPDB, dock_res->distCenters[i],
                                     dock_res->theta[i], dock_res->phi[i],
                                     dalpha, dbeta, dgamma, dock_res, i);
+            }else if(HEX){
+                newPDB = rotate_HEX(pdbL,pdbR,NULL,dock_res->trans_X[i],
+                                    dock_res->trans_Y[i],dock_res->trans_Z[i],
+                                    dock_res->alpha[i],dock_res->beta[i],
+                                    dock_res->gamma[i]);
+            }else{
+                newPDB = rotate_global(pdbL,pdbR,NULL,dock_res->distCenters[i],
+                                    dock_res->theta[i],dock_res->phi[i],
+                                    alpha,beta,gamma);
+            }
 
             /* Write the number of the conformation and get the interface */
             if(!doNotOutputINT){

src/allocate.c

@@ -47,12 +47,57 @@ struct docking_results* allocate_dockingResults(int nbConf){
     dock_res->beta        = NULL;
     dock_res->gamma       = NULL;
     dock_res->t_conf      = NULL;
+    dock_res->trans_X     = NULL;
+    dock_res->trans_Y     = NULL;
+    dock_res->trans_Z     = NULL;
 
     dock_res->listEner = malloc((2*nbConf)*sizeof(float));
     if(dock_res->listEner == NULL){
         perror("malloc");
         exit(EXIT_FAILURE);
     }
+    dock_res->alpha = malloc((2*nbConf)*sizeof(float));
+    if(dock_res->alpha == NULL){
+        perror("malloc");
+        exit(EXIT_FAILURE);
+    }
+    dock_res->beta = malloc((2*nbConf)*sizeof(float));
+    if(dock_res->beta == NULL){
+        perror("malloc");
+        exit(EXIT_FAILURE);
+    }
+    dock_res->gamma = malloc((2*nbConf)*sizeof(float));
+    if(dock_res->gamma == NULL){
+        perror("malloc");
+        exit(EXIT_FAILURE);
+    }
+
+    dock_res->t_conf = malloc((2*nbConf)*sizeof(int));
+    if(dock_res->t_conf == NULL){
+        perror("malloc");
+        exit(EXIT_FAILURE);
+    }
+
+    /* Only HEX uses translation values */
+    if(HEX){
+        dock_res->trans_X = malloc((2*nbConf)*sizeof(int));
+        if(dock_res->trans_X == NULL){
+            perror("malloc");
+            exit(EXIT_FAILURE);
+        }
+
+        dock_res->trans_Y = malloc((2*nbConf)*sizeof(int));
+        if(dock_res->trans_Y == NULL){
+            perror("malloc");
+            exit(EXIT_FAILURE);
+        }
+        
+        dock_res->trans_Z = malloc((2*nbConf)*sizeof(int));
+        if(dock_res->trans_Z == NULL){
+            perror("malloc");
+            exit(EXIT_FAILURE);
+        }
+    }else{
         dock_res->distCenters = malloc((2*nbConf)*sizeof(float));
         if(dock_res->distCenters == NULL){
             perror("malloc");
@@ -68,37 +113,39 @@ struct docking_results* allocate_dockingResults(int nbConf){
             perror("malloc");
             exit(EXIT_FAILURE);
         }
-    dock_res->alpha = malloc((2*nbConf)*sizeof(float));
-    if(dock_res->alpha == NULL){
-        perror("malloc");
-        exit(EXIT_FAILURE);
     }
-    dock_res->beta = malloc((2*nbConf)*sizeof(float));
-    if(dock_res->beta == NULL){
+    dock_res->distCenters = malloc((2*nbConf)*sizeof(float));
+    if(dock_res->distCenters == NULL){
         perror("malloc");
         exit(EXIT_FAILURE);
     }
-    dock_res->gamma = malloc((2*nbConf)*sizeof(float));
-    if(dock_res->gamma == NULL){
+    dock_res->theta = malloc((2*nbConf)*sizeof(float));
+    if(dock_res->theta == NULL){
         perror("malloc");
         exit(EXIT_FAILURE);
     }
-
-    dock_res->t_conf = malloc((2*nbConf)*sizeof(int));
-    if(dock_res->t_conf == NULL){
+    dock_res->phi = malloc((2*nbConf)*sizeof(float));
+    if(dock_res->phi == NULL){
         perror("malloc");
         exit(EXIT_FAILURE);
     }
 
     for(i=0;i<nbConf;i++){
         dock_res->listEner[i]    = 0;
-        dock_res->distCenters[i] = 0;
-        dock_res->theta[i]       = 0;
-        dock_res->phi[i]         = 0;
+
         dock_res->alpha[i]       = 0;
         dock_res->beta[i]        = 0;
         dock_res->gamma[i]       = 0;
         dock_res->t_conf[i]      = 0;
+        if(HEX){
+            dock_res->trans_X    = 0;
+            dock_res->trans_Y    = 0;
+            dock_res->trans_Z    = 0;
+        }else{
+            dock_res->distCenters[i] = 0;
+            dock_res->theta[i]       = 0;
+            dock_res->phi[i]         = 0;
+        }
     }
 
 

src/fileIO.c

@@ -9,21 +9,34 @@ struct docking_results* getDataForComplex(){
     char* buf = NULL;
     FILE* condFile_stream = NULL;
     int numberOfConf = 0, rc = 0;
-    int idConf = 0, buf4 = 0;
+    int idConf = 0, buf4 = 0, i = 0;
     float buf3 = 0;
     size_t len = 2000;
     struct docking_results* dock_res = NULL;
 
     condFile_stream = fopen(dockingFile,"r");
 
+    if(HEX){
+        for(i=0;i<NB_HEADER_LINE_HEX;i++){
+            getline(&buf,&len,condFile_stream);
+        }
+    }
+
     /* Get the number of conformations */
-    rc = fscanf(condFile_stream,"%f",&buf3);
+    if(HEX){
+        rc = fscanf(condFile_stream,"%f %f",&buf3,&buf3);
+    }else{
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
+    }
     if(buf3 > numberOfConf)
         numberOfConf = (int) buf3;
 
     while(!feof(condFile_stream)){
+        if(HEX){
+            rc = fscanf(condFile_stream,"%f %f",&buf3,&buf3);
+        }else{
             rc = fscanf(condFile_stream,"%f",&buf3);
+        }
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
         if(buf3 > numberOfConf)
             numberOfConf = (int) buf3;
@@ -32,8 +45,12 @@ struct docking_results* getDataForComplex(){
     rewind(condFile_stream);
 
     /* This is with the first MAXDo format */
-    if(!HCMD2){
+    if(MAXDo){
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the first line */
+    }else if(HEX){
+        for(i=0;i<NB_HEADER_LINE_HEX;i++){
+            getline(&buf,&len,condFile_stream);
+        }
     }
 
     /* Now that we know how many conformations we are dealing with
@@ -46,6 +63,8 @@ struct docking_results* getDataForComplex(){
 
         if(MAXDo){ /* First MAXDo format */
             rc = fscanf(condFile_stream,"%d %d",&idConf,&buf4);
+        }else if(HEX){
+            rc = fscanf(condFile_stream,"%d %d",&buf4,&idConf);
         }else{
             /* Else, that means that we have a format where
              * values for distance between center follows directly
@@ -64,7 +83,12 @@ struct docking_results* getDataForComplex(){
         /* printf("curConf : %d\n",idConf); */
 
 
+        if(HEX){
+            rc = fscanf(condFile_stream,"%f %f %f %f %f %f",&dock_res->alpha[idConf],&dock_res->beta[idConf],&dock_res->gamma[idConf],&dock_res->trans_X[idConf],&dock_res->trans_Y[idConf],&dock_res->trans_Z[idConf]);
+        }else{
             rc = fscanf(condFile_stream,"%f %f %f %f %f %f",&dock_res->distCenters[idConf],&dock_res->theta[idConf],&dock_res->phi[idConf],&dock_res->alpha[idConf],&dock_res->beta[idConf],&dock_res->gamma[idConf]);
+        }
+        // TODO : Is this necessary for HEX ?
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
 
         /* theta and phi angles are written in Radian already

src/param.c

@@ -123,6 +123,8 @@ struct argLine* parseLineOfArgument(int argc, char** argv){
             HCMD2 = 1;
         }else if(strcmp(argv[i],"-MAXDo") == 0){
             MAXDo = 1;
+        }else if(strcmp(argv[i],"-HEX") == 0){
+            HEX = 1;
         }else if(strcmp(argv[i],"--atom-res") == 0){
             atom_res = 1;
         }else if(strcmp(argv[i],"-noOutput") == 0){

src/rotation.c

@@ -15,6 +15,110 @@ void sphericToxyz(float *x, float *y, float *z, float x0, float y0, float z0, fl
     *z = z0 + R * cos(theta);
 }
 
+struct pdb_values* rotate_HEX(struct pdb_values* pdb, struct pdb_values* pdbR, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha,float beta,float gamma){
+    /* This function aims at rotating the pdb parameter and gets its new coordinates
+     * in the newPDB parameter.
+     * Let the N axis the rotation of the x axis about the z axis by an alpha angle.
+     * Let the Z axis the rotation of the z axis about the N axis by a beta angle.
+     *
+     * The xyz system is centered on the center of mass of the pdb parameter.
+     * We then rotate the pdb by an alpha angle about the axis z, then by a beta
+     * angle about the axis N and finally by a gamma angle about the axis Z.
+     * The following image available on wikipedia can help to understand the scheme :
+     * https://upload.wikimedia.org/wikipedia/commons/thumb/a/a1/Eulerangles.svg/300px-Eulerangles.svg.png
+     */
+
+    float xi = 0, yi = 0, zi = 0;
+    float x1i = 0, y1i = 0, z1i = 0;
+    float x2i = 0, y2i = 0, z2i = 0;
+    float x3i = 0, y3i = 0, z3i = 0;
+
+    float c_a=cos(alpha),s_a=sin(alpha);
+    float c_b=cos(beta),s_b=sin(beta);
+    float c_g=cos(gamma),s_g=sin(gamma);
+
+
+    /* Coordinates of the N axis */
+    float xN = c_a;
+    float yN = s_a;
+    float zN = 0;
+
+    /* The coordinates of the unit vector for the axis Z can be computed as such :
+     *  xZ : c_b
+     *  yZ : c_a * s_b
+     *  zZ : s_a * s_b
+     */
+    float xZ = c_b;
+    float yZ = c_a * s_b;
+    float zZ = s_a * s_b;
+
+    if(newPDB == NULL){
+        newPDB = allocate_pdb(pdb->nbRes);
+        newPDB->nbAtom = pdb->nbAtom;
+    }
+    memcpy(newPDB->residues,pdb->residues,pdb->nbRes*sizeof(struct residue));
+
+
+    float newX = pdbR->centerX + trans_X;
+    float newY = pdbR->centerY + trans_Y;
+    float newZ = pdbR->centerZ + trans_Z;
+
+    /* This is now the part where we actually compute the rotations, given
+     * the coordinates for the rotation axis N and Z
+     */
+    int i = 0, j = 0;
+    for(i=0;i<pdb->nbRes;i++){
+        for(j=0;j<pdb->residues[i].nbAtom;j++){
+
+            xi = pdb->residues[i].x[j] - pdb->centerX;
+            yi = pdb->residues[i].y[j] - pdb->centerY;
+            zi = pdb->residues[i].z[j] - pdb->centerZ;
+
+
+            /* alpha rotation of the residue */
+            x1i=c_a*xi-s_a*yi;
+            y1i=c_a*yi+s_a*xi;
+            z1i=zi;
+
+            /* beta rotation of the residue */
+            x2i = (xN*xN+(1.0-xN*xN)*c_b)*x1i + xN*yN*(1.0-c_b)*y1i         + yN*s_b*z1i;
+            y2i = xN*yN*(1.0-c_b)*x1i         + (yN*yN+(1.0-yN*yN)*c_b)*y1i - xN*s_b*z1i;
+            z2i = -yN*s_b*x1i                 + xN*s_b*y1i                  + c_b*z1i;
+
+            /* gamma rotation of the residue */
+            x3i = (xZ*xZ+(1.0-xZ*xZ)*c_g) *x2i + (xZ*yZ*(1.0-c_g)-zZ*s_g)*y2i + (xZ*zZ*(1.0-c_g)+yZ*s_g)*z2i;
+            y3i = (xZ*yZ*(1.0-c_g)+zZ*s_g)*x2i + (yZ*yZ+(1.0-yZ*yZ)*c_g) *y2i + (yZ*zZ*(1.0-c_g)-xZ*s_g)*z2i;
+            z3i = (xZ*zZ*(1.0-c_g)-yZ*s_g)*x2i + (yZ*zZ*(1.0-c_g)+xZ*s_g)*y2i + (zZ*zZ+(1.0-zZ*zZ)*c_g) *z2i;
+
+
+            /* Now the residue has been rotated, and we need to translate it
+             * to its final destination
+             */
+            newPDB->residues[i].x[j] = x3i + newX;
+            newPDB->residues[i].y[j] = y3i + newY;
+            newPDB->residues[i].z[j] = z3i + newZ;
+
+            /* We also translate the coordinates of the first and fifth CA, even
+             * though they are not used afterwards. This is for the sake of
+             * having the right values at the right places
+             */
+            newPDB->xCA1 = (newPDB->xCA1 - pdb->centerX) + newX;
+            newPDB->yCA1 = (newPDB->yCA1 - pdb->centerY) + newY;
+            newPDB->zCA1 = (newPDB->zCA1 - pdb->centerZ) + newZ;
+
+            newPDB->xCA5 = (newPDB->xCA5 - pdb->centerX) + newX;
+            newPDB->yCA5 = (newPDB->yCA5 - pdb->centerY) + newY;
+            newPDB->zCA5 = (newPDB->zCA5 - pdb->centerZ) + newZ;
+
+            newPDB->centerX = newX;
+            newPDB->centerY = newY;
+            newPDB->centerZ = newZ;
+        }
+    }
+
+    return newPDB;
+}
+
 struct pdb_values* rotate_global(struct pdb_values* pdb, struct pdb_values* pdbR, struct pdb_values* newPDB, float R, float theta, float phi, float alpha,float beta,float gamma){
     /* This function aims at rotating the pdb parameter and gets its new coordinates
      * in the newPDB parameter.

src/rotation.h

@@ -11,4 +11,5 @@ void sphericToxyz(float *x, float *y, float *z, float x0, float y0, float z0, fl
 
 struct pdb_values* rotate_global(struct pdb_values* pdb, struct pdb_values* pdbR, struct pdb_values* newPDB, float R, float theta, float phi, float alpha,float beta,float gamma);
 
+struct pdb_values* rotate_HEX(struct pdb_values* pdb, struct pdb_values* pdbR, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha,float beta,float gamma);
 #endif

src/struct.h

@@ -17,6 +17,7 @@
                                              * two residues are in contact */
 #define DIST_FOR_CLASH 2           /* Distance in Angstrom under which we consider
                                     * the contact is a clash */
+#define NB_HEADER_LINE_HEX 5        /* Number of header lines at the top of a hex docking file */
 
 /* This is the format of a line in a PDB file */
 #define FORMAT_LINE_PDB "ATOM  %5d  %3s %3s %c%5s   %8.3f%8.3f%8.3f\n"
@@ -26,6 +27,7 @@ int     constructPDB;       /* 1 if we should build the PDB */
 int     target_conf;        /* 1 if we are specific to one conformation */
 int     HCMD2;              /* 1 if this is HCMD2 format */
 int     MAXDo;              /* 1 if this is MAXDo format */
+int     HEX;                /* 1 if this is HEX format */
 int     complexPDB;         /* 1 if we are computing the interface of a given complex */
 int     doNotOutputINT;     /* 1 if we don't want to output the interface */
 int     atom_res;           /* 1 if we want a resolution at atom's level */
@@ -81,6 +83,9 @@ struct docking_results {
     float*              alpha;          /* alpha angle */
     float*              beta;           /* beta angle */
     float*              gamma;          /* Gamma angle */
+    float*              trans_X;
+    float*              trans_Y;
+    float*              trans_Z;
     int*                t_conf;         /* 1 for a conformation position if it exists, 0 otherwise */
     int                 nbConf;         /* Number of conformations */
 };