src/tools/hxprops.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <string.h>
  44
  45 #include "typedefs.h"
  46 #include "macros.h"
  47 #include "physics.h"
  48 #include "vec.h"
  49 #include "index.h"
  50 #include "hxprops.h"
  51 #include "smalloc.h"
  52 #include "bondf.h"
  53
  54 real ellipticity(int nres,t_bb bb[])
  55 {
  56   typedef struct {
  57     real phi,psi,w;
  58   } t_ppwstr;
  59
  60   static const t_ppwstr ppw[] = {
  61     {  -67,  -44,  0.31 },
  62     {  -66,  -41,  0.31 },
  63     {  -59,  -44,  0.44 },
  64     {  -57,  -47,  0.56 },
  65     {  -53,  -52,  0.78 },
  66     {  -48,  -57,  1.00 },
  67     {  -70.5,-35.8,0.15 },
  68     {  -57,  -79,  0.23 },
  69     {  -38,  -78,  1.20 },
  70     {  -60,  -30,  0.24 },
  71     {  -54,  -28,  0.46 },
  72     {  -44,  -33,  0.68 }
  73   };
  74 #define NPPW asize(ppw)
  75
  76   int   i,j;
  77   const real Xi=1.0;
  78   real  ell,pp2,phi,psi;
  79
  80   ell=0;
  81   for(i=0; (i<nres); i++) {
  82     phi=bb[i].phi;
  83     psi=bb[i].psi;
  84     for(j=0; (j<NPPW); j++) {
  85       pp2=sqr(phi-ppw[j].phi)+sqr(psi-ppw[j].psi);
  86       if (pp2 < 64) {
  87         bb[i].nhx++;
  88         ell+=ppw[j].w;
  89         break;
  90       }
  91     }
  92   }
  93   return ell;
  94 }
  95
  96 real ahx_len(int gnx,atom_id index[],rvec x[],matrix box)
  97 /* Assume we have a list of Calpha atoms only! */
  98 {
  99   rvec dx;
 100
 101   rvec_sub(x[index[0]],x[index[gnx-1]],dx);
 102
 103   return norm(dx);
 104 }
 105
 106 real radius(FILE *fp,int nca,atom_id ca_index[],rvec x[])
 107 /* Assume we have all the backbone */
 108 {
 109   real dl2,dlt;
 110   int  i,ai;
 111
 112   dlt=0;
 113   for(i=0; (i<nca); i++) {
 114     ai=ca_index[i];
 115     dl2=sqr(x[ai][XX])+sqr(x[ai][YY]);
 116
 117     if (fp)
 118       fprintf(fp,"  %10g",dl2);
 119
 120     dlt+=dl2;
 121   }
 122   if (fp)
 123     fprintf(fp,"\n");
 124
 125   return sqrt(dlt/nca);
 126 }
 127
 128 static real rot(rvec x1,rvec x2)
 129 {
 130   real phi1,dphi,cp,sp;
 131   real xx,yy;
 132
 133   phi1=atan2(x1[YY],x1[XX]);
 134   cp=cos(phi1);
 135   sp=sin(phi1);
 136   xx= cp*x2[XX]+sp*x2[YY];
 137   yy=-sp*x2[XX]+cp*x2[YY];
 138
 139   dphi=RAD2DEG*atan2(yy,xx);
 140
 141   return dphi;
 142 }
 143
 144 real twist(FILE *fp,int nca,atom_id caindex[],rvec x[])
 145 {
 146   real pt,dphi;
 147   int  i,a0,a1;
 148
 149   pt=0;
 150   a0=caindex[0];
 151   for(i=1; (i<nca); i++) {
 152     a1=caindex[i];
 153
 154     dphi=rot(x[a0],x[a1]);
 155     if (dphi < -90)
 156       dphi+=360;
 157     pt+=dphi;
 158     a0=a1;
 159   }
 160
 161   return (pt/(nca-1));
 162 }
 163
 164 real ca_phi(int gnx,atom_id index[],rvec x[],matrix box)
 165 /* Assume we have a list of Calpha atoms only! */
 166 {
 167   real phi,phitot;
 168   int  i,ai,aj,ak,al,t1,t2,t3;
 169   rvec r_ij,r_kj,r_kl,m,n;
 170   real sign;
 171
 172   if (gnx <= 4)
 173     return 0;
 174
 175   phitot=0;
 176   for(i=0; (i<gnx-4); i++) {
 177     ai=index[i+0];
 178     aj=index[i+1];
 179     ak=index[i+2];
 180     al=index[i+3];
 181     phi=RAD2DEG*
 182       dih_angle(x[ai],x[aj],x[ak],x[al],NULL,
 183                 r_ij,r_kj,r_kl,m,n,
 184                 &sign,&t1,&t2,&t3);
 185     phitot+=phi;
 186   }
 187
 188   return (phitot/(gnx-4.0));
 189 }
 190
 191 real dip(int nbb,atom_id bbind[],rvec x[],t_atom atom[])
 192 {
 193   int  i,m,ai;
 194   rvec dipje;
 195   real q;
 196
 197   clear_rvec(dipje);
 198   for(i=0; (i<nbb); i++) {
 199     ai=bbind[i];
 200     q=atom[ai].q;
 201     for(m=0; (m<DIM); m++)
 202       dipje[m]+=x[ai][m]*q;
 203   }
 204   return norm(dipje);
 205 }
 206
 207 real rise(int gnx,atom_id index[],rvec x[])
 208 /* Assume we have a list of Calpha atoms only! */
 209 {
 210   real z,z0,ztot;
 211   int  i,ai;
 212
 213   ai=index[0];
 214   z0=x[ai][ZZ];
 215   ztot=0;
 216   for(i=1; (i<gnx); i++) {
 217     ai=index[i];
 218     z=x[ai][ZZ];
 219     ztot+=(z-z0);
 220     z0=z;
 221   }
 222
 223   return (ztot/(gnx-1.0));
 224 }
 225
 226 void av_hblen(FILE *fp3,FILE *fp3a,
 227               FILE *fp4,FILE *fp4a,
 228               FILE *fp5,FILE *fp5a,
 229               real t,int nres,t_bb bb[])
 230 {
 231   int  i,n3=0,n4=0,n5=0;
 232   real d3=0,d4=0,d5=0;
 233
 234   for(i=0; (i<nres-3); i++)
 235     if (bb[i].bHelix) {
 236       fprintf(fp3a,  "%10g",bb[i].d3);
 237       n3++;
 238       d3+=bb[i].d3;
 239       if (i<nres-4) {
 240         fprintf(fp4a,"%10g",bb[i].d4);
 241         n4++;
 242         d4+=bb[i].d4;
 243       }
 244       if (i<nres-5) {
 245         fprintf(fp5a,"%10g",bb[i].d5);
 246         n5++;
 247         d5+=bb[i].d5;
 248       }
 249     }
 250   fprintf(fp3,"%10g  %10g\n",t,d3/n3);
 251   fprintf(fp4,"%10g  %10g\n",t,d4/n4);
 252   fprintf(fp5,"%10g  %10g\n",t,d5/n5);
 253   fprintf(fp3a,"\n");
 254   fprintf(fp4a,"\n");
 255   fprintf(fp5a,"\n");
 256
 257 }
 258
 259 void av_phipsi(FILE *fphi,FILE *fpsi,FILE *fphi2,FILE *fpsi2,
 260                real t,int nres,t_bb bb[])
 261 {
 262   int  i,n=0;
 263   real phi=0,psi=0;
 264
 265   fprintf(fphi2,"%10g",t);
 266   fprintf(fpsi2,"%10g",t);
 267   for(i=0; (i<nres); i++)
 268     if (bb[i].bHelix) {
 269       phi+=bb[i].phi;
 270       psi+=bb[i].psi;
 271       fprintf(fphi2,"  %10g",bb[i].phi);
 272       fprintf(fpsi2,"  %10g",bb[i].psi);
 273       n++;
 274     }
 275   fprintf(fphi,"%10g  %10g\n",t,(phi/n));
 276   fprintf(fpsi,"%10g  %10g\n",t,(psi/n));
 277   fprintf(fphi2,"\n");
 278   fprintf(fpsi2,"\n");
 279 }
 280
 281 static void set_ahcity(int nbb,t_bb bb[])
 282 {
 283   real pp2;
 284   int  n;
 285
 286   for(n=0; (n<nbb); n++) {
 287     pp2=sqr(bb[n].phi-PHI_AHX)+sqr(bb[n].psi-PSI_AHX);
 288
 289     bb[n].bHelix=FALSE;
 290     if (pp2 < 2500) {
 291       if ((bb[n].d4 < 0.36) || ((n > 0) && bb[n-1].bHelix))
 292         bb[n].bHelix=TRUE;
 293     }
 294   }
 295 }
 296
 297 t_bb *mkbbind(const char *fn,int *nres,int *nbb,int res0,
 298               int *nall,atom_id **index,
 299               char ***atomname,t_atom atom[],
 300               t_resinfo *resinfo)
 301 {
 302   static const char * bb_nm[] = { "N", "H", "CA", "C", "O", "HN" };
 303 #define NBB asize(bb_nm)
 304   t_bb    *bb;
 305   char    *grpname;
 306   int     ai,i,i0,i1,j,k,ri,rnr,gnx,r0,r1;
 307
 308   fprintf(stderr,"Please select a group containing the entire backbone\n");
 309   rd_index(fn,1,&gnx,index,&grpname);
 310   *nall=gnx;
 311   fprintf(stderr,"Checking group %s\n",grpname);
 312   r0=r1=atom[(*index)[0]].resind;
 313   for(i=1; (i<gnx); i++) {
 314     r0=min(r0,atom[(*index)[i]].resind);
 315     r1=max(r1,atom[(*index)[i]].resind);
 316   }
 317   rnr=r1-r0+1;
 318   fprintf(stderr,"There are %d residues\n",rnr);
 319   snew(bb,rnr);
 320   for(i=0; (i<rnr); i++)
 321     bb[i].N=bb[i].H=bb[i].CA=bb[i].C=bb[i].O=-1,bb[i].resno=res0+i;
 322
 323   for(i=j=0; (i<gnx); i++) {
 324     ai=(*index)[i];
 325     ri=atom[ai].resind-r0;
 326     if (strcmp(*(resinfo[ri].name),"PRO") == 0) {
 327       if (strcmp(*(atomname[ai]),"CD") == 0)
 328         bb[ri].H=ai;
 329     }
 330     for(k=0; (k<NBB); k++)
 331       if (strcmp(bb_nm[k],*(atomname[ai])) == 0) {
 332         j++;
 333         break;
 334       }
 335     switch (k) {
 336     case 0:
 337       bb[ri].N=ai;
 338       break;
 339     case 1:
 340     case 5:
 341       /* No attempt to address the case where some weird input has both H and HN atoms in the group */
 342       bb[ri].H=ai;
 343       break;
 344     case 2:
 345       bb[ri].CA=ai;
 346       break;
 347     case 3:
 348       bb[ri].C=ai;
 349       break;
 350     case 4:
 351       bb[ri].O=ai;
 352       break;
 353     default:
 354       break;
 355     }
 356   }
 357
 358   for(i0=0; (i0<rnr); i0++) {
 359     if ((bb[i0].N != -1) && (bb[i0].H != -1) &&
 360         (bb[i0].CA != -1) &&
 361         (bb[i0].C != -1) && (bb[i0].O != -1))
 362       break;
 363   }
 364   for(i1=rnr-1; (i1>=0); i1--) {
 365     if ((bb[i1].N != -1) && (bb[i1].H != -1) &&
 366         (bb[i1].CA != -1) &&
 367         (bb[i1].C != -1) && (bb[i1].O != -1))
 368       break;
 369   }
 370   if (i0 == 0)
 371     i0++;
 372   if (i1 == rnr-1)
 373     i1--;
 374
 375   for(i=i0; (i<i1); i++) {
 376     bb[i].Cprev=bb[i-1].C;
 377     bb[i].Nnext=bb[i+1].N;
 378   }
 379   rnr=max(0,i1-i0+1);
 380   fprintf(stderr,"There are %d complete backbone residues (from %d to %d)\n",
 381           rnr,bb[i0].resno,bb[i1].resno);
 382   if (rnr==0)
 383     gmx_fatal(FARGS,"Zero complete backbone residues were found, cannot proceed");
 384   for(i=0; (i<rnr); i++,i0++)
 385     bb[i]=bb[i0];
 386
 387   /* Set the labels */
 388   for(i=0; (i<rnr); i++) {
 389     ri=atom[bb[i].CA].resind;
 390     sprintf(bb[i].label,"%s%d",*(resinfo[ri].name),ri+res0);
 391   }
 392
 393   *nres=rnr;
 394   *nbb=rnr*asize(bb_nm);
 395
 396   return bb;
 397 }
 398
 399 real pprms(FILE *fp,int nbb,t_bb bb[])
 400 {
 401   int  i,n;
 402   real rms,rmst,rms2;
 403
 404   rmst=rms2=0;
 405   for(i=n=0; (i<nbb); i++) {
 406     if (bb[i].bHelix) {
 407       rms=sqrt(bb[i].pprms2);
 408       rmst+=rms;
 409       rms2+=bb[i].pprms2;
 410       fprintf(fp,"%10g  ",rms);
 411       n++;
 412     }
 413   }
 414   fprintf(fp,"\n");
 415   rms=sqrt(rms2/n-sqr(rmst/n));
 416
 417   return rms;
 418 }
 419
 420 void calc_hxprops(int nres,t_bb bb[],rvec x[],matrix box)
 421 {
 422   int  i,ao,an,t1,t2,t3;
 423   rvec dx,r_ij,r_kj,r_kl,m,n;
 424   real sign;
 425
 426   for(i=0; (i<nres); i++) {
 427     ao=bb[i].O;
 428     bb[i].d4=bb[i].d3=bb[i].d5=0;
 429     if (i < nres-3) {
 430       an=bb[i+3].N;
 431       rvec_sub(x[ao],x[an],dx);
 432       bb[i].d3=norm(dx);
 433     }
 434     if (i < nres-4) {
 435       an=bb[i+4].N;
 436       rvec_sub(x[ao],x[an],dx);
 437       bb[i].d4=norm(dx);
 438     }
 439     if (i < nres-5) {
 440       an=bb[i+5].N;
 441       rvec_sub(x[ao],x[an],dx);
 442       bb[i].d5=norm(dx);
 443     }
 444
 445     bb[i].phi=RAD2DEG*
 446       dih_angle(x[bb[i].Cprev],x[bb[i].N],x[bb[i].CA],x[bb[i].C],NULL,
 447                 r_ij,r_kj,r_kl,m,n,
 448                 &sign,&t1,&t2,&t3);
 449     bb[i].psi=RAD2DEG*
 450       dih_angle(x[bb[i].N],x[bb[i].CA],x[bb[i].C],x[bb[i].Nnext],NULL,
 451                 r_ij,r_kj,r_kl,m,n,
 452                 &sign,&t1,&t2,&t3);
 453     bb[i].pprms2=sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX);
 454
 455     bb[i].jcaha+=
 456       1.4*sin((bb[i].psi+138.0)*DEG2RAD) -
 457         4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
 458           2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0));
 459   }
 460 }
 461
 462 static void check_ahx(int nres,t_bb bb[],rvec x[],
 463                       int *hstart,int *hend)
 464 {
 465   int  h0,h1,h0sav,h1sav;
 466
 467   set_ahcity(nres,bb);
 468   h0=h0sav=h1sav=0;
 469   do {
 470     for(; (!bb[h0].bHelix) && (h0<nres-4); h0++)
 471       ;
 472     for(h1=h0; bb[h1+1].bHelix && (h1<nres-1); h1++)
 473       ;
 474     if (h1 > h0) {
 475       /*fprintf(stderr,"Helix from %d to %d\n",h0,h1);*/
 476       if (h1-h0 > h1sav-h0sav) {
 477         h0sav=h0;
 478         h1sav=h1;
 479       }
 480     }
 481     h0=h1+1;
 482   } while (h1 < nres-1);
 483   *hstart=h0sav;
 484   *hend=h1sav;
 485 }
 486
 487 void do_start_end(int nres,t_bb bb[],rvec x[],int *nbb,atom_id bbindex[],
 488                   int *nca,atom_id caindex[],
 489                   gmx_bool bRange,int rStart,int rEnd)
 490 {
 491   int    i,j,hstart=0,hend=0;
 492
 493   if (bRange) {
 494     for(i=0; (i<nres); i++) {
 495       if ((bb[i].resno >= rStart) && (bb[i].resno <= rEnd))
 496         bb[i].bHelix=TRUE;
 497       if (bb[i].resno == rStart)
 498         hstart=i;
 499       if (bb[i].resno == rEnd)
 500         hend=i;
 501     }
 502   }
 503   else {
 504     /* Find start and end of longest helix fragment */
 505     check_ahx(nres,bb,x,&hstart,&hend);
 506   }
 507   fprintf(stderr,"helix from: %d through %d\n",
 508           bb[hstart].resno,bb[hend].resno);
 509
 510   for(j=0,i=hstart; (i<=hend); i++) {
 511     bbindex[j++]=bb[i].N;
 512     bbindex[j++]=bb[i].H;
 513     bbindex[j++]=bb[i].CA;
 514     bbindex[j++]=bb[i].C;
 515     bbindex[j++]=bb[i].O;
 516     caindex[i-hstart]=bb[i].CA;
 517   }
 518   *nbb=j;
 519   *nca=(hend-hstart+1);
 520 }
 521
 522 void pr_bb(FILE *fp,int nres,t_bb bb[])
 523 {
 524   int i;
 525
 526   fprintf(fp,"\n");
 527   fprintf(fp,"%3s %3s %3s %3s %3s %7s %7s %7s %7s %7s %3s\n",
 528           "AA","N","Ca","C","O","Phi","Psi","D3","D4","D5","Hx?");
 529   for(i=0; (i<nres); i++) {
 530     fprintf(fp,"%3d %3d %3d %3d %3d %7.2f %7.2f %7.3f %7.3f %7.3f %3s\n",
 531             bb[i].resno,bb[i].N,bb[i].CA,bb[i].C,bb[i].O,
 532             bb[i].phi,bb[i].psi,bb[i].d3,bb[i].d4,bb[i].d5,
 533             bb[i].bHelix ? "Yes" : "No");
 534   }
 535   fprintf(fp,"\n");
 536 }
 537