finishing ZDOCK implementation

src/rotation.c

@@ -10,18 +10,61 @@ void init_ZDOCK(struct pdb_values** pdbR_addr, struct pdb_values** pdbL_addr){
     struct pdb_values* pdbL = *pdbL_addr;
     struct pdb_values* newPDB = NULL;
     
-    newPDB = rotate_ZDOCK_init(pdbR,newPDB,ZDOCK_REC_X_TRANS,ZDOCK_REC_Y_TRANS,ZDOCK_REC_Z_TRANS,ZDOCK_REC_A_ROT,ZDOCK_REC_B_ROT,ZDOCK_REC_G_ROT);
-    (*pdbR_addr) = newPDB;
-    newPDB = NULL;
+//    newPDB = rotate_ZDOCK_init(pdbR,newPDB);
+//    (*pdbR_addr) = newPDB;
+//    newPDB = NULL;
     // TODO : free pdbR
 
-    newPDB = rotate_ZDOCK_init(pdbL,newPDB,ZDOCK_LIG_X_TRANS,ZDOCK_LIG_Y_TRANS,ZDOCK_LIG_Z_TRANS,ZDOCK_LIG_A_ROT,ZDOCK_LIG_B_ROT,ZDOCK_LIG_G_ROT);
+    newPDB = rotate_ZDOCK_init(pdbL,newPDB);
     (*pdbL_addr) = newPDB;
     newPDB = NULL;
     // TODO : free pdbL
 }
 
 
+/**********************************************/
+/* Function rotateAtom  */
+/*  rotates around 3 euler angles  */
+/**********************************************/
+void rotateAtom (float oldX, float oldY, float oldZ,
+                 float *newX, float *newY, float *newZ,
+                 float alpha, float beta, float gamma, int rev) {
+    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);
+    float r11, r21, r31, r12, r22, r32, r13, r23, r33;
+    if (rev == 0) {
+        r11 = c_a*c_g  -  s_a*c_b*s_g;
+        r21 = s_a*c_g  +  c_a*c_b*s_g;
+        r31 = s_b*s_g;
+        
+        r12 = -c_a*s_g  -  s_a*c_b*c_g;
+        r22 = -s_a*s_g  +  c_a*c_b*c_g;
+        r32 = s_b*c_g;
+
+        r13 = s_a*s_b;
+        r23 = -c_a*s_b;
+        r33 = c_b;
+    }else { // if we are performing a reverse then need the transpose of the matrix
+        r11 = c_a*c_g  -  s_a*c_b*s_g;
+        r12 = s_a*c_g  +  c_a*c_b*s_g;
+        r13 = s_b*s_g;
+        
+        r21 = -c_a*s_g  -  s_a*c_b*c_g;
+        r22 = -s_a*s_g  +  c_a*c_b*c_g;
+        r23 = s_b*c_g;
+
+        r31 = s_a*s_b;
+        r32 = -c_a*s_b;
+        r33 = c_b;
+    }
+
+    *newX = r11 * oldX + r12 * oldY + r13 * oldZ;
+    *newY = r21 * oldX + r22 * oldY + r23 * oldZ;
+    *newZ = r31 * oldX + r32 * oldY + r33 * oldZ;
+
+} /* rotateAtom */
+
 
 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
@@ -42,7 +85,7 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
      * 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 :
+     * The following image available on wikipedia may help to understand the scheme :
      * https://upload.wikimedia.org/wikipedia/commons/thumb/a/a1/Eulerangles.svg/300px-Eulerangles.svg.png
      */
 
@@ -53,20 +96,6 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
     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 : */
-    float xZ = yN*s_b;
-    float yZ = -xN*s_b;
-    float zZ = c_b;
 
     if(newPDB == NULL){
         newPDB = allocate_pdb(pdb->nbRes);
@@ -74,11 +103,14 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
     }
     memcpy(newPDB->residues,pdb->residues,pdb->nbRes*sizeof(struct residue));
 
