Adding new format reading (distance follows ID conf). New rotate function, some…

Adding new format reading (distance follows ID conf). New rotate function, some TODO in it before it can be used

src/INTBuilder.c

@@ -236,8 +236,8 @@ void getInterface(struct pdb_values* pdbR, struct pdb_values* pdbL, int confID){
             break;
         }
     }
-    
-    
+
+
     /* If the amount of clashes is still ok */
     if(clash < TOO_MANY_CLASHES){
         /* Output the ligand's interface */

src/fileIO.c

@@ -50,8 +50,14 @@ struct docking_results* getDataForComplex(){
              * The rest of the line doesn't interests us here
              */
             rc = fscanf(condFile_stream,"%d %d %d %d",&idConf,&buf4, &buf4, &buf4);
-        }else{ /* First MAXDo format */
+        }else if(MAXDo){ /* First MAXDo format */
             rc = fscanf(condFile_stream,"%d %d",&idConf,&buf4);
+        }else{
+            /* Else, that means that we have a format where
+             * values for distance between center follows directly
+             * the ID of the conformation.
+             */
+            rc = fscanf(condFile_stream,"%d",&idConf);
         }
 
         /* We confirm that this conformation exists in the MAXDo file */
@@ -67,9 +73,6 @@ struct docking_results* getDataForComplex(){
         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]);
         rc = getline(&buf,&len,condFile_stream); /* Get rid of the rest of the line */
 
-        if(idConf == 137)
-        printf("\n%f %f %f %f %f %f\n",dock_res->distCenters[idConf],dock_res->theta[idConf],dock_res->phi[idConf],dock_res->alpha[idConf],dock_res->beta[idConf],dock_res->gamma[idConf]);
-
         /* theta and phi angles are written in Radian already
          * However that is not the case for alpha, beta and gamma angles
          */
@@ -77,9 +80,14 @@ struct docking_results* getDataForComplex(){
         //dock_res->theta[idConf] = dock_res->theta[idConf]*COEF_DEGREE_TO_RADIAN;
         //dock_res->phi[idConf]  = dock_res->phi[idConf]*COEF_DEGREE_TO_RADIAN;
 
-        dock_res->alpha[idConf] = dock_res->alpha[idConf] * COEF_DEGREE_TO_RADIAN;
-        dock_res->beta[idConf]  = dock_res->beta[idConf]  * COEF_DEGREE_TO_RADIAN;
-        dock_res->gamma[idConf] = dock_res->gamma[idConf] * COEF_DEGREE_TO_RADIAN;
+        /* MAXDo will output the values of alpha, beta and gamma in degrees.
+         * We need to convert these to radian for the rest of the program
+         */
+        if(HCMD2 || MAXDo){
+            dock_res->alpha[idConf] = dock_res->alpha[idConf] * COEF_DEGREE_TO_RADIAN;
+            dock_res->beta[idConf]  = dock_res->beta[idConf]  * COEF_DEGREE_TO_RADIAN;
+            dock_res->gamma[idConf] = dock_res->gamma[idConf] * COEF_DEGREE_TO_RADIAN;
+        }
     }
 
     dock_res->nbConf = numberOfConf;

src/param.c

