src/gromacs/gmxana/gmx_disre.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,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstdlib>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/commandline/viewit.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/matio.h"
  49 #include "gromacs/fileio/pdbio.h"
  50 #include "gromacs/fileio/tpxio.h"
  51 #include "gromacs/fileio/trxio.h"
  52 #include "gromacs/fileio/xvgr.h"
  53 #include "gromacs/gmxana/gmx_ana.h"
  54 #include "gromacs/gmxana/gstat.h"
  55 #include "gromacs/gmxlib/nrnb.h"
  56 #include "gromacs/listed_forces/disre.h"
  57 #include "gromacs/math/do_fit.h"
  58 #include "gromacs/math/functions.h"
  59 #include "gromacs/math/vec.h"
  60 #include "gromacs/mdlib/force.h"
  61 #include "gromacs/mdlib/mdatoms.h"
  62 #include "gromacs/mdlib/mdrun.h"
  63 #include "gromacs/mdtypes/fcdata.h"
  64 #include "gromacs/mdtypes/inputrec.h"
  65 #include "gromacs/mdtypes/md_enums.h"
  66 #include "gromacs/pbcutil/ishift.h"
  67 #include "gromacs/pbcutil/mshift.h"
  68 #include "gromacs/pbcutil/pbc.h"
  69 #include "gromacs/pbcutil/rmpbc.h"
  70 #include "gromacs/topology/index.h"
  71 #include "gromacs/topology/mtop_util.h"
  72 #include "gromacs/topology/topology.h"
  73 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  74 #include "gromacs/utility/arraysize.h"
  75 #include "gromacs/utility/fatalerror.h"
  76 #include "gromacs/utility/futil.h"
  77 #include "gromacs/utility/smalloc.h"
  78
  79 typedef struct {
  80     int  n;
  81     real v;
  82 } t_toppop;
  83
  84 static t_toppop *top  = nullptr;
  85 static int       ntop = 0;
  86
  87 typedef struct {
  88     int   nv, nframes;
  89     real  sumv, averv, maxv;
  90     real *aver1, *aver2, *aver_3, *aver_6;
  91 } t_dr_result;
  92
  93 static void init5(int n)
  94 {
  95     ntop = n;
  96     snew(top, ntop);
  97 }
  98
  99 static void reset5()
 100 {
 101     int i;
 102
 103     for (i = 0; (i < ntop); i++)
 104     {
 105         top[i].n = -1;
 106         top[i].v = 0;
 107     }
 108 }
 109
 110 static void add5(int ndr, real viol)
 111 {
 112     int i, mini;
 113
 114     mini = 0;
 115     for (i = 1; (i < ntop); i++)
 116     {
 117         if (top[i].v < top[mini].v)
 118         {
 119             mini = i;
 120         }
 121     }
 122     if (viol > top[mini].v)
 123     {
 124         top[mini].v = viol;
 125         top[mini].n = ndr;
 126     }
 127 }
 128
 129 static void print5(FILE *fp)
 130 {
 131     int i;
 132
 133     std::sort(top, top+ntop, [](const t_toppop &a, const t_toppop &b) {return a.v > b.v; }); //reverse sort
 134     fprintf(fp, "Index:");
 135     for (i = 0; (i < ntop); i++)
 136     {
 137         fprintf(fp, " %6d", top[i].n);
 138     }
 139     fprintf(fp, "\nViol: ");
 140     for (i = 0; (i < ntop); i++)
 141     {
 142         fprintf(fp, " %6.3f", top[i].v);
 143     }
 144     fprintf(fp, "\n");
 145 }
 146
 147 static void check_viol(FILE *log,
 148                        t_ilist *disres, t_iparams forceparams[],
 149                        rvec x[], rvec4 f[],
 150                        t_pbc *pbc, t_graph *g, t_dr_result dr[],
 151                        int clust_id, int isize, const int index[], real vvindex[],
 152                        t_fcdata *fcd)
 153 {
 154     t_iatom         *forceatoms;
 155     int              i, j, nat, n, type, nviol, ndr, label;
 156     real             rt, mviol, tviol, viol, lam, dvdl, drt;
 157     rvec            *fshift;
 158     static  gmx_bool bFirst = TRUE;
 159
 160     lam   = 0;
 161     dvdl  = 0;
 162     tviol = 0;
 163     nviol = 0;
 164     mviol = 0;
 165     ndr   = 0;
 166     if (ntop)
 167     {
 168         reset5();
 169     }
 170     forceatoms = disres->iatoms;
 171     for (j = 0; (j < isize); j++)
 172     {
 173         vvindex[j] = 0;
 174     }
 175     nat = interaction_function[F_DISRES].nratoms+1;
 176     for (i = 0; (i < disres->nr); )
 177     {
 178         type  = forceatoms[i];
 179         n     = 0;
 180         label = forceparams[type].disres.label;
 181         if (debug)
 182         {
 183             fprintf(debug, "DISRE: ndr = %d, label = %d  i=%d, n =%d\n",
 184                     ndr, label, i, n);
 185         }
 186         if (ndr != label)
 187         {
 188             gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
 189         }
 190         do
 191         {
 192             n += nat;
 193         }
 194         while (((i+n) < disres->nr) &&
 195                (forceparams[forceatoms[i+n]].disres.label == label));
 196
 197         calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i],
 198                         x, pbc, fcd, nullptr);
 199
 200         if (fcd->disres.Rt_6[label] <= 0)
 201         {
 202             gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[label]);
 203         }
 204
 205         rt = gmx::invsixthroot(fcd->disres.Rt_6[label]);
 206         dr[clust_id].aver1[ndr]  += rt;
 207         dr[clust_id].aver2[ndr]  += gmx::square(rt);
 208         drt = 1.0/gmx::power3(rt);
 209         dr[clust_id].aver_3[ndr] += drt;
 210         dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label];
 211
 212         snew(fshift, SHIFTS);
 213         ta_disres(n, &forceatoms[i], forceparams,
 214                   x, f, fshift,
 215                   pbc, g, lam, &dvdl, nullptr, fcd, nullptr);
 216         sfree(fshift);
 217         viol = fcd->disres.sumviol;
 218
 219         if (viol > 0)
 220         {
 221             nviol++;
 222             if (ntop)
 223             {
 224                 add5(forceparams[type].disres.label, viol);
 225             }
 226             if (viol > mviol)
 227             {
 228                 mviol = viol;
 229             }
 230             tviol += viol;
 231             for (j = 0; (j < isize); j++)
 232             {
 233                 if (index[j] == forceparams[type].disres.label)
 234                 {
 235                     vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[label]);
 236                 }
 237             }
 238         }
 239         ndr++;
 240         i   += n;
 241     }
 242     dr[clust_id].nv    = nviol;
 243     dr[clust_id].maxv  = mviol;
 244     dr[clust_id].sumv  = tviol;
 245     dr[clust_id].averv = tviol/ndr;
 246     dr[clust_id].nframes++;
 247
 248     if (bFirst)
 249     {
 250         fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
 251                 ndr, disres->nr/nat);
 252         bFirst = FALSE;
 253     }
 254     if (ntop)
 255     {
 256         print5(log);
 257     }
 258 }
 259
 260 typedef struct {
 261     int      index;
 262     gmx_bool bCore;
 263     real     up1, r, rT3, rT6, viol, violT3, violT6;
 264 } t_dr_stats;
 265
 266 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
 267 {
 268     static const char *core[] = { "All restraints", "Core restraints" };
 269     static const char *tp[]   = { "linear", "third power", "sixth power" };
 270     real               viol_tot, viol_max, viol = 0;
 271     gmx_bool           bCore;
 272     int                nviol, nrestr;
 273     int                i, kkk;
 274
 275     for (int iCore = 0; iCore < 2; iCore++)
 276     {
 277         bCore = (iCore == 1);
 278         for (kkk = 0; (kkk < 3); kkk++)
 279         {
 280             viol_tot  = 0;
 281             viol_max  = 0;
 282             nviol     = 0;
 283             nrestr    = 0;
 284             for (i = 0; (i < ndr); i++)
 285             {
 286                 if (!bCore || drs[i].bCore)
 287                 {
 288                     switch (kkk)
 289                     {
 290                         case 0:
 291                             viol = drs[i].viol;
 292                             break;
 293                         case 1:
 294                             viol = drs[i].violT3;
 295                             break;
 296                         case 2:
 297                             viol = drs[i].violT6;
 298                             break;
 299                         default:
 300                             gmx_incons("Dumping violations");
 301                     }
 302                     viol_max     = std::max(viol_max, viol);
 303                     if (viol > 0)
 304                     {
 305                         nviol++;
 306                     }
 307                     viol_tot  += viol;
 308                     nrestr++;
 309                 }
 310             }
 311             if ((nrestr > 0) || (bCore && (nrestr < ndr)))
 312             {
 313                 fprintf(log, "\n");
 314                 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
 315                 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
 316                 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
 317                 if (nrestr > 0)
 318                 {
 319                     fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
 320                 }
 321                 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
 322                 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
 323             }
 324         }
 325     }
 326 }
 327
 328 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
 329 {
 330     int i;
 331
 332     fprintf(log, "Restr. Core     Up1     <r>   <rT3>   <rT6>  <viol><violT3><violT6>\n");
 333     for (i = 0; (i < ndr); i++)
 334     {
 335         if (bLinear  && (drs[i].viol == 0))
 336         {
 337             break;
 338         }
 339         fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
 340                 drs[i].index, yesno_names[drs[i].bCore],
 341                 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
 342                 drs[i].viol, drs[i].violT3, drs[i].violT6);
 343     }
 344 }
 345
 346 static gmx_bool is_core(int i, int isize, const int index[])
 347 {
 348     int      kk;
 349     gmx_bool bIC = FALSE;
 350
 351     for (kk = 0; !bIC && (kk < isize); kk++)
 352     {
 353         bIC = (index[kk] == i);
 354     }
 355
 356     return bIC;
 357 }
 358
 359 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
 360                        t_iparams ip[], t_dr_result *dr,
 361                        int isize, int index[], t_atoms *atoms)
 362 {
 363     int         i, j, nra;
 364     t_dr_stats *drs;
 365
 366     fprintf(log, "\n");
 367     fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
 368     snew(drs, ndr);
 369     j         = 0;
 370     nra       = interaction_function[F_DISRES].nratoms+1;
 371     for (i = 0; (i < ndr); i++)
 372     {
 373         /* Search for the appropriate restraint */
 374         for (; (j < disres->nr); j += nra)
 375         {
 376             if (ip[disres->iatoms[j]].disres.label == i)
 377             {
 378                 break;
 379             }
 380         }
 381         drs[i].index  = i;
 382         drs[i].bCore  = is_core(i, isize, index);
 383         drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
 384         drs[i].r      = dr->aver1[i]/nsteps;
 385         drs[i].rT3    = gmx::invcbrt(dr->aver_3[i]/nsteps);
 386         drs[i].rT6    = gmx::invsixthroot(dr->aver_6[i]/nsteps);
 387         drs[i].viol   = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
 388         drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
 389         drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
 390         if (atoms)
 391         {
 392             int j1 = disres->iatoms[j+1];
 393             int j2 = disres->iatoms[j+2];
 394             atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
 395             atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
 396         }
 397     }
 398     dump_viol(log, ndr, drs, FALSE);
 399
 400     fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
 401     std::sort(drs, drs+ndr, [](const t_dr_stats &a, const t_dr_stats &b)
 402               {return a.viol > b.viol; });            //Reverse sort
 403     dump_viol(log, ndr, drs, TRUE);
 404
 405     dump_dump(log, ndr, drs);
 406
 407     sfree(drs);
 408 }
 409
 410 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
 411                              t_iparams ip[], t_blocka *clust, t_dr_result dr[],
 412                              char *clust_name[], int isize, int index[])
 413 {
 414     int         i, j, k, nra, mmm = 0;
 415     double      sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
 416     t_dr_stats *drs;
 417
 418     fprintf(fp, "\n");
 419     fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
 420     fprintf(fp, "Cluster  NFrames    SumV      MaxV     SumVT     MaxVT     SumVS     MaxVS\n");
 421
 422     snew(drs, ndr);
 423
 424     for (k = 0; (k < clust->nr); k++)
 425     {
 426         if (dr[k].nframes == 0)
 427         {
 428             continue;
 429         }
 430         if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
 431         {
 432             gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
 433                       "Found %d frames in trajectory rather than the expected %d\n",
 434                       clust_name[k], dr[k].nframes,
 435                       clust->index[k+1]-clust->index[k]);
 436         }
 437         if (!clust_name[k])
 438         {
 439             gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
 440         }
 441         j         = 0;
 442         nra       = interaction_function[F_DISRES].nratoms+1;
 443         sumV      = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
 444         for (i = 0; (i < ndr); i++)
 445         {
 446             /* Search for the appropriate restraint */
 447             for (; (j < disres->nr); j += nra)
 448             {
 449                 if (ip[disres->iatoms[j]].disres.label == i)
 450                 {
 451                     break;
 452                 }
 453             }
 454             drs[i].index  = i;
 455             drs[i].bCore  = is_core(i, isize, index);
 456             drs[i].up1    = ip[disres->iatoms[j]].disres.up1;
 457             drs[i].r      = dr[k].aver1[i]/dr[k].nframes;
 458             if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
 459             {
 460                 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
 461             }
 462             drs[i].rT3    = gmx::invcbrt(dr[k].aver_3[i]/dr[k].nframes);
 463             drs[i].rT6    = gmx::invsixthroot(dr[k].aver_6[i]/dr[k].nframes);
 464             drs[i].viol   = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
 465             drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
 466             drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
 467             sumV         += drs[i].viol;
 468             sumVT3       += drs[i].violT3;
 469             sumVT6       += drs[i].violT6;
 470             maxV          = std::max(maxV, static_cast<double>(drs[i].viol));
 471             maxVT3        = std::max(maxVT3, static_cast<double>(drs[i].violT3));
 472             maxVT6        = std::max(maxVT6, static_cast<double>(drs[i].violT6));
 473         }
 474         if (std::strcmp(clust_name[k], "1000") == 0)
 475         {
 476             mmm++;
 477         }
 478         fprintf(fp, "%-10s%6d%8.3f  %8.3f  %8.3f  %8.3f  %8.3f  %8.3f\n",
 479                 clust_name[k],
 480                 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
 481
 482     }
 483     fflush(fp);
 484     sfree(drs);
 485 }
 486
 487 static void init_dr_res(t_dr_result *dr, int ndr)
 488 {
 489     snew(dr->aver1, ndr+1);
 490     snew(dr->aver2, ndr+1);
 491     snew(dr->aver_3, ndr+1);
 492     snew(dr->aver_6, ndr+1);
 493     dr->nv      = 0;
 494     dr->nframes = 0;
 495     dr->sumv    = 0;
 496     dr->maxv    = 0;
 497     dr->averv   = 0;
 498 }
 499
 500 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
 501                               int nsteps, t_idef *idef, const gmx_mtop_t *mtop,
 502                               real max_dr, int nlevels, gmx_bool bThird)
 503 {
 504     FILE      *fp;
 505     int       *resnr;
 506     int        n_res, a_offset, mol, a;
 507     int        i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
 508     int        ai, aj, *ptr;
 509     real     **matrix, *t_res, hi, *w_dr, rav, rviol;
 510     t_rgb      rlo = { 1, 1, 1 };
 511     t_rgb      rhi = { 0, 0, 0 };
 512     if (fn == nullptr)
 513     {
 514         return;
 515     }
 516
 517     snew(resnr, mtop->natoms);
 518     n_res    = 0;
 519     a_offset = 0;
 520     for (const gmx_molblock_t &molb : mtop->molblock)
 521     {
 522         const t_atoms &atoms = mtop->moltype[molb.type].atoms;
 523         for (mol = 0; mol < molb.nmol; mol++)
 524         {
 525             for (a = 0; a < atoms.nr; a++)
 526             {
 527                 resnr[a_offset + a] = n_res + atoms.atom[a].resind;
 528             }
 529             n_res    += atoms.nres;
 530             a_offset += atoms.nr;
 531         }
 532     }
 533
 534     snew(t_res, n_res);
 535     for (i = 0; (i < n_res); i++)
 536     {
 537         t_res[i] = i+1;
 538     }
 539     snew(matrix, n_res);
 540     for (i = 0; (i < n_res); i++)
 541     {
 542         snew(matrix[i], n_res);
 543     }
 544     nratoms = interaction_function[F_DISRES].nratoms;
 545     nra     = (idef->il[F_DISRES].nr/(nratoms+1));
 546     snew(ptr, nra+1);
 547     index   = 0;
 548     nlabel  = 0;
 549     ptr[0]  = 0;
 550     snew(w_dr, ndr);
 551     for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
 552     {
 553         tp       = idef->il[F_DISRES].iatoms[i];
 554         label    = idef->iparams[tp].disres.label;
 555
 556         if (label != index)
 557         {
 558             /* Set index pointer */
 559             ptr[index+1] = i;
 560             if (nlabel <= 0)
 561             {
 562                 gmx_fatal(FARGS, "nlabel is %d, label = %d", nlabel, label);
 563             }
 564             if (index >= ndr)
 565             {
 566                 gmx_fatal(FARGS, "ndr = %d, index = %d", ndr, index);
 567             }
 568             /* Update the weight */
 569             w_dr[index] = 1.0/nlabel;
 570             index       = label;
 571             nlabel      = 1;
 572         }
 573         else
 574         {
 575             nlabel++;
 576         }
 577     }
 578     printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
 579     hi = 0;
 580     for (i = 0; (i < ndr); i++)
 581     {
 582         for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
 583         {
 584             tp  = idef->il[F_DISRES].iatoms[j];
 585             ai  = idef->il[F_DISRES].iatoms[j+1];
 586             aj  = idef->il[F_DISRES].iatoms[j+2];
 587
 588             ri = resnr[ai];
 589             rj = resnr[aj];
 590             if (bThird)
 591             {
 592                 rav = gmx::invcbrt(dr->aver_3[i]/nsteps);
 593             }
 594             else
 595             {
 596                 rav = dr->aver1[i]/nsteps;
 597             }
 598             if (debug)
 599             {
 600                 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
 601             }
 602             rviol           = std::max(0.0_real, rav-idef->iparams[tp].disres.up1);
 603             matrix[ri][rj] += w_dr[i]*rviol;
 604             matrix[rj][ri] += w_dr[i]*rviol;
 605             hi              = std::max(hi, matrix[ri][rj]);
 606             hi              = std::max(hi, matrix[rj][ri]);
 607         }
 608     }
 609
 610     sfree(resnr);
 611
 612     if (max_dr > 0)
 613     {
 614         if (hi > max_dr)
 615         {
 616             printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
 617         }
 618         hi = max_dr;
 619     }
 620     printf("Highest level in the matrix will be %g\n", hi);
 621     fp = gmx_ffopen(fn, "w");
 622     write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
 623               n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
 624     gmx_ffclose(fp);
 625 }
 626
 627 int gmx_disre(int argc, char *argv[])
 628 {
 629     const char       *desc[] = {
 630         "[THISMODULE] computes violations of distance restraints.",
 631         "The program always",
 632         "computes the instantaneous violations rather than time-averaged,",
 633         "because this analysis is done from a trajectory file afterwards",
 634         "it does not make sense to use time averaging. However,",
 635         "the time averaged values per restraint are given in the log file.[PAR]",
 636         "An index file may be used to select specific restraints for",
 637         "printing.[PAR]",
 638         "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
 639         "amount of average violations.[PAR]",
 640         "When the [TT]-c[tt] option is given, an index file will be read",
 641         "containing the frames in your trajectory corresponding to the clusters",
 642         "(defined in another manner) that you want to analyze. For these clusters",
 643         "the program will compute average violations using the third power",
 644         "averaging algorithm and print them in the log file."
 645     };
 646     static int        ntop      = 0;
 647     static int        nlevels   = 20;
 648     static real       max_dr    = 0;
 649     static gmx_bool   bThird    = TRUE;
 650     t_pargs           pa[]      = {
 651         { "-ntop", FALSE, etINT,  {&ntop},
 652           "Number of large violations that are stored in the log file every step" },
 653         { "-maxdr", FALSE, etREAL, {&max_dr},
 654           "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
 655         { "-nlevels", FALSE, etINT, {&nlevels},
 656           "Number of levels in the matrix output" },
 657         { "-third", FALSE, etBOOL, {&bThird},
 658           "Use inverse third power averaging or linear for matrix output" }
 659     };
 660
 661     FILE             *out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
 662     gmx_localtop_t    top;
 663     t_fcdata          fcd;
 664     t_nrnb            nrnb;
 665     t_graph          *g;
 666     int               i, j, kkk;
 667     t_trxstatus      *status;
 668     real              t;
 669     rvec             *x, *xav = nullptr;
 670     rvec4            *f;
 671     matrix            box;
 672     gmx_bool          bPDB;
 673     int               isize;
 674     int              *index = nullptr, *ind_fit = nullptr;
 675     char             *grpname;
 676     t_cluster_ndx    *clust = nullptr;
 677     t_dr_result       dr, *dr_clust = nullptr;
 678     char            **leg;
 679     real             *vvindex = nullptr, *w_rls = nullptr;
 680     t_pbc             pbc, *pbc_null;
 681     int               my_clust;
 682     FILE             *fplog;
 683     gmx_output_env_t *oenv;
 684     gmx_rmpbc_t       gpbc = nullptr;
 685
 686     t_filenm          fnm[] = {
 687         { efTPR, nullptr, nullptr, ffREAD },
 688         { efTRX, "-f", nullptr, ffREAD },
 689         { efXVG, "-ds", "drsum",  ffWRITE },
 690         { efXVG, "-da", "draver", ffWRITE },
 691         { efXVG, "-dn", "drnum",  ffWRITE },
 692         { efXVG, "-dm", "drmax",  ffWRITE },
 693         { efXVG, "-dr", "restr",  ffWRITE },
 694         { efLOG, "-l",  "disres", ffWRITE },
 695         { efNDX, nullptr,  "viol",   ffOPTRD },
 696         { efPDB, "-q",  "viol",   ffOPTWR },
 697         { efNDX, "-c",  "clust",  ffOPTRD },
 698         { efXPM, "-x",  "matrix", ffOPTWR }
 699     };
 700 #define NFILE asize(fnm)
 701
 702     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 703                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 704     {
 705         return 0;
 706     }
 707
 708     fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
 709
 710     if (ntop)
 711     {
 712         init5(ntop);
 713     }
 714
 715     t_inputrec               irInstance;
 716     t_inputrec              *ir = &irInstance;
 717
 718     gmx::TopologyInformation topInfo;
 719     topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
 720     int                      ntopatoms = topInfo.mtop()->natoms;
 721     AtomsDataPtr             atoms;
 722     bPDB = opt2bSet("-q", NFILE, fnm);
 723     if (bPDB)
 724     {
 725         snew(xav, ntopatoms);
 726         snew(ind_fit, ntopatoms);
 727         snew(w_rls, ntopatoms);
 728         for (kkk = 0; (kkk < ntopatoms); kkk++)
 729         {
 730             w_rls[kkk]   = 1;
 731             ind_fit[kkk] = kkk;
 732         }
 733
 734         atoms = topInfo.copyAtoms();
 735
 736         if (atoms->pdbinfo == nullptr)
 737         {
 738             snew(atoms->pdbinfo, atoms->nr);
 739         }
 740         atoms->havePdbInfo = TRUE;
 741     }
 742
 743     gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
 744
 745     g        = nullptr;
 746     pbc_null = nullptr;
 747     if (ir->ePBC != epbcNONE)
 748     {
 749         if (ir->bPeriodicMols)
 750         {
 751             pbc_null = &pbc;
 752         }
 753         else
 754         {
 755             g = mk_graph(fplog, &top.idef, 0, ntopatoms, FALSE, FALSE);
 756         }
 757     }
 758
 759     if (ftp2bSet(efNDX, NFILE, fnm))
 760     {
 761         /* TODO: Nothing is written to this file if -c is provided, but it is
 762          * still opened... */
 763         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 764         xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
 765                        "nm", oenv);
 766         snew(vvindex, isize);
 767         snew(leg, isize);
 768         for (i = 0; (i < isize); i++)
 769         {
 770             index[i]++;
 771             snew(leg[i], 12);
 772             sprintf(leg[i], "index %d", index[i]);
 773         }
 774         xvgr_legend(xvg, isize, leg, oenv);
 775     }
 776     else
 777     {
 778         isize = 0;
 779     }
 780
 781     ir->dr_tau = 0.0;
 782     init_disres(fplog, topInfo.mtop(), ir, nullptr, nullptr, &fcd, nullptr, FALSE);
 783
 784     int natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 785     snew(f, 5*natoms);
 786
 787     init_dr_res(&dr, fcd.disres.nres);
 788     if (opt2bSet("-c", NFILE, fnm))
 789     {
 790         clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
 791         snew(dr_clust, clust->clust->nr+1);
 792         for (i = 0; (i <= clust->clust->nr); i++)
 793         {
 794             init_dr_res(&dr_clust[i], fcd.disres.nres);
 795         }
 796     }
 797     else
 798     {
 799         out = xvgropen(opt2fn("-ds", NFILE, fnm),
 800                        "Sum of Violations", "Time (ps)", "nm", oenv);
 801         aver = xvgropen(opt2fn("-da", NFILE, fnm),
 802                         "Average Violation", "Time (ps)", "nm", oenv);
 803         numv = xvgropen(opt2fn("-dn", NFILE, fnm),
 804                         "# Violations", "Time (ps)", "#", oenv);
 805         maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
 806                          "Largest Violation", "Time (ps)", "nm", oenv);
 807     }
 808
 809     auto mdAtoms = gmx::makeMDAtoms(fplog, *topInfo.mtop(), *ir, false);
 810     atoms2md(topInfo.mtop(), ir, -1, nullptr, ntopatoms, mdAtoms.get());
 811     update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
 812     init_nrnb(&nrnb);
 813     if (ir->ePBC != epbcNONE)
 814     {
 815         gpbc = gmx_rmpbc_init(&top.idef, ir->ePBC, natoms);
 816     }
 817
 818     j = 0;
 819     do
 820     {
 821         if (ir->ePBC != epbcNONE)
 822         {
 823             if (ir->bPeriodicMols)
 824             {
 825                 set_pbc(&pbc, ir->ePBC, box);
 826             }
 827             else
 828             {
 829                 gmx_rmpbc(gpbc, natoms, box, x);
 830             }
 831         }
 832
 833         if (clust)
 834         {
 835             if (j > clust->maxframe)
 836             {
 837                 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
 838             }
 839             my_clust = clust->inv_clust[j];
 840             range_check(my_clust, 0, clust->clust->nr);
 841             check_viol(fplog, &(top.idef.il[F_DISRES]),
 842                        top.idef.iparams,
 843                        x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
 844         }
 845         else
 846         {
 847             check_viol(fplog, &(top.idef.il[F_DISRES]),
 848                        top.idef.iparams,
 849                        x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
 850         }
 851         if (bPDB)
 852         {
 853             reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
 854             do_fit(atoms->nr, w_rls, x, x);
 855             if (j == 0)
 856             {
 857                 /* Store the first frame of the trajectory as 'characteristic'
 858                  * for colouring with violations.
 859                  */
 860                 for (kkk = 0; (kkk < atoms->nr); kkk++)
 861                 {
 862                     copy_rvec(x[kkk], xav[kkk]);
 863                 }
 864             }
 865         }
 866         if (!clust)
 867         {
 868             if (isize > 0)
 869             {
 870                 fprintf(xvg, "%10g", t);
 871                 for (i = 0; (i < isize); i++)
 872                 {
 873                     fprintf(xvg, "  %10g", vvindex[i]);
 874                 }
 875                 fprintf(xvg, "\n");
 876             }
 877             fprintf(out,  "%10g  %10g\n", t, dr.sumv);
 878             fprintf(aver, "%10g  %10g\n", t, dr.averv);
 879             fprintf(maxxv, "%10g  %10g\n", t, dr.maxv);
 880             fprintf(numv, "%10g  %10d\n", t, dr.nv);
 881         }
 882         j++;
 883     }
 884     while (read_next_x(oenv, status, &t, x, box));
 885     close_trx(status);
 886     if (ir->ePBC != epbcNONE)
 887     {
 888         gmx_rmpbc_done(gpbc);
 889     }
 890
 891     if (clust)
 892     {
 893         dump_clust_stats(fplog, fcd.disres.nres, &(top.idef.il[F_DISRES]),
 894                          top.idef.iparams, clust->clust, dr_clust,
 895                          clust->grpname, isize, index);
 896     }
 897     else
 898     {
 899         dump_stats(fplog, j, fcd.disres.nres, &(top.idef.il[F_DISRES]),
 900                    top.idef.iparams, &dr, isize, index,
 901                    bPDB ? atoms.get() : nullptr);
 902         if (bPDB)
 903         {
 904             write_sto_conf(opt2fn("-q", NFILE, fnm),
 905                            "Coloured by average violation in Angstrom",
 906                            atoms.get(), xav, nullptr, ir->ePBC, box);
 907         }
 908         dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
 909                           j, &top.idef, topInfo.mtop(), max_dr, nlevels, bThird);
 910         xvgrclose(out);
 911         xvgrclose(aver);
 912         xvgrclose(numv);
 913         xvgrclose(maxxv);
 914         do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
 915         do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
 916         do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
 917         do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
 918     }
 919     if (isize > 0)
 920     {
 921         xvgrclose(xvg);
 922         if (!clust)
 923         {
 924             do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
 925         }
 926     }
 927
 928     gmx_ffclose(fplog);
 929
 930     return 0;
 931 }