src/tools/hxprops.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 <math.h>
  40 #include <string.h>
  41
  42 #include "typedefs.h"
  43 #include "macros.h"
  44 #include "physics.h"
  45 #include "vec.h"
  46 #include "index.h"
  47 #include "hxprops.h"
  48 #include "smalloc.h"
  49 #include "bondf.h"
  50
  51 int nhelix(int nres,t_bb bb[])
  52 {
  53   int i,n;
  54
  55   for(i=n=0; (i<nres); i++)
  56     if (bb[i].bHelix)
  57       n++;
  58   return n;
  59 }
  60
  61 real ellipticity(int nres,t_bb bb[])
  62 {
  63   typedef struct {
  64     real phi,psi,w;
  65   } t_ppwstr;
  66
  67   static const t_ppwstr ppw[] = {
  68     {  -67,  -44,  0.31 },
  69     {  -66,  -41,  0.31 },
  70     {  -59,  -44,  0.44 },
  71     {  -57,  -47,  0.56 },
  72     {  -53,  -52,  0.78 },
  73     {  -48,  -57,  1.00 },
  74     {  -70.5,-35.8,0.15 },
  75     {  -57,  -79,  0.23 },
  76     {  -38,  -78,  1.20 },
  77     {  -60,  -30,  0.24 },
  78     {  -54,  -28,  0.46 },
  79     {  -44,  -33,  0.68 }
  80   };
  81 #define NPPW asize(ppw)
  82
  83   int   i,j;
  84   const real Xi=1.0;
  85   real  ell,pp2,phi,psi;
  86
  87   ell=0;
  88   for(i=0; (i<nres); i++) {
  89     phi=bb[i].phi;
  90     psi=bb[i].psi;
  91     for(j=0; (j<NPPW); j++) {
  92       pp2=sqr(phi-ppw[j].phi)+sqr(psi-ppw[j].psi);
  93       if (pp2 < 64) {
  94         bb[i].nhx++;
  95         ell+=ppw[j].w;
  96         break;
  97       }
  98     }
  99   }
 100   return ell;
 101 }
 102
 103 real ahx_len(int gnx,atom_id index[],rvec x[],matrix box)
 104 /* Assume we have a list of Calpha atoms only! */
 105 {
 106   rvec dx;
 107
 108   rvec_sub(x[index[0]],x[index[gnx-1]],dx);
 109
 110   return norm(dx);
 111 }
 112
 113 real radius(FILE *fp,int nca,atom_id ca_index[],rvec x[])
 114 /* Assume we have all the backbone */
 115 {
 116   real dl2,dlt;
 117   int  i,ai;
 118
 119   dlt=0;
 120   for(i=0; (i<nca); i++) {
 121     ai=ca_index[i];
 122     dl2=sqr(x[ai][XX])+sqr(x[ai][YY]);
 123
 124     if (fp)
 125       fprintf(fp,"  %10g",dl2);
 126
 127     dlt+=dl2;
 128   }
 129   if (fp)
 130     fprintf(fp,"\n");
 131
 132   return sqrt(dlt/nca);
 133 }
 134
 135 static real rot(rvec x1,rvec x2)
 136 {
 137   real phi1,dphi,cp,sp;
 138   real xx,yy;
 139
 140   phi1=atan2(x1[YY],x1[XX]);
 141   cp=cos(phi1);
 142   sp=sin(phi1);
 143   xx= cp*x2[XX]+sp*x2[YY];
 144   yy=-sp*x2[XX]+cp*x2[YY];
 145
 146   dphi=RAD2DEG*atan2(yy,xx);
 147
 148   return dphi;
 149 }
 150
 151 real twist(FILE *fp,int nca,atom_id caindex[],rvec x[])
 152 {
 153   real pt,dphi;
 154   int  i,a0,a1;
 155
 156   pt=0;
 157   a0=caindex[0];
 158   for(i=1; (i<nca); i++) {
 159     a1=caindex[i];
 160
 161     dphi=rot(x[a0],x[a1]);
 162     if (dphi < -90)
 163       dphi+=360;
 164     pt+=dphi;
 165     a0=a1;
 166   }
 167
 168   return (pt/(nca-1));
 169 }
 170
 171 real ca_phi(int gnx,atom_id index[],rvec x[],matrix box)
 172 /* Assume we have a list of Calpha atoms only! */
 173 {
 174   real phi,phitot;
 175   int  i,ai,aj,ak,al,t1,t2,t3;
 176   rvec r_ij,r_kj,r_kl,m,n;
 177   real sign;
 178
 179   if (gnx <= 4)
 180     return 0;
 181
 182   phitot=0;
 183   for(i=0; (i<gnx-4); i++) {
 184     ai=index[i+0];
 185     aj=index[i+1];
 186     ak=index[i+2];
 187     al=index[i+3];
 188     phi=RAD2DEG*
 189       dih_angle(x[ai],x[aj],x[ak],x[al],NULL,
 190                 r_ij,r_kj,r_kl,m,n,
 191                 &sign,&t1,&t2,&t3);
 192     phitot+=phi;
 193   }
 194
 195   return (phitot/(gnx-4.0));
 196 }
 197
 198 real dip(int nbb,atom_id bbind[],rvec x[],t_atom atom[])
 199 {
 200   int  i,m,ai;
 201   rvec dipje;
 202   real q;
 203
 204   clear_rvec(dipje);
 205   for(i=0; (i<nbb); i++) {
 206     ai=bbind[i];
 207     q=atom[ai].q;
 208     for(m=0; (m<DIM); m++)
 209       dipje[m]+=x[ai][m]*q;
 210   }
 211   return norm(dipje);
 212 }
 213
 214 real rise(int gnx,atom_id index[],rvec x[])
 215 /* Assume we have a list of Calpha atoms only! */
 216 {
 217   real z,z0,ztot;
 218   int  i,ai;
 219
 220   ai=index[0];
 221   z0=x[ai][ZZ];
 222   ztot=0;
 223   for(i=1; (i<gnx); i++) {
 224     ai=index[i];
 225     z=x[ai][ZZ];
 226     ztot+=(z-z0);
 227     z0=z;
 228   }
 229
 230   return (ztot/(gnx-1.0));
 231 }
 232
 233 void av_hblen(FILE *fp3,FILE *fp3a,
 234               FILE *fp4,FILE *fp4a,
 235               FILE *fp5,FILE *fp5a,
 236               real t,int nres,t_bb bb[])
 237 {
 238   int  i,n3=0,n4=0,n5=0;
 239   real d3=0,d4=0,d5=0;
 240
 241   for(i=0; (i<nres-3); i++)
 242     if (bb[i].bHelix) {
 243       fprintf(fp3a,  "%10g",bb[i].d3);
 244       n3++;
 245       d3+=bb[i].d3;
 246       if (i<nres-4) {
 247         fprintf(fp4a,"%10g",bb[i].d4);
 248         n4++;
 249         d4+=bb[i].d4;
 250       }
 251       if (i<nres-5) {
 252         fprintf(fp5a,"%10g",bb[i].d5);
 253         n5++;
 254         d5+=bb[i].d5;
 255       }
 256     }
 257   fprintf(fp3,"%10g  %10g\n",t,d3/n3);
 258   fprintf(fp4,"%10g  %10g\n",t,d4/n4);
 259   fprintf(fp5,"%10g  %10g\n",t,d5/n5);
 260   fprintf(fp3a,"\n");
 261   fprintf(fp4a,"\n");
 262   fprintf(fp5a,"\n");
 263
 264 }
 265
 266 void av_phipsi(FILE *fphi,FILE *fpsi,FILE *fphi2,FILE *fpsi2,
 267                real t,int nres,t_bb bb[])
 268 {
 269   int  i,n=0;
 270   real phi=0,psi=0;
 271
 272   fprintf(fphi2,"%10g",t);
 273   fprintf(fpsi2,"%10g",t);
 274   for(i=0; (i<nres); i++)
 275     if (bb[i].bHelix) {
 276       phi+=bb[i].phi;
 277       psi+=bb[i].psi;
 278       fprintf(fphi2,"  %10g",bb[i].phi);
 279       fprintf(fpsi2,"  %10g",bb[i].psi);
 280       n++;
 281     }
 282   fprintf(fphi,"%10g  %10g\n",t,(phi/n));
 283   fprintf(fpsi,"%10g  %10g\n",t,(psi/n));
 284   fprintf(fphi2,"\n");
 285   fprintf(fpsi2,"\n");
 286 }
 287
 288 static void set_ahcity(int nbb,t_bb bb[])
 289 {
 290   real pp2;
 291   int  n;
 292
 293   for(n=0; (n<nbb); n++) {
 294     pp2=sqr(bb[n].phi-PHI_AHX)+sqr(bb[n].psi-PSI_AHX);
 295
 296     bb[n].bHelix=FALSE;
 297     if (pp2 < 2500) {
 298       if ((bb[n].d4 < 0.36) || ((n > 0) && bb[n-1].bHelix))
 299         bb[n].bHelix=TRUE;
 300     }
 301   }
 302 }
 303
 304 t_bb *mkbbind(const char *fn,int *nres,int *nbb,int res0,
 305               int *nall,atom_id **index,
 306               char ***atomname,t_atom atom[],
 307               t_resinfo *resinfo)
 308 {
 309   static const char * bb_nm[] = { "N", "H", "CA", "C", "O" };
 310 #define NBB asize(bb_nm)
 311   t_bb    *bb;
 312   char    *grpname;
 313   int     ai,i,i0,i1,j,k,ri,rnr,gnx,r0,r1;
 314
 315   fprintf(stderr,"Please select a group containing the entire backbone\n");
 316   rd_index(fn,1,&gnx,index,&grpname);
 317   *nall=gnx;
 318   fprintf(stderr,"Checking group %s\n",grpname);
 319   r0=r1=atom[(*index)[0]].resind;
 320   for(i=1; (i<gnx); i++) {
 321     r0=min(r0,atom[(*index)[i]].resind);
 322     r1=max(r1,atom[(*index)[i]].resind);
 323   }
 324   rnr=r1-r0+1;
 325   fprintf(stderr,"There are %d residues\n",rnr);
 326   snew(bb,rnr);
 327   for(i=0; (i<rnr); i++)
 328     bb[i].N=bb[i].H=bb[i].CA=bb[i].C=bb[i].O=-1,bb[i].resno=res0+i;
 329
 330   for(i=j=0; (i<gnx); i++) {
 331     ai=(*index)[i];
 332     ri=atom[ai].resind-r0;
 333     if (strcmp(*(resinfo[ri].name),"PRO") == 0) {
 334       if (strcmp(*(atomname[ai]),"CD") == 0)
 335         bb[ri].H=ai;
 336     }
 337     for(k=0; (k<NBB); k++)
 338       if (strcmp(bb_nm[k],*(atomname[ai])) == 0) {
 339         j++;
 340         break;
 341       }
 342     switch (k) {
 343     case 0:
 344       bb[ri].N=ai;
 345       break;
 346     case 1:
 347       bb[ri].H=ai;
 348       break;
 349     case 2:
 350       bb[ri].CA=ai;
 351       break;
 352     case 3:
 353       bb[ri].C=ai;
 354       break;
 355     case 4:
 356       bb[ri].O=ai;
 357       break;
 358     default:
 359       break;
 360     }
 361   }
 362
 363   for(i0=0; (i0<rnr); i0++) {
 364     if ((bb[i0].N != -1) && (bb[i0].H != -1) &&
 365         (bb[i0].CA != -1) &&
 366         (bb[i0].C != -1) && (bb[i0].O != -1))
 367       break;
 368   }
 369   for(i1=rnr-1; (i1>=0); i1--) {
 370     if ((bb[i1].N != -1) && (bb[i1].H != -1) &&
 371         (bb[i1].CA != -1) &&
 372         (bb[i1].C != -1) && (bb[i1].O != -1))
 373       break;
 374   }
 375   if (i0 == 0)
 376     i0++;
 377   if (i1 == rnr-1)
 378     i1--;
 379
 380   for(i=i0; (i<i1); i++) {
 381     bb[i].Cprev=bb[i-1].C;
 382     bb[i].Nnext=bb[i+1].N;
 383   }
 384   rnr=max(0,i1-i0+1);
 385   fprintf(stderr,"There are %d complete backbone residues (from %d to %d)\n",
 386           rnr,bb[i0].resno,bb[i1].resno);
 387   if (rnr==0)
 388     gmx_fatal(FARGS,"rnr==0");
 389   for(i=0; (i<rnr); i++,i0++)
 390     bb[i]=bb[i0];
 391
 392   /* Set the labels */
 393   for(i=0; (i<rnr); i++) {
 394     ri=atom[bb[i].CA].resind;
 395     sprintf(bb[i].label,"%s%d",*(resinfo[ri].name),ri+res0);
 396   }
 397
 398   *nres=rnr;
 399   *nbb=rnr*asize(bb_nm);
 400
 401   return bb;
 402 }
 403
 404 real pprms(FILE *fp,int nbb,t_bb bb[])
 405 {
 406   int  i,n;
 407   real rms,rmst,rms2;
 408
 409   rmst=rms2=0;
 410   for(i=n=0; (i<nbb); i++) {
 411     if (bb[i].bHelix) {
 412       rms=sqrt(bb[i].pprms2);
 413       rmst+=rms;
 414       rms2+=bb[i].pprms2;
 415       fprintf(fp,"%10g  ",rms);
 416       n++;
 417     }
 418   }
 419   fprintf(fp,"\n");
 420   rms=sqrt(rms2/n-sqr(rmst/n));
 421
 422   return rms;
 423 }
 424
 425 void calc_hxprops(int nres,t_bb bb[],rvec x[],matrix box)
 426 {
 427   int  i,ao,an,t1,t2,t3;
 428   rvec dx,r_ij,r_kj,r_kl,m,n;
 429   real sign;
 430
 431   for(i=0; (i<nres); i++) {
 432     ao=bb[i].O;
 433     bb[i].d4=bb[i].d3=bb[i].d5=0;
 434     if (i < nres-3) {
 435       an=bb[i+3].N;
 436       rvec_sub(x[ao],x[an],dx);
 437       bb[i].d3=norm(dx);
 438     }
 439     if (i < nres-4) {
 440       an=bb[i+4].N;
 441       rvec_sub(x[ao],x[an],dx);
 442       bb[i].d4=norm(dx);
 443     }
 444     if (i < nres-5) {
 445       an=bb[i+5].N;
 446       rvec_sub(x[ao],x[an],dx);
 447       bb[i].d5=norm(dx);
 448     }
 449
 450     bb[i].phi=RAD2DEG*
 451       dih_angle(x[bb[i].Cprev],x[bb[i].N],x[bb[i].CA],x[bb[i].C],NULL,
 452                 r_ij,r_kj,r_kl,m,n,
 453                 &sign,&t1,&t2,&t3);
 454     bb[i].psi=RAD2DEG*
 455       dih_angle(x[bb[i].N],x[bb[i].CA],x[bb[i].C],x[bb[i].Nnext],NULL,
 456                 r_ij,r_kj,r_kl,m,n,
 457                 &sign,&t1,&t2,&t3);
 458     bb[i].pprms2=sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX);
 459
 460     bb[i].jcaha+=
 461       1.4*sin((bb[i].psi+138.0)*DEG2RAD) -
 462         4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
 463           2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0));
 464   }
 465 }
 466
 467 static void check_ahx(int nres,t_bb bb[],rvec x[],
 468                       int *hstart,int *hend)
 469 {
 470   int  h0,h1,h0sav,h1sav;
 471
 472   set_ahcity(nres,bb);
 473   h0=h0sav=h1sav=0;
 474   do {
 475     for(; (!bb[h0].bHelix) && (h0<nres-4); h0++)
 476       ;
 477     for(h1=h0; bb[h1+1].bHelix && (h1<nres-1); h1++)
 478       ;
 479     if (h1 > h0) {
 480       /*fprintf(stderr,"Helix from %d to %d\n",h0,h1);*/
 481       if (h1-h0 > h1sav-h0sav) {
 482         h0sav=h0;
 483         h1sav=h1;
 484       }
 485     }
 486     h0=h1+1;
 487   } while (h1 < nres-1);
 488   *hstart=h0sav;
 489   *hend=h1sav;
 490 }
 491
 492 void do_start_end(int nres,t_bb bb[],rvec x[],int *nbb,atom_id bbindex[],
 493                   int *nca,atom_id caindex[],
 494                   gmx_bool bRange,int rStart,int rEnd)
 495 {
 496   int    i,j,hstart=0,hend=0;
 497
 498   if (bRange) {
 499     for(i=0; (i<nres); i++) {
 500       if ((bb[i].resno >= rStart) && (bb[i].resno <= rEnd))
 501         bb[i].bHelix=TRUE;
 502       if (bb[i].resno == rStart)
 503         hstart=i;
 504       if (bb[i].resno == rEnd)
 505         hend=i;
 506     }
 507   }
 508   else {
 509     /* Find start and end of longest helix fragment */
 510     check_ahx(nres,bb,x,&hstart,&hend);
 511   }
 512   fprintf(stderr,"helix from: %d through %d\n",
 513           bb[hstart].resno,bb[hend].resno);
 514
 515   for(j=0,i=hstart; (i<=hend); i++) {
 516     bbindex[j++]=bb[i].N;
 517     bbindex[j++]=bb[i].H;
 518     bbindex[j++]=bb[i].CA;
 519     bbindex[j++]=bb[i].C;
 520     bbindex[j++]=bb[i].O;
 521     caindex[i-hstart]=bb[i].CA;
 522   }
 523   *nbb=j;
 524   *nca=(hend-hstart+1);
 525 }
 526
 527 void pr_bb(FILE *fp,int nres,t_bb bb[])
 528 {
 529   int i;
 530
 531   fprintf(fp,"\n");
 532   fprintf(fp,"%3s %3s %3s %3s %3s %7s %7s %7s %7s %7s %3s\n",
 533           "AA","N","Ca","C","O","Phi","Psi","D3","D4","D5","Hx?");
 534   for(i=0; (i<nres); i++) {
 535     fprintf(fp,"%3d %3d %3d %3d %3d %7.2f %7.2f %7.3f %7.3f %7.3f %3s\n",
 536             bb[i].resno,bb[i].N,bb[i].CA,bb[i].C,bb[i].O,
 537             bb[i].phi,bb[i].psi,bb[i].d3,bb[i].d4,bb[i].d5,
 538             bb[i].bHelix ? "Yes" : "No");
 539   }
 540   fprintf(fp,"\n");
 541 }
 542