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, 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, 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) || !std::isfinite(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_pbc             pbc, *pbc_null;
 715     int               my_clust;
 716     FILE             *fplog;
 717     gmx_output_env_t *oenv;
 718     gmx_rmpbc_t       gpbc = nullptr;
 719
 720     t_filenm          fnm[] = {
 721         { efTPR, nullptr, nullptr, ffREAD },
 722         { efTRX, "-f", nullptr, ffREAD },
 723         { efXVG, "-ds", "drsum",  ffWRITE },
 724         { efXVG, "-da", "draver", ffWRITE },
 725         { efXVG, "-dn", "drnum",  ffWRITE },
 726         { efXVG, "-dm", "drmax",  ffWRITE },
 727         { efXVG, "-dr", "restr",  ffWRITE },
 728         { efLOG, "-l",  "disres", ffWRITE },
 729         { efNDX, nullptr,  "viol",   ffOPTRD },
 730         { efPDB, "-q",  "viol",   ffOPTWR },
 731         { efNDX, "-c",  "clust",  ffOPTRD },
 732         { efXPM, "-x",  "matrix", ffOPTWR }
 733     };
 734 #define NFILE asize(fnm)
 735
 736     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 737                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 738     {
 739         return 0;
 740     }
 741
 742     fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
 743
 744     if (ntop)
 745     {
 746         init5(ntop);
 747     }
 748
 749     t_inputrec      irInstance;
 750     t_inputrec     *ir = &irInstance;
 751
 752     read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE);
 753     snew(xtop, header.natoms);
 754     read_tpx(ftp2fn(efTPR, NFILE, fnm), ir, box, &ntopatoms, xtop, nullptr, &mtop);
 755     bPDB = opt2bSet("-q", NFILE, fnm);
 756     if (bPDB)
 757     {
 758         snew(xav, ntopatoms);
 759         snew(ind_fit, ntopatoms);
 760         snew(w_rls, ntopatoms);
 761         for (kkk = 0; (kkk < ntopatoms); kkk++)
 762         {
 763             w_rls[kkk]   = 1;
 764             ind_fit[kkk] = kkk;
 765         }
 766
 767         snew(atoms, 1);
 768         *atoms = gmx_mtop_global_atoms(&mtop);
 769
 770         if (atoms->pdbinfo == nullptr)
 771         {
 772             snew(atoms->pdbinfo, atoms->nr);
 773         }
 774         atoms->havePdbInfo = TRUE;
 775     }
 776
 777     top = gmx_mtop_generate_local_top(&mtop, ir->efep != efepNO);
 778
 779     g        = nullptr;
 780     pbc_null = nullptr;
 781     if (ir->ePBC != epbcNONE)
 782     {
 783         if (ir->bPeriodicMols)
 784         {
 785             pbc_null = &pbc;
 786         }
 787         else
 788         {
 789             g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
 790         }
 791     }
 792
 793     if (ftp2bSet(efNDX, NFILE, fnm))
 794     {
 795         /* TODO: Nothing is written to this file if -c is provided, but it is
 796          * still opened... */
 797         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 798         xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
 799                        "nm", oenv);
 800         snew(vvindex, isize);
 801         snew(leg, isize);
 802         for (i = 0; (i < isize); i++)
 803         {
 804             index[i]++;
 805             snew(leg[i], 12);
 806             sprintf(leg[i], "index %d", index[i]);
 807         }
 808         xvgr_legend(xvg, isize, (const char**)leg, oenv);
 809     }
 810     else
 811     {
 812         isize = 0;
 813     }
 814
 815     ir->dr_tau = 0.0;
 816     init_disres(fplog, &mtop, ir, nullptr, nullptr, &fcd, nullptr, FALSE);
 817
 818     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 819     snew(f, 5*natoms);
 820
 821     init_dr_res(&dr, fcd.disres.nres);
 822     if (opt2bSet("-c", NFILE, fnm))
 823     {
 824         clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
 825         snew(dr_clust, clust->clust->nr+1);
 826         for (i = 0; (i <= clust->clust->nr); i++)
 827         {
 828             init_dr_res(&dr_clust[i], fcd.disres.nres);
 829         }
 830     }
 831     else
 832     {
 833         out = xvgropen(opt2fn("-ds", NFILE, fnm),
 834                        "Sum of Violations", "Time (ps)", "nm", oenv);
 835         aver = xvgropen(opt2fn("-da", NFILE, fnm),
 836                         "Average Violation", "Time (ps)", "nm", oenv);
 837         numv = xvgropen(opt2fn("-dn", NFILE, fnm),
 838                         "# Violations", "Time (ps)", "#", oenv);
 839         maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
 840                          "Largest Violation", "Time (ps)", "nm", oenv);
 841     }
 842
 843     auto mdAtoms = gmx::makeMDAtoms(fplog, mtop, *ir, false);
 844     atoms2md(&mtop, ir, -1, nullptr, mtop.natoms, mdAtoms.get());
 845     update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
 846     init_nrnb(&nrnb);
 847     if (ir->ePBC != epbcNONE)
 848     {
 849         gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
 850     }
 851
 852     j = 0;
 853     do
 854     {
 855         if (ir->ePBC != epbcNONE)
 856         {
 857             if (ir->bPeriodicMols)
 858             {
 859                 set_pbc(&pbc, ir->ePBC, box);
 860             }
 861             else
 862             {
 863                 gmx_rmpbc(gpbc, natoms, box, x);
 864             }
 865         }
 866
 867         if (clust)
 868         {
 869             if (j > clust->maxframe)
 870             {
 871                 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
 872             }
 873             my_clust = clust->inv_clust[j];
 874             range_check(my_clust, 0, clust->clust->nr);
 875             check_viol(fplog, &(top->idef.il[F_DISRES]),
 876                        top->idef.iparams,
 877                        x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
 878         }
 879         else
 880         {
 881             check_viol(fplog, &(top->idef.il[F_DISRES]),
 882                        top->idef.iparams,
 883                        x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
 884         }
 885         if (bPDB)
 886         {
 887             reset_x(atoms->nr, ind_fit, atoms->nr, nullptr, x, w_rls);
 888             do_fit(atoms->nr, w_rls, x, x);
 889             if (j == 0)
 890             {
 891                 /* Store the first frame of the trajectory as 'characteristic'
 892                  * for colouring with violations.
 893                  */
 894                 for (kkk = 0; (kkk < atoms->nr); kkk++)
 895                 {
 896                     copy_rvec(x[kkk], xav[kkk]);
 897                 }
 898             }
 899         }
 900         if (!clust)
 901         {
 902             if (isize > 0)
 903             {
 904                 fprintf(xvg, "%10g", t);
 905                 for (i = 0; (i < isize); i++)
 906                 {
 907                     fprintf(xvg, "  %10g", vvindex[i]);
 908                 }
 909                 fprintf(xvg, "\n");
 910             }
 911             fprintf(out,  "%10g  %10g\n", t, dr.sumv);
 912             fprintf(aver, "%10g  %10g\n", t, dr.averv);
 913             fprintf(maxxv, "%10g  %10g\n", t, dr.maxv);
 914             fprintf(numv, "%10g  %10d\n", t, dr.nv);
 915         }
 916         j++;
 917     }
 918     while (read_next_x(oenv, status, &t, x, box));
 919     close_trx(status);
 920     if (ir->ePBC != epbcNONE)
 921     {
 922         gmx_rmpbc_done(gpbc);
 923     }
 924
 925     if (clust)
 926     {
 927         dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 928                          top->idef.iparams, clust->clust, dr_clust,
 929                          clust->grpname, isize, index);
 930     }
 931     else
 932     {
 933         dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 934                    top->idef.iparams, &dr, isize, index,
 935                    bPDB ? atoms : nullptr);
 936         if (bPDB)
 937         {
 938             write_sto_conf(opt2fn("-q", NFILE, fnm),
 939                            "Coloured by average violation in Angstrom",
 940                            atoms, xav, nullptr, ir->ePBC, box);
 941         }
 942         dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
 943                           j, &top->idef, &mtop, max_dr, nlevels, bThird);
 944         xvgrclose(out);
 945         xvgrclose(aver);
 946         xvgrclose(numv);
 947         xvgrclose(maxxv);
 948         do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
 949         do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
 950         do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
 951         do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
 952     }
 953     if (isize > 0)
 954     {
 955         xvgrclose(xvg);
 956         if (!clust)
 957         {
 958             do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
 959         }
 960     }
 961
 962     gmx_ffclose(fplog);
 963
 964     return 0;
 965 }