src/Tools/pdbUtil/pdb2mrc2d/src/pdb2mrc2d.c

   1 /*
   2 # pdb2mrc2d.c  1.5
   3 # The latest update : 11/02/96 at 23:15:49
   4 #
   5 #@(#) pdb2mrc2d ver 1.5
   6 #@(#) Created by
   7 #@(#)
   8 #@(#) Usage : pdb2mrc2d
   9 #@(#) Attention
  10 #@(#)
  11 */
  12
  13 static char __sccs_id[] = "@(#)pdb2mrc2d ver1.5; Date:96/11/02 @(#)";
  14
  15 /*
  16
  17       y
  18       |
  19       |
  20       |
  21       |
  22       |____________ x
  23      /
  24     /
  25    /
  26  z
  27
  28         Projection : x-y plane along z-axis
  29
  30         RotationRule    : YOYS
  31         First Rotation  : y-axis : Y: Y
  32         Second Rotation : x-axis : O: Odd
  33         Last Rotation   : z-axis : Y: Same
  34         v1 = A v0                : S: Staic
  35
  36 */
  37 #include <stdio.h>
  38 #include <stdlib.h>
  39 #include <string.h>
  40 #include <math.h>
  41 #undef DEBUG
  42 #include "mrcImage.h"
  43 #include "../inc/config.h"
  44 #include "pdbFile.h"
  45 #include "lpdb2mrcInfo.h"
  46 #include "genUtil.h"
  47 #include "Matrix3D.h"
  48 #include "Map2D.h"
  49
  50 #ifdef CUDA
  51 #include "eosCuda.h"
  52 #endif
  53
  54 #define HowToCreateImageBit (0x1)
  55 #define HowToDrawMapBit     (0x2)
  56
  57 int
  58 main(int argc, char* argv[])
  59 {
  60     pdb2mrc2dInfo info;
  61     lpdb2mrc2dInfo linfo;
  62     pdbFile pdb;
  63     mrcImage mrc;
  64     mrcImage out;
  65     mrcImage out3D;
  66         mrcImageParaTypeRealCoord to;
  67         mrcImageParaTypeReal ix, iy;
  68         pdbFileParaTypeReal rotx, roty, rotz, droty;
  69         Map2DParaTypeInteger nrotx, nroty;
  70         mrcImageParaTypeReal zsection;
  71
  72     init0(&info);
  73     argCheck(&info, argc, argv);
  74     init1(&info);
  75
  76 #ifdef CUDA
  77         if(info.flagcudaDeviceID) {
  78                 eosCudaInit(info.cudaDeviceID);
  79         }
  80 #endif
  81
  82         linfo.dx = info.dx; linfo.dy = info.dy;
  83         linfo.nx = info.nx; linfo.ny = info.ny;
  84
  85         if(info.flagsx && info.flagsy ) {
  86                 linfo.sx = info.sx; linfo.sy = info.sy;
  87         } else {
  88                 linfo.sx = -1.0*info.dx*info.nx/2.0 ;
  89                 linfo.sy = -1.0*info.dy*info.ny/2.0 ;
  90         }
  91         linfo.Weight = info.Weight;
  92         linfo.Sigma = info.Sigma;
  93
  94         if(!(info.flagOut|info.flagOut3D)) {
  95                 fprintf(stderr, "This program has no work for your request\n");
  96                 fprintf(stderr, "-o or -O option is required.\n");
  97                 usage(argv[0]);
  98         }
  99
 100     pdbFileRead(info.fptIn, &pdb);
 101         if(info.flagsrotx) {
 102                 pdbRotationZXY(&pdb, info.srotx*RADIAN, -info.sroty*RADIAN, info.srotz*RADIAN);
 103         }
 104         if(info.flagsRotMode) {
 105                 pdbRotationFollowingEulerAngle(&pdb, info.sRotMode, info.sRot1*RADIAN, -info.sRot2*RADIAN, info.sRot3*RADIAN);
 106         }
 107         if(info.flagdRot1 && !info.flagrotnx) {
 108                 info.rotnx = (int)((info.maxRot1 - info.minRot1)/info.dRot1+0.5);
 109         }
 110         if(info.flagdRot2 && !info.flagrotny) {
 111                 info.rotny = (int)((info.maxRot2 - info.minRot2)/info.dRot2+0.5)*2;
 112         }
 113         if(info.flagOut) {
 114                 out.HeaderN.x = info.rotnx*info.nx;
 115                 out.HeaderN.y = info.rotny*info.ny;
 116                 out.HeaderN.z = 1.0;
 117                 out.HeaderMode =        mrcFloatImage;
 118                 out.HeaderLength.x = info.dx;
 119                 out.HeaderLength.y = info.dy;
 120                 out.HeaderLength.z = 0.0;
 121             mrcInit(&out, (char*)NULL);
 122         }
 123         if(info.flagOut3D) {
 124                 out3D.HeaderN.x = info.nx;
 125                 out3D.HeaderN.y = info.ny;
 126                 out3D.HeaderN.z = info.rotnx*info.rotny;
 127                 out3D.HeaderMode =      mrcFloatImage;
 128                 out3D.HeaderLength.x = info.dx;
 129                 out3D.HeaderLength.y = info.dy;
 130                 out3D.HeaderLength.z = 0.0;
 131             mrcInit(&out3D, NULL);
 132                 out3D.numTailer = out3D.HeaderN.z;
 133                 mrcTailerInit(&out3D, 0);
 134         }
 135         rotz = 0.0;
 136         nrotx = info.rotny; /* Attention: x - y change */
 137
 138         zsection = 0.0;
 139         for(ix=0; ix < nrotx ; ix++) {
 140                 if(nrotx/4 < ix && ix < nrotx/4*3) {
 141                         /* Skip */
 142                 } else {
 143                         if(!info.flagRotMode) {
 144                                 rotx = ((double)ix)/nrotx*2.0*M_PI;
 145                         } else {
 146                                 rotx = ((double)ix)/nrotx*2.0*M_PI;
 147                         }
 148                         map2DParallelInfo(&nroty, &droty,
 149                                                           rotx, info.rotnx,
 150                                                           (info.mode&HowToCreateImageBit));
 151                         for(iy=0; iy < info.rotny; iy+= droty) {
 152                                 DEBUGPRINT2("%g %g\n", ix, iy);
 153
 154 #ifdef DEBUG
 155                                 {
 156                                         pdbFileTop(&pdb);
 157                                         pdbFileNextAtom(&pdb);
 158                                         DEBUGPRINT3("(x, y, z)=(%g, %g %g)\n", pdb.PDB->Coord.x, pdb.PDB->Coord.y, pdb.PDB->Coord.z);
 159                                 }
 160 #endif
 161
 162                                 if(!info.flagRotMode) {
 163                                         roty = ((double)iy)/info.rotny*2.0*M_PI;
 164                                         pdbRotationZXY(&pdb, rotx, -roty, -rotz);
 165                                 } else {
 166                                         roty = ((double)iy)/info.rotny*2.0*M_PI;
 167                                         pdbRotationFollowingEulerAngle(&pdb, info.RotMode, -roty, rotx, -rotz);
 168                                 }
 169
 170 #ifdef DEBUG
 171                                 {
 172                                         pdbFileTop(&pdb);
 173                                         pdbFileNextAtom(&pdb);
 174                                         DEBUGPRINT3("(x, y, z)=(%g, %g %g)\n", pdb.PDB->Coord.x, pdb.PDB->Coord.y, pdb.PDB->Coord.z);
 175                                 }
 176 #endif
 177
 178                                 DEBUGPRINT3("rot: (%g %g %g) \n", rotx*DEGREE, roty*DEGREE, rotz*DEGREE);
 179
 180                         lpdb2mrc2d(&mrc, &pdb,  &linfo, 'z', 0);
 181 #ifdef DEBUG
 182                                 {
 183                                         char s[1024];
 184
 185                                         sprintf(s, "/tmp/test.mrc.%03d", (int)zsection);
 186
 187                                         mrcFileWrite(&mrc, s, "tmp", 0);
 188                                 }
 189 #endif
 190
 191                                 if(!info.flagRotMode) {
 192                                         pdbRotationYXZ(&pdb, -rotx, roty, rotz);
 193                                 } else {
 194                                         pdbRotationFollowingEulerAngleInverse(&pdb, info.RotMode, roty, -rotx, -rotz);
 195                                 }
 196
 197                                 DEBUGPRINT3("%g %g %g\n", -roty*DEGREE, rotx*DEGREE, -rotz*DEGREE);
 198                                 if(info.flagOut) {
 199                                         map2DCoordGet(&(to.x), &(to.y), roty, rotx, info.rotnx, info.rotny, ((info.mode&HowToDrawMapBit)>>1));
 200                                         to.x = fmod(to.x*info.nx, out.HeaderN.x);
 201                                         to.y = fmod(to.y*info.ny, out.HeaderN.y);
 202                                         DEBUGPRINT2("map %f %f\n", to.x, to.y);
 203                                         to.z = 0.0;
 204                                         DEBUGPRINT3("to: %f %f %f\n", to.x, to.y, to.z);
 205                                         lmrcImageCopy(&out, &mrc, to);
 206 #ifdef DEBUG
 207                                 {
 208                                         char s[1024];
 209                                         sprintf(s, "/tmp/test.2d.%03d", (int)zsection);
 210                                         mrcFileWrite(&out, s, "tmp", 0);
 211                                 }
 212 #endif
 213                                 }
 214                                 if(info.flagOut3D) {
 215                                         to.x = 0.0;
 216                                         to.y = 0.0;
 217                                         to.z = zsection;
 218                                         out3D.Tailer[(int)zsection].Cont.Mode = mrcImageTailerMode2DProjection;
 219                                         out3D.Tailer[(int)zsection].Cont.EulerAngleMode[0] = info.RotMode[0]; /* X or Y or Z */
 220                                         out3D.Tailer[(int)zsection].Cont.EulerAngleMode[1] = info.RotMode[1]; /* E[ven] or O[dd] */
 221                                         out3D.Tailer[(int)zsection].Cont.EulerAngleMode[2] = info.RotMode[2]; /* S[ame] or D[iff] */
 222                                         out3D.Tailer[(int)zsection].Cont.EulerAngleMode[3] = info.RotMode[3]; /* S[tatic] or R[otating] */
 223                                         out3D.Tailer[(int)zsection].Cont.Rot1 = -roty;
 224                                         out3D.Tailer[(int)zsection].Cont.Rot2 =  rotx;
 225                                         out3D.Tailer[(int)zsection].Cont.Rot3 = -rotz;
 226                                         DEBUGPRINT3("to: %f %f %f\n", to.x, to.y, to.z);
 227                                         lmrcImageCopy(&out3D, &mrc, to);
 228                                         zsection += 1.0;
 229                                 }
 230                                 mrcImageFree(&mrc, 0);
 231                         }
 232                 }
 233         }
 234         if(info.flagOut) {
 235            mrcStatDataSet(&out, 0);
 236            mrcFileWrite(&out, info.Out, "in pdb2mrc2d Main Routine", 0);
 237         }
 238         if(info.flagOut3D) {
 239            out3D.HeaderN.z = zsection;
 240            out3D.numTailer = zsection;
 241            mrcHiddenDataSet(&out3D, 0);
 242            mrcStatDataSet(&out3D, 0);
 243            mrcFileWrite(&out3D, info.Out3D, "in pdb2mrc2d Main Routine", 0);
 244         }
 245         return 0;
 246 }
 247
 248 void
 249 additionalUsage()
 250 {
 251         fprintf(stderr, "----- Attention1 -----\n");
 252         fprintf(stderr, "If both of -s and -startEA, first -s and second -startEA will be performed\n");
 253         fprintf(stderr, "----- Attention2 -----\n");
 254         fprintf(stderr, "Bug fixed: rotation around x-axis. + is changed to - \n");
 255         fprintf(stderr, "----- Mode -----\n");
 256         fprintf(stderr, "%d: 0: equal angle 1: equal area\n", HowToCreateImageBit);
 257         fprintf(stderr, "%d: 0: Mercatol    1: Morwide \n", HowToDrawMapBit);
 258 }
 259