src/tools/pp2shift.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdlib.h>
  40 #include <math.h>
  41 #include "typedefs.h"
  42 #include "futil.h"
  43 #include "macros.h"
  44 #include "physics.h"
  45 #include "smalloc.h"
  46 #include "statutil.h"
  47 #include "gstat.h"
  48 #include "matio.h"
  49 #include "copyrite.h"
  50 #include "gmx_fatal.h"
  51
  52 typedef struct {
  53   int  nx,ny;
  54   real dx,dy;
  55   real **data;
  56 } t_shiftdata;
  57
  58 static real interpolate(real phi,real psi,t_shiftdata *sd)
  59 {
  60   int  iphi,ipsi,iphi1,ipsi1;
  61   real fphi,fpsi,wx0,wx1,wy0,wy1;
  62
  63   /*phi  += M_PI;
  64   if (phi > 2*M_PI) phi -= 2*M_PI;
  65   psi  += M_PI;
  66   if (psi > 2*M_PI) psi -= 2*M_PI;
  67   */
  68   while(phi < 0)
  69     phi += 2*M_PI;
  70   while(psi < 0)
  71     psi += 2*M_PI;
  72   phi    = 2*M_PI-phi;
  73
  74   fphi  = phi*sd->dx;
  75   fpsi  = psi*sd->dy;
  76
  77   iphi  = (int)fphi;
  78   ipsi  = (int)fpsi;
  79   fphi -= iphi; /* Fraction (offset from gridpoint) */
  80   fpsi -= ipsi;
  81
  82   wx0   = 1.0-fphi;
  83   wx1   = fphi;
  84   wy0   = 1.0-fpsi;
  85   wy1   = fpsi;
  86   iphi  = iphi % sd->nx;
  87   ipsi  = ipsi % sd->ny;
  88   iphi1 = (iphi+1) % sd->nx;
  89   ipsi1 = (ipsi+1) % sd->ny;
  90
  91   return (sd->data[iphi]  [ipsi]  * wx0*wy0 +
  92           sd->data[iphi1] [ipsi]  * wx1*wy0 +
  93           sd->data[iphi]  [ipsi1] * wx0*wy1 +
  94           sd->data[iphi1] [ipsi1] * wx1*wy1);
  95 }
  96
  97 static void dump_sd(const char *fn,t_shiftdata *sd)
  98 {
  99   FILE  *fp;
 100   int   i,j;
 101   char  buf[256];
 102   int   nnx,nny,nfac=4,nlevels=20;
 103   real  phi,psi,*x_phi,*y_psi,**newdata;
 104   real  lo,hi;
 105   t_rgb rlo = { 1, 0, 0 }, rhi = { 0, 0, 1 };
 106
 107   nnx = sd->nx*nfac+1;
 108   nny = sd->ny*nfac+1;
 109   snew(x_phi,nnx);
 110   snew(y_psi,nny);
 111   snew(newdata,nnx);
 112   lo =  100000;
 113   hi = -100000;
 114   for(i=0; (i<nnx); i++) {
 115     snew(newdata[i],nny);
 116     phi = i*2*M_PI/(nnx-1);
 117     x_phi[i] =phi*RAD2DEG-180;
 118     for(j=0; (j<nny); j++) {
 119       psi = j*2*M_PI/(nny-1);
 120       if (i == 0)
 121         y_psi[j] = psi*RAD2DEG-180;
 122       /*if (((i % nfac) == 0) && ((j % nfac) == 0))
 123         newdata[i][j] = sd->data[i/nfac][j/nfac];
 124       else*/
 125         newdata[i][j] = interpolate(phi,psi,sd);
 126       lo = min(lo,newdata[i][j]);
 127       hi = max(hi,newdata[i][j]);
 128     }
 129   }
 130   sprintf(buf,"%s.xpm",fn);
 131   fp=ffopen(buf,"w");
 132   write_xpm(fp,0,fn,fn,"Phi","Psi",nnx,nny,
 133             x_phi,y_psi,newdata,lo,hi,rlo,rhi,&nlevels);
 134   for(i=0; (i<nnx); i++)
 135     sfree(newdata[i]);
 136   sfree(newdata);
 137   sfree(x_phi);
 138   sfree(y_psi);
 139 }
 140
 141 static t_shiftdata *read_shifts(const char *fn)
 142 {
 143   FILE        *fp;
 144   double      xx;
 145   int         i,j,nx,ny;
 146   t_shiftdata *sd;
 147
 148   snew(sd,1);
 149   fp=libopen(fn);
 150   if(2 != fscanf(fp,"%d%d",&nx,&ny))
 151   {
 152     gmx_fatal(FARGS,"Error reading from file %s",fn);
 153   }
 154
 155   sd->nx = nx;
 156   sd->ny = ny;
 157   sd->dx = nx/(2*M_PI);
 158   sd->dy = ny/(2*M_PI);
 159   snew(sd->data,nx+1);
 160   for(i=0; (i<=nx); i++) {
 161     snew(sd->data[i],ny+1);
 162     for(j=0; (j<ny); j++) {
 163       if (i == nx)
 164         sd->data[i][j] = sd->data[0][j];
 165       else {
 166         if(1 != fscanf(fp,"%lf",&xx))
 167         {
 168             gmx_fatal(FARGS,"Error reading from file %s",fn);
 169         }
 170         sd->data[i][j] = xx;
 171       }
 172     }
 173     sd->data[i][j] = sd->data[i][0];
 174   }
 175   ffclose(fp);
 176
 177   if (bDebugMode())
 178     dump_sd(fn,sd);
 179
 180   return sd;
 181 }
 182
 183
 184 static void done_shifts(t_shiftdata *sd)
 185 {
 186   int i;
 187
 188   for(i=0; (i<=sd->nx); i++)
 189     sfree(sd->data[i]);
 190   sfree(sd->data);
 191   sfree(sd);
 192 }
 193
 194 void do_pp2shifts(FILE *fp,int nf,int nlist,t_dlist dlist[],real **dih)
 195 {
 196   t_shiftdata *ca_sd,*co_sd,*ha_sd,*cb_sd;
 197   int         i,j,Phi,Psi;
 198   real        phi,psi;
 199   real        ca,co,ha,cb;
 200
 201   /* Read the shift files */
 202   ca_sd = read_shifts("ca-shift.dat");
 203   cb_sd = read_shifts("cb-shift.dat");
 204   ha_sd = read_shifts("ha-shift.dat");
 205   co_sd = read_shifts("co-shift.dat");
 206
 207   fprintf(fp,"\n *** Chemical shifts from the chemical shift index ***\n");
 208   please_cite(fp,"Wishart98a");
 209   fprintf(fp,"%12s  %10s  %10s  %10s  %10s\n",
 210           "Residue","delta Ca","delta Ha","delta CO","delta Cb");
 211   for(i=0; (i<nlist); i++) {
 212     if ((has_dihedral(edPhi,&(dlist[i]))) &&
 213         (has_dihedral(edPsi,&(dlist[i])))) {
 214       Phi = dlist[i].j0[edPhi];
 215       Psi = dlist[i].j0[edPsi];
 216       ca = cb = co = ha = 0;
 217       for(j=0; (j<nf); j++) {
 218         phi = dih[Phi][j];
 219         psi = dih[Psi][j];
 220
 221         ca += interpolate(phi,psi,ca_sd);
 222         cb += interpolate(phi,psi,cb_sd);
 223         co += interpolate(phi,psi,co_sd);
 224         ha += interpolate(phi,psi,ha_sd);
 225       }
 226       fprintf(fp,"%12s  %10g  %10g  %10g  %10g\n",
 227               dlist[i].name,ca/nf,ha/nf,co/nf,cb/nf);
 228     }
 229   }
 230   fprintf(fp,"\n");
 231
 232   /* Free memory */
 233   done_shifts(ca_sd);
 234   done_shifts(cb_sd);
 235   done_shifts(co_sd);
 236   done_shifts(ha_sd);
 237 }
 238