src/gromacs/gmxana/hxprops.cpp

   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  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source 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 #include "gmxpre.h"
  39
  40 #include "hxprops.h"
  41
  42 #include <cmath>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/listed_forces/bonded.h"
  48 #include "gromacs/math/functions.h"
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/topology/index.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/utility/arraysize.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 real ellipticity(int nres, t_bb bb[])
  59 {
  60     typedef struct
  61     {
  62         real phi, psi, w;
  63     } t_ppwstr;
  64     // Avoid warnings about narrowing conversions from double to real
  65 #ifdef _MSC_VER
  66 #    pragma warning(disable : 4838)
  67 #endif
  68     static const t_ppwstr ppw[] = { { -67, -44, 0.31 },     { -66, -41, 0.31 }, { -59, -44, 0.44 },
  69                                     { -57, -47, 0.56 },     { -53, -52, 0.78 }, { -48, -57, 1.00 },
  70                                     { -70.5, -35.8, 0.15 }, { -57, -79, 0.23 }, { -38, -78, 1.20 },
  71                                     { -60, -30, 0.24 },     { -54, -28, 0.46 }, { -44, -33, 0.68 } };
  72 #ifdef _MSC_VER
  73 #    pragma warning(default : 4838)
  74 #endif
  75 #define NPPW asize(ppw)
  76
  77     int  i, j;
  78     real ell, pp2, phi, psi;
  79
  80     ell = 0;
  81     for (i = 0; (i < nres); i++)
  82     {
  83         phi = bb[i].phi;
  84         psi = bb[i].psi;
  85         for (j = 0; (j < NPPW); j++)
  86         {
  87             pp2 = gmx::square(phi - ppw[j].phi) + gmx::square(psi - ppw[j].psi);
  88             if (pp2 < 64)
  89             {
  90                 bb[i].nhx++;
  91                 ell += ppw[j].w;
  92                 break;
  93             }
  94         }
  95     }
  96     return ell;
  97 }
  98
  99 real ahx_len(int gnx, const int index[], rvec x[])
 100 /* Assume we have a list of Calpha atoms only! */
 101 {
 102     rvec dx;
 103
 104     rvec_sub(x[index[0]], x[index[gnx - 1]], dx);
 105
 106     return norm(dx);
 107 }
 108
 109 real radius(FILE* fp, int nca, const int ca_index[], rvec x[])
 110 /* Assume we have all the backbone */
 111 {
 112     real dl2, dlt;
 113     int  i, ai;
 114
 115     dlt = 0;
 116     for (i = 0; (i < nca); i++)
 117     {
 118         ai  = ca_index[i];
 119         dl2 = gmx::square(x[ai][XX]) + gmx::square(x[ai][YY]);
 120
 121         if (fp)
 122         {
 123             fprintf(fp, "  %10g", dl2);
 124         }
 125
 126         dlt += dl2;
 127     }
 128     if (fp)
 129     {
 130         fprintf(fp, "\n");
 131     }
 132
 133     return std::sqrt(dlt / nca);
 134 }
 135
 136 static real rot(rvec x1, const rvec x2)
 137 {
 138     real phi1, dphi, cp, sp;
 139     real xx, yy;
 140
 141     phi1 = std::atan2(x1[YY], x1[XX]);
 142     cp   = std::cos(phi1);
 143     sp   = std::sin(phi1);
 144     xx   = cp * x2[XX] + sp * x2[YY];
 145     yy   = -sp * x2[XX] + cp * x2[YY];
 146
 147     dphi = gmx::c_rad2Deg * std::atan2(yy, xx);
 148
 149     return dphi;
 150 }
 151
 152 real twist(int nca, const int caindex[], rvec x[])
 153 {
 154     real pt, dphi;
 155     int  i, a0, a1;
 156
 157     pt = 0;
 158     a0 = caindex[0];
 159     for (i = 1; (i < nca); i++)
 160     {
 161         a1 = caindex[i];
 162
 163         dphi = rot(x[a0], x[a1]);
 164         if (dphi < -90)
 165         {
 166             dphi += 360;
 167         }
 168         pt += dphi;
 169         a0 = a1;
 170     }
 171
 172     return (pt / (nca - 1));
 173 }
 174
 175 real ca_phi(int gnx, const int index[], rvec x[])
 176 /* Assume we have a list of Calpha atoms only! */
 177 {
 178     real phi, phitot;
 179     int  i, ai, aj, ak, al, t1, t2, t3;
 180     rvec r_ij, r_kj, r_kl, m, n;
 181
 182     if (gnx <= 4)
 183     {
 184         return 0;
 185     }
 186
 187     phitot = 0;
 188     for (i = 0; (i < gnx - 4); i++)
 189     {
 190         ai  = index[i + 0];
 191         aj  = index[i + 1];
 192         ak  = index[i + 2];
 193         al  = index[i + 3];
 194         phi = gmx::c_rad2Deg
 195               * dih_angle(x[ai], x[aj], x[ak], x[al], nullptr, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3);
 196         phitot += phi;
 197     }
 198
 199     return (phitot / (gnx - 4.0));
 200 }
 201
 202 real dip(int nbb, int const bbind[], const rvec x[], const 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, const int 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, FILE* fp4, FILE* fp4a, FILE* fp5, FILE* fp5a, real t, int nres, t_bb bb[])
 242 {
 243     int  i, n3 = 0, n4 = 0, n5 = 0;
 244     real d3 = 0, d4 = 0, d5 = 0;
 245
 246     for (i = 0; (i < nres - 3); i++)
 247     {
 248         if (bb[i].bHelix)
 249         {
 250             fprintf(fp3a, "%10g", bb[i].d3);
 251             n3++;
 252             d3 += bb[i].d3;
 253             if (i < nres - 4)
 254             {
 255                 fprintf(fp4a, "%10g", bb[i].d4);
 256                 n4++;
 257                 d4 += bb[i].d4;
 258             }
 259             if (i < nres - 5)
 260             {
 261                 fprintf(fp5a, "%10g", bb[i].d5);
 262                 n5++;
 263                 d5 += bb[i].d5;
 264             }
 265         }
 266     }
 267     fprintf(fp3, "%10g  %10g\n", t, d3 / n3);
 268     fprintf(fp4, "%10g  %10g\n", t, d4 / n4);
 269     fprintf(fp5, "%10g  %10g\n", t, d5 / n5);
 270     fprintf(fp3a, "\n");
 271     fprintf(fp4a, "\n");
 272     fprintf(fp5a, "\n");
 273 }
 274
 275 void av_phipsi(FILE* fphi, FILE* fpsi, FILE* fphi2, FILE* fpsi2, real t, int nres, t_bb bb[])
 276 {
 277     int  i, n = 0;
 278     real phi = 0, psi = 0;
 279
 280     fprintf(fphi2, "%10g", t);
 281     fprintf(fpsi2, "%10g", t);
 282     for (i = 0; (i < nres); i++)
 283     {
 284         if (bb[i].bHelix)
 285         {
 286             phi += bb[i].phi;
 287             psi += bb[i].psi;
 288             fprintf(fphi2, "  %10g", bb[i].phi);
 289             fprintf(fpsi2, "  %10g", bb[i].psi);
 290             n++;
 291         }
 292     }
 293     fprintf(fphi, "%10g  %10g\n", t, (phi / n));
 294     fprintf(fpsi, "%10g  %10g\n", t, (psi / n));
 295     fprintf(fphi2, "\n");
 296     fprintf(fpsi2, "\n");
 297 }
 298
 299 static void set_ahcity(int nbb, t_bb bb[])
 300 {
 301     real pp2;
 302     int  n;
 303
 304     for (n = 0; (n < nbb); n++)
 305     {
 306         pp2 = gmx::square(bb[n].phi - PHI_AHX) + gmx::square(bb[n].psi - PSI_AHX);
 307
 308         bb[n].bHelix = FALSE;
 309         if (pp2 < 2500)
 310         {
 311             if ((bb[n].d4 < 0.36) || ((n > 0) && bb[n - 1].bHelix))
 312             {
 313                 bb[n].bHelix = TRUE;
 314             }
 315         }
 316     }
 317 }
 318
 319 t_bb* mkbbind(const char* fn,
 320               int*        nres,
 321               int*        nbb,
 322               int         res0,
 323               int*        nall,
 324               int**       index,
 325               char***     atomname,
 326               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   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 (int i = 1; (i < gnx); i++)
 341     {
 342         r0 = std::min(r0, atom[(*index)[i]].resind);
 343         r1 = std::max(r1, atom[(*index)[i]].resind);
 344     }
 345     int rnr = r1 - r0 + 1;
 346     fprintf(stderr, "There are %d residues\n", rnr);
 347     snew(bb, rnr);
 348     for (int i = 0; (i < rnr); i++)
 349     {
 350         bb[i].N = bb[i].H = bb[i].CA = bb[i].C = bb[i].O = -1;
 351         bb[i].resno                                      = res0 + i;
 352     }
 353
 354     for (int i = 0; (i < gnx); i++)
 355     {
 356         int ai = (*index)[i];
 357         // Create an index into the residue index for the topology.
 358         int resindex = atom[ai].resind;
 359         // Create an index into the residues present in the selected
 360         // index group.
 361         int bbindex = resindex - r0;
 362         if (std::strcmp(*(resinfo[resindex].name), "PRO") == 0)
 363         {
 364             // For PRO in a peptide, there is no H bound to backbone
 365             // N, so use CD instead.
 366             if (std::strcmp(*(atomname[ai]), "CD") == 0)
 367             {
 368                 bb[bbindex].H = ai;
 369             }
 370         }
 371         int k = 0;
 372         for (; (k < NBB); k++)
 373         {
 374             if (std::strcmp(bb_nm[k], *(atomname[ai])) == 0)
 375             {
 376                 break;
 377             }
 378         }
 379         switch (k)
 380         {
 381             case 0: bb[bbindex].N = ai; break;
 382             case 1:
 383             case 5:
 384                 /* No attempt to address the case where some weird input has both H and HN atoms in the group */
 385                 bb[bbindex].H = ai;
 386                 break;
 387             case 2: bb[bbindex].CA = ai; break;
 388             case 3: bb[bbindex].C = ai; break;
 389             case 4: bb[bbindex].O = ai; break;
 390             default: break;
 391         }
 392     }
 393
 394     int i0 = 0;
 395     for (; (i0 < rnr); i0++)
 396     {
 397         if ((bb[i0].N != -1) && (bb[i0].H != -1) && (bb[i0].CA != -1) && (bb[i0].C != -1)
 398             && (bb[i0].O != -1))
 399         {
 400             break;
 401         }
 402     }
 403     int i1 = rnr - 1;
 404     for (; (i1 >= 0); i1--)
 405     {
 406         if ((bb[i1].N != -1) && (bb[i1].H != -1) && (bb[i1].CA != -1) && (bb[i1].C != -1)
 407             && (bb[i1].O != -1))
 408         {
 409             break;
 410         }
 411     }
 412     if (i0 == 0)
 413     {
 414         i0++;
 415     }
 416     if (i1 == rnr - 1)
 417     {
 418         i1--;
 419     }
 420
 421     for (int i = i0; (i < i1); i++)
 422     {
 423         bb[i].Cprev = bb[i - 1].C;
 424         bb[i].Nnext = bb[i + 1].N;
 425     }
 426     rnr = std::max(0, i1 - i0 + 1);
 427     fprintf(stderr, "There are %d complete backbone residues (from %d to %d)\n", rnr, bb[i0].resno, bb[i1].resno);
 428     if (rnr == 0)
 429     {
 430         gmx_fatal(FARGS, "Zero complete backbone residues were found, cannot proceed");
 431     }
 432     for (int i = 0; (i < rnr); i++, i0++)
 433     {
 434         bb[i] = bb[i0];
 435     }
 436
 437     /* Set the labels */
 438     for (int i = 0; (i < rnr); i++)
 439     {
 440         int resindex = atom[bb[i].CA].resind;
 441         sprintf(bb[i].label, "%s%d", *(resinfo[resindex].name), resinfo[resindex].nr);
 442     }
 443
 444     *nres = rnr;
 445     *nbb  = rnr * asize(bb_nm);
 446
 447     return bb;
 448 }
 449
 450 real pprms(FILE* fp, int nbb, t_bb bb[])
 451 {
 452     int  i, n;
 453     real rms, rmst, rms2;
 454
 455     rmst = rms2 = 0;
 456     for (i = n = 0; (i < nbb); i++)
 457     {
 458         if (bb[i].bHelix)
 459         {
 460             rms = std::sqrt(bb[i].pprms2);
 461             rmst += rms;
 462             rms2 += bb[i].pprms2;
 463             fprintf(fp, "%10g  ", rms);
 464             n++;
 465         }
 466     }
 467     fprintf(fp, "\n");
 468     rms = std::sqrt(rms2 / n - gmx::square(rmst / n));
 469
 470     return rms;
 471 }
 472
 473 void calc_hxprops(int nres, t_bb bb[], const rvec x[])
 474 {
 475     int  i, ao, an, t1, t2, t3;
 476     rvec dx, r_ij, r_kj, r_kl, m, n;
 477
 478     for (i = 0; (i < nres); i++)
 479     {
 480         ao       = bb[i].O;
 481         bb[i].d4 = bb[i].d3 = bb[i].d5 = 0;
 482         if (i < nres - 3)
 483         {
 484             an = bb[i + 3].N;
 485             rvec_sub(x[ao], x[an], dx);
 486             bb[i].d3 = norm(dx);
 487         }
 488         if (i < nres - 4)
 489         {
 490             an = bb[i + 4].N;
 491             rvec_sub(x[ao], x[an], dx);
 492             bb[i].d4 = norm(dx);
 493         }
 494         if (i < nres - 5)
 495         {
 496             an = bb[i + 5].N;
 497             rvec_sub(x[ao], x[an], dx);
 498             bb[i].d5 = norm(dx);
 499         }
 500
 501         bb[i].phi = gmx::c_rad2Deg
 502                     * dih_angle(x[bb[i].Cprev],
 503                                 x[bb[i].N],
 504                                 x[bb[i].CA],
 505                                 x[bb[i].C],
 506                                 nullptr,
 507                                 r_ij,
 508                                 r_kj,
 509                                 r_kl,
 510                                 m,
 511                                 n,
 512                                 &t1,
 513                                 &t2,
 514                                 &t3);
 515         bb[i].psi = gmx::c_rad2Deg
 516                     * dih_angle(x[bb[i].N],
 517                                 x[bb[i].CA],
 518                                 x[bb[i].C],
 519                                 x[bb[i].Nnext],
 520                                 nullptr,
 521                                 r_ij,
 522                                 r_kj,
 523                                 r_kl,
 524                                 m,
 525                                 n,
 526                                 &t1,
 527                                 &t2,
 528                                 &t3);
 529         bb[i].pprms2 = gmx::square(bb[i].phi - PHI_AHX) + gmx::square(bb[i].psi - PSI_AHX);
 530
 531         bb[i].jcaha += 1.4 * std::sin((bb[i].psi + 138.0) * gmx::c_deg2Rad)
 532                        - 4.1 * std::cos(2.0 * gmx::c_deg2Rad * (bb[i].psi + 138.0))
 533                        + 2.0 * std::cos(2.0 * gmx::c_deg2Rad * (bb[i].phi + 30.0));
 534     }
 535 }
 536
 537 static void check_ahx(int nres, t_bb bb[], int* hstart, int* hend)
 538 {
 539     int h0, h1, h0sav, h1sav;
 540
 541     set_ahcity(nres, bb);
 542     h0 = h0sav = h1sav = 0;
 543     do
 544     {
 545         for (; (!bb[h0].bHelix) && (h0 < nres - 4); h0++) {}
 546         for (h1 = h0; bb[h1 + 1].bHelix && (h1 < nres - 1); h1++) {}
 547         if (h1 > h0)
 548         {
 549             /*fprintf(stderr,"Helix from %d to %d\n",h0,h1);*/
 550             if (h1 - h0 > h1sav - h0sav)
 551             {
 552                 h0sav = h0;
 553                 h1sav = h1;
 554             }
 555         }
 556         h0 = h1 + 1;
 557     } while (h1 < nres - 1);
 558     *hstart = h0sav;
 559     *hend   = h1sav;
 560 }
 561
 562 void do_start_end(int nres, t_bb bb[], int* nbb, int bbindex[], int* nca, int caindex[], gmx_bool bRange, int rStart, int rEnd)
 563 {
 564     int i, j, hstart = 0, hend = 0;
 565
 566     if (bRange)
 567     {
 568         for (i = 0; (i < nres); i++)
 569         {
 570             if ((bb[i].resno >= rStart) && (bb[i].resno <= rEnd))
 571             {
 572                 bb[i].bHelix = TRUE;
 573             }
 574             if (bb[i].resno == rStart)
 575             {
 576                 hstart = i;
 577             }
 578             if (bb[i].resno == rEnd)
 579             {
 580                 hend = i;
 581             }
 582         }
 583     }
 584     else
 585     {
 586         /* Find start and end of longest helix fragment */
 587         check_ahx(nres, bb, &hstart, &hend);
 588     }
 589     fprintf(stderr, "helix from: %d through %d\n", bb[hstart].resno, bb[hend].resno);
 590
 591     for (j = 0, i = hstart; (i <= hend); i++)
 592     {
 593         bbindex[j++]        = bb[i].N;
 594         bbindex[j++]        = bb[i].H;
 595         bbindex[j++]        = bb[i].CA;
 596         bbindex[j++]        = bb[i].C;
 597         bbindex[j++]        = bb[i].O;
 598         caindex[i - hstart] = bb[i].CA;
 599     }
 600     *nbb = j;
 601     *nca = (hend - hstart + 1);
 602 }
 603
 604 void pr_bb(FILE* fp, int nres, t_bb bb[])
 605 {
 606     int i;
 607
 608     fprintf(fp, "\n");
 609     fprintf(fp,
 610             "%3s %3s %3s %3s %3s %7s %7s %7s %7s %7s %3s\n",
 611             "AA",
 612             "N",
 613             "Ca",
 614             "C",
 615             "O",
 616             "Phi",
 617             "Psi",
 618             "D3",
 619             "D4",
 620             "D5",
 621             "Hx?");
 622     for (i = 0; (i < nres); i++)
 623     {
 624         fprintf(fp,
 625                 "%3d %3d %3d %3d %3d %7.2f %7.2f %7.3f %7.3f %7.3f %3s\n",
 626                 bb[i].resno,
 627                 bb[i].N,
 628                 bb[i].CA,
 629                 bb[i].C,
 630                 bb[i].O,
 631                 bb[i].phi,
 632                 bb[i].psi,
 633                 bb[i].d3,
 634                 bb[i].d4,
 635                 bb[i].d5,
 636                 bb[i].bHelix ? "Yes" : "No");
 637     }
 638     fprintf(fp, "\n");
 639 }