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