+    if(trans_X >= ZDOCK_GRID_SIZE/2) trans_X -= ZDOCK_GRID_SIZE;
+    if(trans_Y >= ZDOCK_GRID_SIZE/2) trans_Y -= ZDOCK_GRID_SIZE;
+    if(trans_Z >= ZDOCK_GRID_SIZE/2) trans_Z -= ZDOCK_GRID_SIZE;
 
     /* New center of the protein */
-    float newX = pdb->centerX + (trans_X-(ZDOCK_GRID_SIZE/2))*ZDOCK_GRID_UNIT;
-    float newY = pdb->centerY + (trans_Y-(ZDOCK_GRID_SIZE/2))*ZDOCK_GRID_UNIT;
-    float newZ = pdb->centerZ + (trans_Z-(ZDOCK_GRID_SIZE/2))*ZDOCK_GRID_UNIT;
+    float newX = 0;
+    float newY = 0;
+    float newZ = 0;
 
 
     /* This is now the part where we actually compute the rotations, given
@@ -88,47 +120,61 @@ struct pdb_values* rotate_ZDOCK(struct pdb_values* pdb, struct pdb_values* newPD
     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;
+            xi = pdb->residues[i].x[j];
+            yi = pdb->residues[i].y[j];
+            zi = pdb->residues[i].z[j];
 
+            // rotate for initial randomization and the pose rotation
+            if(ZDOCK_PROT_FIXED == 1){ /* The receptor is fixed */
+                xi += trans_X * ZDOCK_GRID_UNIT;
+                yi += trans_Y * ZDOCK_GRID_UNIT;
+                zi += trans_Z * ZDOCK_GRID_UNIT;
+                rotateAtom(xi, yi, zi, &x1i, &y1i, &z1i, alpha, beta, gamma, ZDOCK_PROT_FIXED);
+                rotateAtom(x1i, y1i, z1i, &x2i, &y2i, &z2i, ZDOCK_LIG_A_ROT, ZDOCK_LIG_B_ROT, ZDOCK_LIG_G_ROT, ZDOCK_PROT_FIXED);
+
+                x2i += ZDOCK_LIG_X_TRANS;
+                y2i += ZDOCK_LIG_Y_TRANS;
+                z2i += ZDOCK_LIG_Z_TRANS;
+
+                newPDB->residues[i].x[j] = x2i;
+                newPDB->residues[i].y[j] = y2i;
+                newPDB->residues[i].z[j] = z2i;
+
+                newX += x2i;
+                newY += y2i;
+                newZ += z2i;
+            }else{ /* The ligand is fixed */
+                rotateAtom(xi, yi, zi, &x1i, &y1i, &z1i, alpha, beta, gamma, ZDOCK_PROT_FIXED);
+                x2i = x1i - trans_X * ZDOCK_GRID_UNIT + ZDOCK_REC_X_TRANS;
+                y2i = y1i - trans_Y * ZDOCK_GRID_UNIT + ZDOCK_REC_Y_TRANS;
+                z2i = z1i - trans_Z * ZDOCK_GRID_UNIT + ZDOCK_REC_Z_TRANS;
+
+                rotateAtom(x2i, y2i, z2i, &x3i, &y3i, &z3i, ZDOCK_REC_A_ROT, ZDOCK_REC_B_ROT, ZDOCK_REC_G_ROT, 1);
+
+                newPDB->residues[i].x[j] = x3i;
+                newPDB->residues[i].y[j] = y3i;
+                newPDB->residues[i].z[j] = z3i;
+
+                newX += x3i;
+                newY += y3i;
+                newZ += z3i;
+            }
         }
     }
 
+    newX /= newPDB->nbAtom;
+    newY /= newPDB->nbAtom;
+    newZ /= newPDB->nbAtom;
+
     newPDB->centerX = newX;
     newPDB->centerY = newY;
     newPDB->centerZ = newZ;
     printf("NEW : %f %f %f\n",newPDB->centerX,newPDB->centerY,newPDB->centerZ);
-    printf("trans : %f %f %f\n",trans_X,trans_Y,trans_Z);
 
     return newPDB;
 }
 
-struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha, float beta, float gamma){
+struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values* newPDB){
     /* 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.
@@ -143,26 +189,6 @@ struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values*
 
     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;
-    printf("N AXIS : %f %f %f\n",xN,yN,zN);
-
-    /* The coordinates of the unit vector for the axis Z can be computed as such :
-     */
-    float xZ = yN*s_b;
-    float yZ = -xN*s_b;
-    float zZ = c_b;
-    printf("Z AXIS : %f %f %f\n",xZ,yZ,zZ);
 
     if(newPDB == NULL){
         newPDB = allocate_pdb(pdb->nbRes);
@@ -176,9 +202,8 @@ struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values*
     float newY = 0;
     float newZ = 0;
 
-
-    /* This is now the part where we actually compute the rotations, given
-     * the coordinates for the rotation axis N and Z
+    /* This is now the part where we actually compute the rotations using
+     * the given initial informations from ZDOCK
      */
     int i = 0, j = 0;
     for(i=0;i<pdb->nbRes;i++){
@@ -188,38 +213,32 @@ struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values*
             yi = pdb->residues[i].y[j];
             zi = pdb->residues[i].z[j];
 
+            if(ZDOCK_PROT_FIXED == 1){ /* means that the receptor is fixed in this case */
+                rotateAtom(xi, yi, zi, &x1i, &y1i, &z1i, ZDOCK_REC_A_ROT, ZDOCK_REC_B_ROT, ZDOCK_REC_G_ROT, 0);
+                x1i -= ZDOCK_REC_X_TRANS;
+                y1i -= ZDOCK_REC_Y_TRANS;
+                z1i -= ZDOCK_REC_Z_TRANS;
+            }else{ /* Here the ligand is fixed */
+                xi -= ZDOCK_LIG_X_TRANS;
+                yi -= ZDOCK_LIG_Y_TRANS;
+                zi -= ZDOCK_LIG_Z_TRANS;
+                rotateAtom(xi, yi, zi, &x1i, &y1i, &z1i, ZDOCK_LIG_A_ROT, ZDOCK_LIG_B_ROT, ZDOCK_LIG_G_ROT, 0);
+            }
             
-            /* 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;
-
+            newPDB->residues[i].x[j] = x1i;
+            newPDB->residues[i].y[j] = y1i;
+            newPDB->residues[i].z[j] = z1i;
 
-            /* Now the residue has been rotated, and we need to translate it
-             * to its final destination
-             */
-            newPDB->residues[i].x[j] = x3i - trans_X;
-            newPDB->residues[i].y[j] = y3i - trans_Y;
-            newPDB->residues[i].z[j] = z3i - trans_Z;
-            newX += newPDB->residues[i].x[j];
-            newY += newPDB->residues[i].y[j];
-            newZ += newPDB->residues[i].z[j];
+            newX += x1i;
+            newY += y1i;
+            newZ += z1i;
         }
     }
 
-    newX /= newPDB->nbAtom;
-    newY /= newPDB->nbAtom;
-    newZ /= newPDB->nbAtom;
+    newX /= pdb->nbAtom;
+    newY /= pdb->nbAtom;
+    newZ /= pdb->nbAtom;
+
     newPDB->centerX = newX;
     newPDB->centerY = newY;
     newPDB->centerZ = newZ;

src/rotation.h

@@ -13,9 +13,13 @@ struct pdb_values* rotate_global(struct pdb_values* pdb, struct pdb_values* pdbR
 
 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);
 
-struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values* newPDB, float trans_X, float trans_Y, float trans_Z, float alpha, float beta, float gamma);
+struct pdb_values* rotate_ZDOCK_init(struct pdb_values* pdb, struct pdb_values* newPDB);
 
 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);
 
 void init_ZDOCK(struct pdb_values** pdbR_addr, struct pdb_values** pdbL_addr);
+
+void rotateAtom (float oldX, float oldY, float oldZ,
+                 float *newX, float *newY, float *newZ,
+                 float alpha, float beta, float gamma, int rev);
 #endif