src/tools/hxprops.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_hxprops_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "macros.h"
  34 #include "physics.h"
  35 #include "vec.h"
  36 #include "rdgroup.h"
  37 #include "hxprops.h"
  38 #include "smalloc.h"
  39 #include "bondf.h"
  40
  41 int nhelix(int nres,t_bb bb[])
  42 {
  43   int i,n;
  44
  45   for(i=n=0; (i<nres); i++)
  46     if (bb[i].bHelix)
  47       n++;
  48   return n;
  49 }
  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+=Xi;
  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;
 166   rvec r_ij,r_kj,r_kl,m,n;
 167   real cos_phi,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(box,
 180                 x[ai],x[aj],x[ak],x[al],
 181                 r_ij,r_kj,r_kl,m,n,
 182                 &cos_phi,&sign);
 183     phitot+=phi;
 184   }
 185
 186   return (phitot/(gnx-4.0));
 187 }
 188
 189 real dip(int nbb,atom_id bbind[],rvec x[],t_atom atom[])
 190 {
 191   int  i,m,ai;
 192   rvec dipje;
 193   real q;
 194
 195   clear_rvec(dipje);
 196   for(i=0; (i<nbb); i++) {
 197     ai=bbind[i];
 198     q=atom[ai].q;
 199     for(m=0; (m<DIM); m++)
 200       dipje[m]+=x[ai][m]*q;
 201   }
 202   return norm(dipje);
 203 }
 204
 205 real rise(int gnx,atom_id index[],rvec x[])
 206 /* Assume we have a list of Calpha atoms only! */
 207 {
 208   real z,z0,ztot;
 209   int  i,ai;
 210
 211   ai=index[0];
 212   z0=x[ai][ZZ];
 213   ztot=0;
 214   for(i=1; (i<gnx); i++) {
 215     ai=index[i];
 216     z=x[ai][ZZ];
 217     ztot+=(z-z0);
 218     z0=z;
 219   }
 220
 221   return (ztot/(gnx-1.0));
 222 }
 223
 224 void av_hblen(FILE *fp3,FILE *fp3a,
 225               FILE *fp4,FILE *fp4a,
 226               FILE *fp5,FILE *fp5a,
 227               real t,int nres,t_bb bb[])
 228 {
 229   int  i,n3=0,n4=0,n5=0;
 230   real d3=0,d4=0,d5=0;
 231
 232   for(i=0; (i<nres-3); i++)
 233     if (bb[i].bHelix) {
 234       fprintf(fp3a,  "%10g",bb[i].d3);
 235       n3++;
 236       d3+=bb[i].d3;
 237       if (i<nres-4) {
 238         fprintf(fp4a,"%10g",bb[i].d4);
 239         n4++;
 240         d4+=bb[i].d4;
 241       }
 242       if (i<nres-5) {
 243         fprintf(fp5a,"%10g",bb[i].d5);
 244         n5++;
 245         d5+=bb[i].d5;
 246       }
 247     }
 248   fprintf(fp3,"%10g  %10g\n",t,d3/n3);
 249   fprintf(fp4,"%10g  %10g\n",t,d4/n4);
 250   fprintf(fp5,"%10g  %10g\n",t,d5/n5);
 251   fprintf(fp3a,"\n");
 252   fprintf(fp4a,"\n");
 253   fprintf(fp5a,"\n");
 254
 255 }
 256
 257 void av_phipsi(FILE *fphi,FILE *fpsi,FILE *fphi2,FILE *fpsi2,
 258                real t,int nres,t_bb bb[])
 259 {
 260   int  i,n=0;
 261   real phi=0,psi=0;
 262
 263   fprintf(fphi2,"%10g",t);
 264   fprintf(fpsi2,"%10g",t);
 265   for(i=0; (i<nres); i++)
 266     if (bb[i].bHelix) {
 267       phi+=bb[i].phi;
 268       psi+=bb[i].psi;
 269       fprintf(fphi2,"  %10g",bb[i].phi);
 270       fprintf(fpsi2,"  %10g",bb[i].psi);
 271       n++;
 272     }
 273   fprintf(fphi,"%10g  %10g\n",t,(phi/n));
 274   fprintf(fpsi,"%10g  %10g\n",t,(psi/n));
 275   fprintf(fphi2,"\n");
 276   fprintf(fpsi2,"\n");
 277 }
 278
 279 static void set_ahcity(int nbb,t_bb bb[])
 280 {
 281   real pp2;
 282   int  n;
 283
 284   for(n=0; (n<nbb); n++) {
 285     pp2=sqr(bb[n].phi-PHI_AHX)+sqr(bb[n].psi-PSI_AHX);
 286
 287     bb[n].bHelix=FALSE;
 288     if (pp2 < 2500) {
 289       if ((bb[n].d4 < 0.36) || ((n > 0) && bb[n-1].bHelix))
 290         bb[n].bHelix=TRUE;
 291     }
 292   }
 293 }
 294
 295 t_bb *mkbbind(char *fn,int *nres,int *nbb,int res0,
 296               int *nall,atom_id **index,
 297               char ***atomname,t_atom atom[],
 298               char ***resname)
 299 {
 300   static char * bb_nm[] = { "N", "H", "CA", "C", "O" };
 301 #define NBB asize(bb_nm)
 302   t_bb    *bb;
 303   char    *grpname;
 304   int     ai,i,i0,i1,j,k,ri,rnr,gnx,r0,r1;
 305
 306   fprintf(stderr,"Please select a group containing the entire backbone\n");
 307   rd_index(fn,1,&gnx,index,&grpname);
 308   *nall=gnx;
 309   fprintf(stderr,"Checking group %s\n",grpname);
 310   r0=r1=atom[(*index)[0]].resnr;
 311   for(i=1; (i<gnx); i++) {
 312     r0=min(r0,atom[(*index)[i]].resnr);
 313     r1=max(r1,atom[(*index)[i]].resnr);
 314   }
 315   rnr=r1-r0+1;
 316   fprintf(stderr,"There are %d residues\n",rnr);
 317   snew(bb,rnr);
 318   for(i=0; (i<rnr); i++)
 319     bb[i].N=bb[i].H=bb[i].CA=bb[i].C=bb[i].O=-1,bb[i].resno=res0+i;
 320
 321   for(i=j=0; (i<gnx); i++) {
 322     ai=(*index)[i];
 323     ri=atom[ai].resnr-r0;
 324     if (strcmp(*(resname[ri]),"PRO") == 0) {
 325       if (strcmp(*(atomname[ai]),"CD") == 0)
 326         bb[ri].H=ai;
 327     }
 328     for(k=0; (k<NBB); k++)
 329       if (strcmp(bb_nm[k],*(atomname[ai])) == 0) {
 330         j++;
 331         break;
 332       }
 333     switch (k) {
 334     case 0:
 335       bb[ri].N=ai;
 336       break;
 337     case 1:
 338       bb[ri].H=ai;
 339       break;
 340     case 2:
 341       bb[ri].CA=ai;
 342       break;
 343     case 3:
 344       bb[ri].C=ai;
 345       break;
 346     case 4:
 347       bb[ri].O=ai;
 348       break;
 349     default:
 350       break;
 351     }
 352   }
 353
 354   for(i0=0; (i0<rnr); i0++) {
 355     if ((bb[i0].N != -1) && (bb[i0].H != -1) &&
 356         (bb[i0].CA != -1) &&
 357         (bb[i0].C != -1) && (bb[i0].O != -1))
 358       break;
 359   }
 360   for(i1=rnr-1; (i1>=0); i1--) {
 361     if ((bb[i1].N != -1) && (bb[i1].H != -1) &&
 362         (bb[i1].CA != -1) &&
 363         (bb[i1].C != -1) && (bb[i1].O != -1))
 364       break;
 365   }
 366   if (i0 == 0)
 367     i0++;
 368   if (i1 == rnr-1)
 369     i1--;
 370
 371   for(i=i0; (i<i1); i++) {
 372     bb[i].Cprev=bb[i-1].C;
 373     bb[i].Nnext=bb[i+1].N;
 374   }
 375   rnr=max(0,i1-i0+1);
 376   fprintf(stderr,"There are %d complete backbone residues (from %d to %d)\n",
 377           rnr,bb[i0].resno,bb[i1].resno);
 378   if (rnr==0)
 379     fatal_error(0,"rnr==0");
 380   for(i=0; (i<rnr); i++,i0++)
 381     bb[i]=bb[i0];
 382
 383   /* Set the labels */
 384   for(i=0; (i<rnr); i++) {
 385     ri=atom[bb[i].CA].resnr;
 386     sprintf(bb[i].label,"%s%d",*(resname[ri]),ri+res0);
 387   }
 388
 389   *nres=rnr;
 390   *nbb=rnr*asize(bb_nm);
 391
 392   return bb;
 393 }
 394
 395 real pprms(FILE *fp,int nbb,t_bb bb[])
 396 {
 397   int  i,n;
 398   real rms,rmst,rms2;
 399
 400   rmst=rms2=0;
 401   for(i=n=0; (i<nbb); i++) {
 402     if (bb[i].bHelix) {
 403       rms=sqrt(bb[i].pprms2);
 404       rmst+=rms;
 405       rms2+=bb[i].pprms2;
 406       fprintf(fp,"%10g  ",rms);
 407       n++;
 408     }
 409   }
 410   fprintf(fp,"\n");
 411   rms=sqrt(rms2/n-sqr(rmst/n));
 412
 413   return rms;
 414 }
 415
 416 void calc_hxprops(int nres,t_bb bb[],rvec x[],matrix box)
 417 {
 418   int  i,ao,an;
 419   rvec dx,r_ij,r_kj,r_kl,m,n;
 420   real cos_phi,sign;
 421
 422   for(i=0; (i<nres); i++) {
 423     ao=bb[i].O;
 424     bb[i].d4=bb[i].d3=bb[i].d5=0;
 425     if (i < nres-3) {
 426       an=bb[i+3].N;
 427       rvec_sub(x[ao],x[an],dx);
 428       bb[i].d3=norm(dx);
 429     }
 430     if (i < nres-4) {
 431       an=bb[i+4].N;
 432       rvec_sub(x[ao],x[an],dx);
 433       bb[i].d4=norm(dx);
 434     }
 435     if (i < nres-5) {
 436       an=bb[i+5].N;
 437       rvec_sub(x[ao],x[an],dx);
 438       bb[i].d5=norm(dx);
 439     }
 440
 441     bb[i].phi=RAD2DEG*
 442       dih_angle(box,
 443                 x[bb[i].Cprev],x[bb[i].N],x[bb[i].CA],x[bb[i].C],
 444                 r_ij,r_kj,r_kl,m,n,
 445                 &cos_phi,&sign);
 446     bb[i].psi=RAD2DEG*
 447       dih_angle(box,
 448                 x[bb[i].N],x[bb[i].CA],x[bb[i].C],x[bb[i].Nnext],
 449                 r_ij,r_kj,r_kl,m,n,
 450                 &cos_phi,&sign);
 451     bb[i].pprms2=sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX);
 452
 453     bb[i].jcaha+=
 454       1.4*sin((bb[i].psi+138.0)*DEG2RAD) -
 455         4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
 456           2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0));
 457   }
 458 }
 459
 460 static void check_ahx(int nres,t_bb bb[],rvec x[],
 461                       int *hstart,int *hend)
 462 {
 463   int  h0,h1,h0sav,h1sav;
 464
 465   set_ahcity(nres,bb);
 466   h0=h0sav=h1sav=0;
 467   do {
 468     for(; (!bb[h0].bHelix) && (h0<nres-4); h0++)
 469       ;
 470     for(h1=h0; bb[h1+1].bHelix && (h1<nres-1); h1++)
 471       ;
 472     if (h1 > h0) {
 473       /*fprintf(stderr,"Helix from %d to %d\n",h0,h1);*/
 474       if (h1-h0 > h1sav-h0sav) {
 475         h0sav=h0;
 476         h1sav=h1;
 477       }
 478     }
 479     h0=h1+1;
 480   } while (h1 < nres-1);
 481   *hstart=h0sav;
 482   *hend=h1sav;
 483 }
 484
 485 void do_start_end(int nres,t_bb bb[],rvec x[],int *nbb,atom_id bbindex[],
 486                   int *nca,atom_id caindex[],
 487                   bool bRange,int rStart,int rEnd)
 488 {
 489   int    i,j,hstart,hend;
 490
 491   if (bRange) {
 492     for(i=0; (i<nres); i++) {
 493       if ((bb[i].resno >= rStart) && (bb[i].resno <= rEnd))
 494         bb[i].bHelix=TRUE;
 495       if (bb[i].resno == rStart)
 496         hstart=i;
 497       if (bb[i].resno == rEnd)
 498         hend=i;
 499     }
 500   }
 501   else {
 502     /* Find start and end of longest helix fragment */
 503     check_ahx(nres,bb,x,&hstart,&hend);
 504   }
 505   fprintf(stderr,"helix from: %d thru %d\n",
 506           bb[hstart].resno,bb[hend].resno);
 507
 508   for(j=0,i=hstart; (i<=hend); i++) {
 509     bbindex[j++]=bb[i].N;
 510     bbindex[j++]=bb[i].H;
 511     bbindex[j++]=bb[i].CA;
 512     bbindex[j++]=bb[i].C;
 513     bbindex[j++]=bb[i].O;
 514     caindex[i-hstart]=bb[i].CA;
 515   }
 516   *nbb=j;
 517   *nca=(hend-hstart+1);
 518 }
 519
 520 void pr_bb(FILE *fp,int nres,t_bb bb[])
 521 {
 522   int i;
 523
 524   fprintf(fp,"\n");
 525   fprintf(fp,"%3s %3s %3s %3s %3s %7s %7s %7s %7s %7s %3s\n",
 526           "AA","N","Ca","C","O","Phi","Psi","D3","D4","D5","Hx?");
 527   for(i=0; (i<nres); i++) {
 528     fprintf(fp,"%3d %3d %3d %3d %3d %7.2f %7.2f %7.3f %7.3f %7.3f %3s\n",
 529             bb[i].resno,bb[i].N,bb[i].CA,bb[i].C,bb[i].O,
 530             bb[i].phi,bb[i].psi,bb[i].d3,bb[i].d4,bb[i].d5,
 531             bb[i].bHelix ? "Yes" : "No");
 532   }
 533   fprintf(fp,"\n");
 534 }
 535