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