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