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