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