src/gromacs/gmxana/gmx_disre.c

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