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, 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  = NULL;
  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 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(NULL, n, &forceatoms[i],
 208                         (const rvec*)x, pbc, fcd, NULL);
 209
 210         if (fcd->disres.Rt_6[0] <= 0)
 211         {
 212             gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
 213         }
 214
 215         rt = gmx::invsixthroot(fcd->disres.Rt_6[0]);
 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[0];
 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, NULL, fcd, NULL);
 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[0]);
 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 == NULL)
 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 = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
 692     t_tpxheader       header;
 693     t_inputrec        ir;
 694     gmx_mtop_t        mtop;
 695     rvec             *xtop;
 696     gmx_localtop_t   *top;
 697     t_atoms          *atoms = NULL;
 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 = NULL;
 705     rvec4            *f;
 706     matrix            box;
 707     gmx_bool          bPDB;
 708     int               isize;
 709     int              *index = NULL, *ind_fit = NULL;
 710     char             *grpname;
 711     t_cluster_ndx    *clust = NULL;
 712     t_dr_result       dr, *dr_clust = NULL;
 713     char            **leg;
 714     real             *vvindex = NULL, *w_rls = NULL;
 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 = NULL;
 721
 722     t_filenm          fnm[] = {
 723         { efTPR, NULL, NULL, ffREAD },
 724         { efTRX, "-f", NULL, 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, NULL,  "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, NULL, &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     read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE);
 752     snew(xtop, header.natoms);
 753     read_tpx(ftp2fn(efTPR, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, &mtop);
 754     bPDB = opt2bSet("-q", NFILE, fnm);
 755     if (bPDB)
 756     {
 757         snew(xav, ntopatoms);
 758         snew(ind_fit, ntopatoms);
 759         snew(w_rls, ntopatoms);
 760         for (kkk = 0; (kkk < ntopatoms); kkk++)
 761         {
 762             w_rls[kkk]   = 1;
 763             ind_fit[kkk] = kkk;
 764         }
 765
 766         snew(atoms, 1);
 767         *atoms = gmx_mtop_global_atoms(&mtop);
 768
 769         if (atoms->pdbinfo == NULL)
 770         {
 771             snew(atoms->pdbinfo, atoms->nr);
 772         }
 773     }
 774
 775     top = gmx_mtop_generate_local_top(&mtop, ir.efep != efepNO);
 776
 777     g        = NULL;
 778     pbc_null = NULL;
 779     if (ir.ePBC != epbcNONE)
 780     {
 781         if (ir.bPeriodicMols)
 782         {
 783             pbc_null = &pbc;
 784         }
 785         else
 786         {
 787             g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
 788         }
 789     }
 790
 791     if (ftp2bSet(efNDX, NFILE, fnm))
 792     {
 793         /* TODO: Nothing is written to this file if -c is provided, but it is
 794          * still opened... */
 795         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 796         xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
 797                        "nm", oenv);
 798         snew(vvindex, isize);
 799         snew(leg, isize);
 800         for (i = 0; (i < isize); i++)
 801         {
 802             index[i]++;
 803             snew(leg[i], 12);
 804             sprintf(leg[i], "index %d", index[i]);
 805         }
 806         xvgr_legend(xvg, isize, (const char**)leg, oenv);
 807     }
 808     else
 809     {
 810         isize = 0;
 811     }
 812
 813     ir.dr_tau = 0.0;
 814     init_disres(fplog, &mtop, &ir, NULL, &fcd, NULL, FALSE);
 815
 816     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 817     snew(f, 5*natoms);
 818
 819     init_dr_res(&dr, fcd.disres.nres);
 820     if (opt2bSet("-c", NFILE, fnm))
 821     {
 822         clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
 823         snew(dr_clust, clust->clust->nr+1);
 824         for (i = 0; (i <= clust->clust->nr); i++)
 825         {
 826             init_dr_res(&dr_clust[i], fcd.disres.nres);
 827         }
 828     }
 829     else
 830     {
 831         out = xvgropen(opt2fn("-ds", NFILE, fnm),
 832                        "Sum of Violations", "Time (ps)", "nm", oenv);
 833         aver = xvgropen(opt2fn("-da", NFILE, fnm),
 834                         "Average Violation", "Time (ps)", "nm", oenv);
 835         numv = xvgropen(opt2fn("-dn", NFILE, fnm),
 836                         "# Violations", "Time (ps)", "#", oenv);
 837         maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
 838                          "Largest Violation", "Time (ps)", "nm", oenv);
 839     }
 840
 841     mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
 842     atoms2md(&mtop, &ir, -1, NULL, mtop.natoms, mdatoms);
 843     update_mdatoms(mdatoms, ir.fepvals->init_lambda);
 844     init_nrnb(&nrnb);
 845     if (ir.ePBC != epbcNONE)
 846     {
 847         gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
 848     }
 849
 850     j = 0;
 851     do
 852     {
 853         if (ir.ePBC != epbcNONE)
 854         {
 855             if (ir.bPeriodicMols)
 856             {
 857                 set_pbc(&pbc, ir.ePBC, box);
 858             }
 859             else
 860             {
 861                 gmx_rmpbc(gpbc, natoms, box, x);
 862             }
 863         }
 864
 865         if (clust)
 866         {
 867             if (j > clust->maxframe)
 868             {
 869                 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
 870             }
 871             my_clust = clust->inv_clust[j];
 872             range_check(my_clust, 0, clust->clust->nr);
 873             check_viol(fplog, &(top->idef.il[F_DISRES]),
 874                        top->idef.iparams,
 875                        x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
 876         }
 877         else
 878         {
 879             check_viol(fplog, &(top->idef.il[F_DISRES]),
 880                        top->idef.iparams,
 881                        x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
 882         }
 883         if (bPDB)
 884         {
 885             reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
 886             do_fit(atoms->nr, w_rls, x, x);
 887             if (j == 0)
 888             {
 889                 /* Store the first frame of the trajectory as 'characteristic'
 890                  * for colouring with violations.
 891                  */
 892                 for (kkk = 0; (kkk < atoms->nr); kkk++)
 893                 {
 894                     copy_rvec(x[kkk], xav[kkk]);
 895                 }
 896             }
 897         }
 898         if (!clust)
 899         {
 900             if (isize > 0)
 901             {
 902                 fprintf(xvg, "%10g", t);
 903                 for (i = 0; (i < isize); i++)
 904                 {
 905                     fprintf(xvg, "  %10g", vvindex[i]);
 906                 }
 907                 fprintf(xvg, "\n");
 908             }
 909             fprintf(out,  "%10g  %10g\n", t, dr.sumv);
 910             fprintf(aver, "%10g  %10g\n", t, dr.averv);
 911             fprintf(maxxv, "%10g  %10g\n", t, dr.maxv);
 912             fprintf(numv, "%10g  %10d\n", t, dr.nv);
 913         }
 914         j++;
 915     }
 916     while (read_next_x(oenv, status, &t, x, box));
 917     close_trj(status);
 918     if (ir.ePBC != epbcNONE)
 919     {
 920         gmx_rmpbc_done(gpbc);
 921     }
 922
 923     if (clust)
 924     {
 925         dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 926                          top->idef.iparams, clust->clust, dr_clust,
 927                          clust->grpname, isize, index);
 928     }
 929     else
 930     {
 931         dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
 932                    top->idef.iparams, &dr, isize, index,
 933                    bPDB ? atoms : NULL);
 934         if (bPDB)
 935         {
 936             write_sto_conf(opt2fn("-q", NFILE, fnm),
 937                            "Coloured by average violation in Angstrom",
 938                            atoms, xav, NULL, ir.ePBC, box);
 939         }
 940         dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
 941                           j, &top->idef, &mtop, max_dr, nlevels, bThird);
 942         xvgrclose(out);
 943         xvgrclose(aver);
 944         xvgrclose(numv);
 945         xvgrclose(maxxv);
 946         do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
 947         do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
 948         do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
 949         do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
 950     }
 951     if (isize > 0)
 952     {
 953         xvgrclose(xvg);
 954         if (!clust)
 955         {
 956             do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
 957         }
 958     }
 959
 960     gmx_ffclose(fplog);
 961
 962     return 0;
 963 }