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