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