ZDOCK : bug correction and rotation function

src/fileIO.c

@@ -24,10 +24,10 @@ struct docking_results* getDataForComplex(){
 
     if(ZDOCK){
         getline(&buf,&len,condFile_stream);
-        fscanf(condFile_stream,"%f %f %f",&ZDOCK_REC_X_ROT,&ZDOCK_REC_X_ROT,&ZDOCK_REC_X_ROT);
-        fscanf(condFile_stream,"%f %f %f",&ZDOCK_LIG_X_ROT,&ZDOCK_LIG_X_ROT,&ZDOCK_LIG_X_ROT);
-        fscanf(condFile_stream,"%s %f %f %f",buf,&ZDOCK_REC_X_TRANS,&ZDOCK_REC_X_TRANS,&ZDOCK_REC_X_TRANS);
-        fscanf(condFile_stream,"%s %f %f %f",buf,&ZDOCK_LIG_X_TRANS,&ZDOCK_LIG_X_TRANS,&ZDOCK_LIG_X_TRANS);
+        fscanf(condFile_stream,"%f %f %f",&ZDOCK_REC_A_ROT,&ZDOCK_REC_B_ROT,&ZDOCK_REC_G_ROT);
+        fscanf(condFile_stream,"%f %f %f",&ZDOCK_LIG_A_ROT,&ZDOCK_LIG_B_ROT,&ZDOCK_LIG_G_ROT);
+        fscanf(condFile_stream,"%s %f %f %f",buf,&ZDOCK_REC_X_TRANS,&ZDOCK_REC_Y_TRANS,&ZDOCK_REC_Z_TRANS);
+        fscanf(condFile_stream,"%s %f %f %f",buf,&ZDOCK_LIG_X_TRANS,&ZDOCK_LIG_Y_TRANS,&ZDOCK_LIG_Z_TRANS);
     }
 
     /* Get the number of conformations */

src/rotation.c

@@ -6,6 +6,8 @@
 
 
 
+
+
 void sphericToxyz(float *x, float *y, float *z, float x0, float y0, float z0, float R, float theta, float phi){
     /* Sets the values of x, y and z to the carthesian equivalent of the spherical
      * coordinate given by R, theta, phi with a center x0, y0 and z0
@@ -15,6 +17,102 @@ void sphericToxyz(float *x, float *y, float *z, float x0, float y0, float z0, fl
     *z = z0 + R * cos(theta);
 }
 
+
+struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, 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));
+
+
+    /* New center of the protein */
+    float newX = pdb->centerX + trans_X;
+    float newY = pdb->centerY + trans_Y;
+    float newZ = pdb->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;
+
+        }
+    }
+
+    newPDB->centerX = newX;
+    newPDB->centerY = newY;
+    newPDB->centerZ = newZ;
+
+    return newPDB;
+}
+
 struct pdb_values* rotate_HEX(struct pdb_values* pdb, 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.
@@ -64,11 +162,11 @@ struct pdb_values* rotate_HEX(struct pdb_values* pdb, struct pdb_values* newPDB,
             newPDB->residues[i].y[j] = y3i + trans_Y;
             newPDB->residues[i].z[j] = z3i + trans_Z;
 
-            newPDB->centerX += trans_X;
-            newPDB->centerY += trans_X;
-            newPDB->centerZ += trans_Z;
         }
     }
+    newPDB->centerX += trans_X;
+    newPDB->centerY += trans_Y;
+    newPDB->centerZ += trans_Z;
 
     return newPDB;
 }

src/struct.h

@@ -33,12 +33,12 @@ int     atom_res;           /* 1 if we want a resolution at atom's level */
 int     clash;              /* Maybe be incremented in case of clashes */
 int     force;              /* Will output even if clash, may increase computation time */
 float   DIST_FOR_CONTACT;   /* Distance under which we consider two residues are in contact */
-float   ZDOCK_LIG_X_ROT;    /* Initial rotations and translations if using ZDOCK */
-float   ZDOCK_LIG_Y_ROT;
-float   ZDOCK_LIG_Z_ROT;
-float   ZDOCK_REC_X_ROT;
-float   ZDOCK_REC_Y_ROT;
-float   ZDOCK_REC_Z_ROT;
+float   ZDOCK_LIG_A_ROT;    /* Initial rotations and translations if using ZDOCK */
+float   ZDOCK_LIG_B_ROT;
+float   ZDOCK_LIG_G_ROT;
+float   ZDOCK_REC_A_ROT;
+float   ZDOCK_REC_B_ROT;
+float   ZDOCK_REC_G_ROT;
 float   ZDOCK_LIG_X_TRANS;
 float   ZDOCK_LIG_Y_TRANS;
 float   ZDOCK_LIG_Z_TRANS;