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