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