src/gromacs/gmxana/gmx_disre.c

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