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