@@ -29,6 +29,7 @@ void print_usage() {
     "<-pose>............Desired conformation. Default will compute all conformations\n"
     "<--atom-res>.......If specified, will computed the interface at atom level\n"
     "<-HCMD2>...........Means that the docking file is in the HCMD2 MAXDo format\n"
+    "<-MAXDo>...........Means that the docking file is in the original MAXDo format\n"
     "<-noOutput>........Will prevent the program to output the docking interface\n"
     "                   Useful if you only want to use the '-outputPDB' option\n"
     "<--force/-f>.......Will force the computation of the pause even in case of\n"
@@ -52,6 +53,7 @@ struct argLine* parseLineOfArgument(int argc, char** argv){
     constructPDB    = 0;
     target_conf     = -1;
     HCMD2           = 0;
+    MAXDo           = 0;
     complexPDB      = 0;
     doNotOutputINT  = 0;
     clash           = 0;
@@ -119,6 +121,8 @@ struct argLine* parseLineOfArgument(int argc, char** argv){
             }
         }else if(strcmp(argv[i],"-HCMD2") == 0){
             HCMD2 = 1;
+        }else if(strcmp(argv[i],"-MAXDo") == 0){
+            MAXDo = 1;
         }else if(strcmp(argv[i],"--atom-res") == 0){
             atom_res = 1;
         }else if(strcmp(argv[i],"-noOutput") == 0){
@@ -162,6 +166,12 @@ struct argLine* parseLineOfArgument(int argc, char** argv){
             print_usage();
         exit(EXIT_SUCCESS);
     }
+    if(HCMD2 && MAXDo){
+        fprintf(stderr,"You can't choose both HCMD2 and MAXDo format\n\n");
+        if(!print_help)
+            print_usage();
+        exit(EXIT_SUCCESS);
+    }
 
     if(print_help)
         print_usage();

src/rotation.c

@@ -15,9 +15,116 @@ void sphericToxyz(float *x, float *y, float *z, float x0, float y0, float z0, fl
     *z = z0 + R * cos(theta);
 }
 
+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.
+     * Let the N axis the rotation of the x axis around the z axis by an alpha angle.
+     * Let the Z axis the rotation of the z axis around 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 around the axis z, then by a beta
+     * angle around the axis N and finally by a gamma angle around 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);
+
+
+    /* TODO : This needs to be updated to match global settings */
+    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 :
+     * TODO : To be updated as well !
+     *  xZ : -s_a0 * c_b0
+     *  yZ :  c_a0 * c_b0
+     *  zZ :  s_b0
+     */
+    float xZ = 0;
+    float yZ = 0;
+    float zZ = 0;
+
+    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 + R*sin(theta)*cos(phi);
+    float newY = pdbR->centerY + R*sin(theta)*sin(phi);
+    float newZ = pdbR->centerZ + R*cos(theta);
+
+
+    /* 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_sophie(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 docking_results* dock_res, int nConf){
     /* This function rotates the PDB pdb given in argument and sets the new coordinates in newPDB.
-     * The final center of the repair is the geometric center of the PDB pdbR.
+     * The final center of the xyz system is the geometric center of the PDB pdbR.
      * The angles alpha, beta, gamma in arguments correspond to the rotation angles that have to be applied to
      * each residue.
      * In the struct dock_res, we can find the angles alpha0, beta0 and gamma0.

src/rotation.h

@@ -4,7 +4,11 @@
 #define ROTAT_HEADER
 
 struct pdb_values* rotate_sophie(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 docking_results* dock_res, int nConf);
+
 void getAnglesFromCA1andCA5_sophie(struct pdb_values* pdb, float* alpha, float* beta, float* gamma);
+
 void sphericToxyz(float *x, float *y, float *z, float x0, float y0, float z0, float R, float theta, float phi);
 
+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);
+
 #endif

src/struct.h

@@ -3,28 +3,29 @@
 #ifndef STRUCT_HEADER
 #define STRUCT_HEADER
 
-/* Global variables 
- * They are never modified, and are specified 
+/* Global variables
+ * They are never modified, and are specified
  * when the parameters are parsed.
  * They are used to alter the way the program goes
  */
 
 #define COEF_DEGREE_TO_RADIAN 0.017453293   /* Coefficient to convert degrees to radians */
 #define NB_MAX_ATOM_PER_RES 100              /* No residue has more atoms than that */
-#define TOO_MANY_CLASHES 20                /* If not forcing, program will skip conformation after 
+#define TOO_MANY_CLASHES 20                /* If not forcing, program will skip conformation after
                                             * this number of clashes */
-#define DIST_FOR_CONTACT 6                 /* Distance under which we consider 
+#define DIST_FOR_CONTACT 6                 /* Distance under which we consider
                                              * two residues are in contact */
-#define DIST_FOR_CLASH 2           /* Distance in Angstrom under which we consider 
+#define DIST_FOR_CLASH 2           /* Distance in Angstrom under which we consider
                                     * the contact is a clash */
 
 /* 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" 
+#define FORMAT_LINE_PDB "ATOM  %5d  %3s %3s %c%5s   %8.3f%8.3f%8.3f\n"
 
 int     verbose;            /* 1 if we are in verbose mode */
 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     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 */