src/tools/pp2shift.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_pp2shift_c = "$Id$";
  30
  31 #include <stdlib.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "futil.h"
  35 #include "macros.h"
  36 #include "physics.h"
  37 #include "smalloc.h"
  38 #include "statutil.h"
  39 #include "pp2shift.h"
  40 #include "matio.h"
  41 #include "copyrite.h"
  42
  43 typedef struct {
  44   int  nx,ny;
  45   real dx,dy;
  46   real **data;
  47 } t_shiftdata;
  48
  49 static real interpolate(real phi,real psi,t_shiftdata *sd)
  50 {
  51   int  iphi,ipsi,iphi1,ipsi1;
  52   real fphi,fpsi,wx0,wx1,wy0,wy1;
  53
  54   /*phi  += M_PI;
  55   if (phi > 2*M_PI) phi -= 2*M_PI;
  56   psi  += M_PI;
  57   if (psi > 2*M_PI) psi -= 2*M_PI;
  58   */
  59   while(phi < 0)
  60     phi += 2*M_PI;
  61   while(psi < 0)
  62     psi += 2*M_PI;
  63   phi    = 2*M_PI-phi;
  64
  65   fphi  = phi*sd->dx;
  66   fpsi  = psi*sd->dy;
  67
  68   iphi  = (int)fphi;
  69   ipsi  = (int)fpsi;
  70   fphi -= iphi; /* Fraction (offset from gridpoint) */
  71   fpsi -= ipsi;
  72
  73   wx0   = 1.0-fphi;
  74   wx1   = fphi;
  75   wy0   = 1.0-fpsi;
  76   wy1   = fpsi;
  77   iphi  = iphi % sd->nx;
  78   ipsi  = ipsi % sd->ny;
  79   iphi1 = (iphi+1) % sd->nx;
  80   ipsi1 = (ipsi+1) % sd->ny;
  81
  82   return (sd->data[iphi]  [ipsi]  * wx0*wy0 +
  83           sd->data[iphi1] [ipsi]  * wx1*wy0 +
  84           sd->data[iphi]  [ipsi1] * wx0*wy1 +
  85           sd->data[iphi1] [ipsi1] * wx1*wy1);
  86 }
  87
  88 static void dump_sd(char *fn,t_shiftdata *sd)
  89 {
  90   FILE  *fp;
  91   int   i,j;
  92   char  buf[256];
  93   int   nnx,nny,nfac=4,nlevels=20;
  94   real  phi,psi,*x_phi,*y_psi,**newdata;
  95   real  lo,hi;
  96   t_rgb rlo = { 1, 0, 0 }, rhi = { 0, 0, 1 };
  97
  98   nnx = sd->nx*nfac+1;
  99   nny = sd->ny*nfac+1;
 100   snew(x_phi,nnx);
 101   snew(y_psi,nny);
 102   snew(newdata,nnx);
 103   lo =  100000;
 104   hi = -100000;
 105   for(i=0; (i<nnx); i++) {
 106     snew(newdata[i],nny);
 107     phi = i*2*M_PI/(nnx-1);
 108     x_phi[i] =phi*RAD2DEG-180;
 109     for(j=0; (j<nny); j++) {
 110       psi = j*2*M_PI/(nny-1);
 111       if (i == 0)
 112         y_psi[j] = psi*RAD2DEG-180;
 113       /*if (((i % nfac) == 0) && ((j % nfac) == 0))
 114         newdata[i][j] = sd->data[i/nfac][j/nfac];
 115       else*/
 116         newdata[i][j] = interpolate(phi,psi,sd);
 117       lo = min(lo,newdata[i][j]);
 118       hi = max(hi,newdata[i][j]);
 119     }
 120   }
 121   sprintf(buf,"%s.xpm",fn);
 122   fp=ffopen(buf,"w");
 123   write_xpm(fp,fn,fn,"Phi","Psi",nnx,nny,
 124             x_phi,y_psi,newdata,lo,hi,rlo,rhi,&nlevels);
 125   for(i=0; (i<nnx); i++)
 126     sfree(newdata[i]);
 127   sfree(newdata);
 128   sfree(x_phi);
 129   sfree(y_psi);
 130 }
 131
 132 static t_shiftdata *read_shifts(char *fn)
 133 {
 134   FILE        *fp;
 135   double      xx;
 136   int         i,j,nx,ny;
 137   t_shiftdata *sd;
 138
 139   snew(sd,1);
 140   fp=libopen(fn);
 141   fscanf(fp,"%d%d",&nx,&ny);
 142   sd->nx = nx;
 143   sd->ny = ny;
 144   sd->dx = nx/(2*M_PI);
 145   sd->dy = ny/(2*M_PI);
 146   snew(sd->data,nx+1);
 147   for(i=0; (i<=nx); i++) {
 148     snew(sd->data[i],ny+1);
 149     for(j=0; (j<ny); j++) {
 150       if (i == nx)
 151         sd->data[i][j] = sd->data[0][j];
 152       else {
 153         fscanf(fp,"%lf",&xx);
 154         sd->data[i][j] = xx;
 155       }
 156     }
 157     sd->data[i][j] = sd->data[i][0];
 158   }
 159   fclose(fp);
 160
 161   if (bDebugMode())
 162     dump_sd(fn,sd);
 163
 164   return sd;
 165 }
 166
 167
 168 static void done_shifts(t_shiftdata *sd)
 169 {
 170   int i;
 171
 172   for(i=0; (i<=sd->nx); i++)
 173     sfree(sd->data[i]);
 174   sfree(sd->data);
 175   sfree(sd);
 176 }
 177
 178 void do_pp2shifts(FILE *fp,int nf,int nlist,t_dlist dlist[],real **dih)
 179 {
 180   t_shiftdata *ca_sd,*co_sd,*ha_sd,*cb_sd;
 181   int         i,j,Phi,Psi;
 182   real        phi,psi;
 183   real        ca,co,ha,cb;
 184
 185   /* Read the shift files */
 186   ca_sd = read_shifts("ca-shift.dat");
 187   cb_sd = read_shifts("cb-shift.dat");
 188   ha_sd = read_shifts("ha-shift.dat");
 189   co_sd = read_shifts("co-shift.dat");
 190
 191   fprintf(fp,"\n *** Chemical shifts from the chemical shift index ***\n");
 192   please_cite(fp,"Wishart98a");
 193   fprintf(fp,"%12s  %10s  %10s  %10s  %10s\n",
 194           "Residue","delta Ca","delta Ha","delta CO","delta Cb");
 195   for(i=0; (i<nlist); i++) {
 196     if ((has_dihedral(edPhi,&(dlist[i]))) &&
 197         (has_dihedral(edPsi,&(dlist[i])))) {
 198       Phi = dlist[i].j0[edPhi];
 199       Psi = dlist[i].j0[edPsi];
 200       ca = cb = co = ha = 0;
 201       for(j=0; (j<nf); j++) {
 202         phi = dih[Phi][j];
 203         psi = dih[Psi][j];
 204
 205         ca += interpolate(phi,psi,ca_sd);
 206         cb += interpolate(phi,psi,cb_sd);
 207         co += interpolate(phi,psi,co_sd);
 208         ha += interpolate(phi,psi,ha_sd);
 209       }
 210       fprintf(fp,"%12s  %10g  %10g  %10g  %10g\n",
 211               dlist[i].name,ca/nf,ha/nf,co/nf,cb/nf);
 212     }
 213   }
 214   fprintf(fp,"\n");
 215
 216   /* Free memory */
 217   done_shifts(ca_sd);
 218   done_shifts(cb_sd);
 219   done_shifts(co_sd);
 220   done_shifts(ha_sd);
 221 }
 222