2 * This file is part of the GROMACS molecular simulation package.
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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
45 #include <unordered_map>
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/commandline/viewit.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/fileio/matio.h"
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/fileio/tpxio.h"
53 #include "gromacs/fileio/trxio.h"
54 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/gmxana/gmx_ana.h"
56 #include "gromacs/gmxana/gstat.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/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/trajectoryanalysis/topologyinformation.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"
85 static t_toppop
* top
= nullptr;
91 real sumv
, averv
, maxv
;
92 real
*aver1
, *aver2
, *aver_3
, *aver_6
;
95 static void init5(int n
)
105 for (i
= 0; (i
< ntop
); i
++)
112 static void add5(int ndr
, real viol
)
117 for (i
= 1; (i
< ntop
); i
++)
119 if (top
[i
].v
< top
[mini
].v
)
124 if (viol
> top
[mini
].v
)
131 static void print5(FILE* fp
)
135 std::sort(top
, top
+ ntop
, [](const t_toppop
& a
, const t_toppop
& b
) { return a
.v
> b
.v
; }); // reverse sort
136 fprintf(fp
, "Index:");
137 for (i
= 0; (i
< ntop
); i
++)
139 fprintf(fp
, " %6d", top
[i
].n
);
141 fprintf(fp
, "\nViol: ");
142 for (i
= 0; (i
< ntop
); i
++)
144 fprintf(fp
, " %6.3f", top
[i
].v
);
149 static void check_viol(FILE* log
,
151 t_iparams forceparams
[],
164 int i
, j
, nat
, n
, type
, nviol
, ndr
, label
;
165 real rt
, mviol
, tviol
, viol
, lam
, dvdl
, drt
;
167 static gmx_bool bFirst
= TRUE
;
179 forceatoms
= disres
->iatoms
;
180 for (j
= 0; (j
< isize
); j
++)
184 nat
= interaction_function
[F_DISRES
].nratoms
+ 1;
185 for (i
= 0; (i
< disres
->nr
);)
187 type
= forceatoms
[i
];
189 label
= forceparams
[type
].disres
.label
;
192 fprintf(debug
, "DISRE: ndr = %d, label = %d i=%d, n =%d\n", ndr
, label
, i
, n
);
197 } while (((i
+ n
) < disres
->nr
) && (forceparams
[forceatoms
[i
+ n
]].disres
.label
== label
));
199 calc_disres_R_6(nullptr, nullptr, n
, &forceatoms
[i
], x
, pbc
, fcd
, nullptr);
201 if (fcd
->disres
.Rt_6
[label
] <= 0)
203 gmx_fatal(FARGS
, "ndr = %d, rt_6 = %f", ndr
, fcd
->disres
.Rt_6
[label
]);
206 rt
= gmx::invsixthroot(fcd
->disres
.Rt_6
[label
]);
207 dr
[clust_id
].aver1
[ndr
] += rt
;
208 dr
[clust_id
].aver2
[ndr
] += gmx::square(rt
);
209 drt
= 1.0 / gmx::power3(rt
);
210 dr
[clust_id
].aver_3
[ndr
] += drt
;
211 dr
[clust_id
].aver_6
[ndr
] += fcd
->disres
.Rt_6
[label
];
213 snew(fshift
, SHIFTS
);
214 ta_disres(n
, &forceatoms
[i
], forceparams
, x
, f
, fshift
, pbc
, g
, lam
, &dvdl
, nullptr, fcd
, nullptr);
216 viol
= fcd
->disres
.sumviol
;
223 add5(forceparams
[type
].disres
.label
, viol
);
230 for (j
= 0; (j
< isize
); j
++)
232 if (index
[j
] == forceparams
[type
].disres
.label
)
234 vvindex
[j
] = gmx::invsixthroot(fcd
->disres
.Rt_6
[label
]);
241 dr
[clust_id
].nv
= nviol
;
242 dr
[clust_id
].maxv
= mviol
;
243 dr
[clust_id
].sumv
= tviol
;
244 dr
[clust_id
].averv
= tviol
/ ndr
;
245 dr
[clust_id
].nframes
++;
249 fprintf(stderr
, "\nThere are %d restraints and %d pairs\n", ndr
, disres
->nr
/ nat
);
262 real up1
, r
, rT3
, rT6
, viol
, violT3
, violT6
;
265 static void dump_dump(FILE* log
, int ndr
, t_dr_stats drs
[])
267 static const char* core
[] = { "All restraints", "Core restraints" };
268 static const char* tp
[] = { "linear", "third power", "sixth power" };
269 real viol_tot
, viol_max
, viol
= 0;
274 for (int iCore
= 0; iCore
< 2; iCore
++)
276 bCore
= (iCore
== 1);
277 for (kkk
= 0; (kkk
< 3); kkk
++)
283 for (i
= 0; (i
< ndr
); i
++)
285 if (!bCore
|| drs
[i
].bCore
)
289 case 0: viol
= drs
[i
].viol
; break;
290 case 1: viol
= drs
[i
].violT3
; break;
291 case 2: viol
= drs
[i
].violT6
; break;
292 default: gmx_incons("Dumping violations");
294 viol_max
= std::max(viol_max
, viol
);
303 if ((nrestr
> 0) || (bCore
&& (nrestr
< ndr
)))
306 fprintf(log
, "+++++++ %s ++++++++\n", core
[bCore
]);
307 fprintf(log
, "+++++++ Using %s averaging: ++++++++\n", tp
[kkk
]);
308 fprintf(log
, "Sum of violations: %8.3f nm\n", viol_tot
);
311 fprintf(log
, "Average violation: %8.3f nm\n", viol_tot
/ nrestr
);
313 fprintf(log
, "Largest violation: %8.3f nm\n", viol_max
);
314 fprintf(log
, "Number of violated restraints: %d/%d\n", nviol
, nrestr
);
320 static void dump_viol(FILE* log
, int ndr
, t_dr_stats
* drs
, gmx_bool bLinear
)
324 fprintf(log
, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
325 for (i
= 0; (i
< ndr
); i
++)
327 if (bLinear
&& (drs
[i
].viol
== 0))
331 fprintf(log
, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n", drs
[i
].label
,
332 yesno_names
[drs
[i
].bCore
], drs
[i
].up1
, drs
[i
].r
, drs
[i
].rT3
, drs
[i
].rT6
,
333 drs
[i
].viol
, drs
[i
].violT3
, drs
[i
].violT6
);
337 static gmx_bool
is_core(int i
, int isize
, const int index
[])
340 gmx_bool bIC
= FALSE
;
342 for (kk
= 0; !bIC
&& (kk
< isize
); kk
++)
344 bIC
= (index
[kk
] == i
);
350 static void dump_stats(FILE* log
,
352 const t_disresdata
& dd
,
353 const t_ilist
* disres
,
363 fprintf(log
, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
365 const int nra
= interaction_function
[F_DISRES
].nratoms
+ 1;
366 for (int j
= 0; j
< disres
->nr
; j
+= nra
)
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");
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
));
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;
389 dump_viol(log
, dd
.nres
, drs
, FALSE
);
391 fprintf(log
, "+++ Sorted by linear averaged violations: +++\n");
392 std::sort(drs
, drs
+ dd
.nres
,
393 [](const t_dr_stats
& a
, const t_dr_stats
& b
) { return a
.viol
> b
.viol
; }); // Reverse sort
394 dump_viol(log
, dd
.nres
, drs
, TRUE
);
396 dump_dump(log
, dd
.nres
, drs
);
401 static void dump_clust_stats(FILE* fp
,
402 const t_disresdata
& dd
,
403 const t_ilist
* disres
,
412 double sumV
, maxV
, sumVT3
, sumVT6
, maxVT3
, maxVT6
;
416 fprintf(fp
, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
417 fprintf(fp
, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
421 for (k
= 0; (k
< clust
->nr
); k
++)
423 if (dr
[k
].nframes
== 0)
427 if (dr
[k
].nframes
!= (clust
->index
[k
+ 1] - clust
->index
[k
]))
430 "Inconsistency in cluster %s.\n"
431 "Found %d frames in trajectory rather than the expected %d\n",
432 clust_name
[k
], dr
[k
].nframes
, clust
->index
[k
+ 1] - clust
->index
[k
]);
436 gmx_fatal(FARGS
, "Inconsistency with cluster %d. Invalid name", k
);
438 nra
= interaction_function
[F_DISRES
].nratoms
+ 1;
439 sumV
= sumVT3
= sumVT6
= maxV
= maxVT3
= maxVT6
= 0;
441 // Use a map to process each restraint only once while looping over all pairs
442 std::unordered_map
<int, bool> restraintHasBeenProcessed
;
443 for (int j
= 0; j
< dd
.nres
; j
+= nra
)
445 // Note that the restraint i can be used by multiple pairs
446 const int i
= disres
->iatoms
[j
] - dd
.type_min
;
448 if (restraintHasBeenProcessed
[i
])
453 drs
[i
].label
= ip
[disres
->iatoms
[j
]].disres
.label
;
454 drs
[i
].bCore
= is_core(drs
[i
].label
, isize
, index
);
455 drs
[i
].up1
= ip
[disres
->iatoms
[j
]].disres
.up1
;
456 drs
[i
].r
= dr
[k
].aver1
[i
] / dr
[k
].nframes
;
457 if ((dr
[k
].aver_3
[i
] <= 0) || !std::isfinite(dr
[k
].aver_3
[i
]))
459 gmx_fatal(FARGS
, "dr[%d].aver_3[%d] = %f", k
, i
, dr
[k
].aver_3
[i
]);
461 drs
[i
].rT3
= gmx::invcbrt(dr
[k
].aver_3
[i
] / dr
[k
].nframes
);
462 drs
[i
].rT6
= gmx::invsixthroot(dr
[k
].aver_6
[i
] / dr
[k
].nframes
);
463 drs
[i
].viol
= std::max(0.0, static_cast<double>(drs
[i
].r
- drs
[i
].up1
));
464 drs
[i
].violT3
= std::max(0.0, static_cast<double>(drs
[i
].rT3
- drs
[i
].up1
));
465 drs
[i
].violT6
= std::max(0.0, static_cast<double>(drs
[i
].rT6
- drs
[i
].up1
));
467 sumVT3
+= drs
[i
].violT3
;
468 sumVT6
+= drs
[i
].violT6
;
469 maxV
= std::max(maxV
, static_cast<double>(drs
[i
].viol
));
470 maxVT3
= std::max(maxVT3
, static_cast<double>(drs
[i
].violT3
));
471 maxVT6
= std::max(maxVT6
, static_cast<double>(drs
[i
].violT6
));
473 // We have processed restraint i, mark it as such
474 restraintHasBeenProcessed
[i
] = true;
476 if (std::strcmp(clust_name
[k
], "1000") == 0)
480 fprintf(fp
, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n", clust_name
[k
],
481 dr
[k
].nframes
, sumV
, maxV
, sumVT3
, maxVT3
, sumVT6
, maxVT6
);
487 static void init_dr_res(t_dr_result
* dr
, int ndr
)
489 snew(dr
->aver1
, ndr
+ 1);
490 snew(dr
->aver2
, ndr
+ 1);
491 snew(dr
->aver_3
, ndr
+ 1);
492 snew(dr
->aver_6
, ndr
+ 1);
500 static void dump_disre_matrix(const char* fn
,
505 const gmx_mtop_t
* mtop
,
512 int n_res
, a_offset
, mol
, a
;
513 int i
, j
, nra
, nratoms
, tp
, ri
, rj
, index
, nlabel
, label
;
515 real
**matrix
, *t_res
, hi
, *w_dr
, rav
, rviol
;
516 t_rgb rlo
= { 1, 1, 1 };
517 t_rgb rhi
= { 0, 0, 0 };
523 snew(resnr
, mtop
->natoms
);
526 for (const gmx_molblock_t
& molb
: mtop
->molblock
)
528 const t_atoms
& atoms
= mtop
->moltype
[molb
.type
].atoms
;
529 for (mol
= 0; mol
< molb
.nmol
; mol
++)
531 for (a
= 0; a
< atoms
.nr
; a
++)
533 resnr
[a_offset
+ a
] = n_res
+ atoms
.atom
[a
].resind
;
536 a_offset
+= atoms
.nr
;
541 for (i
= 0; (i
< n_res
); i
++)
546 for (i
= 0; (i
< n_res
); i
++)
548 snew(matrix
[i
], n_res
);
550 nratoms
= interaction_function
[F_DISRES
].nratoms
;
551 nra
= (idef
->il
[F_DISRES
].nr
/ (nratoms
+ 1));
557 for (i
= 0; (i
< idef
->il
[F_DISRES
].nr
); i
+= nratoms
+ 1)
559 tp
= idef
->il
[F_DISRES
].iatoms
[i
];
560 label
= idef
->iparams
[tp
].disres
.label
;
564 /* Set index pointer */
568 gmx_fatal(FARGS
, "nlabel is %d, label = %d", nlabel
, label
);
572 gmx_fatal(FARGS
, "ndr = %d, index = %d", ndr
, index
);
574 /* Update the weight */
575 w_dr
[index
] = 1.0 / nlabel
;
584 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel
, index
, ndr
);
586 for (i
= 0; (i
< ndr
); i
++)
588 for (j
= ptr
[i
]; (j
< ptr
[i
+ 1]); j
+= nratoms
+ 1)
590 tp
= idef
->il
[F_DISRES
].iatoms
[j
];
591 ai
= idef
->il
[F_DISRES
].iatoms
[j
+ 1];
592 aj
= idef
->il
[F_DISRES
].iatoms
[j
+ 2];
598 rav
= gmx::invcbrt(dr
->aver_3
[i
] / nsteps
);
602 rav
= dr
->aver1
[i
] / nsteps
;
606 fprintf(debug
, "DR %d, atoms %d, %d, distance %g\n", i
, ai
, aj
, rav
);
608 rviol
= std::max(0.0_real
, rav
- idef
->iparams
[tp
].disres
.up1
);
609 matrix
[ri
][rj
] += w_dr
[i
] * rviol
;
610 matrix
[rj
][ri
] += w_dr
[i
] * rviol
;
611 hi
= std::max(hi
, matrix
[ri
][rj
]);
612 hi
= std::max(hi
, matrix
[rj
][ri
]);
622 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest "
623 "value in your simulation (%g)\n",
628 printf("Highest level in the matrix will be %g\n", hi
);
629 fp
= gmx_ffopen(fn
, "w");
630 write_xpm(fp
, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue", n_res
, n_res
, t_res
,
631 t_res
, matrix
, 0, hi
, rlo
, rhi
, &nlevels
);
635 int gmx_disre(int argc
, char* argv
[])
637 const char* desc
[] = {
638 "[THISMODULE] computes violations of distance restraints.",
639 "The program always",
640 "computes the instantaneous violations rather than time-averaged,",
641 "because this analysis is done from a trajectory file afterwards",
642 "it does not make sense to use time averaging. However,",
643 "the time averaged values per restraint are given in the log file.[PAR]",
644 "An index file may be used to select specific restraints by index group label for",
646 "When the optional [TT]-q[tt] flag is given a [REF].pdb[ref] file coloured by the",
647 "amount of average violations.[PAR]",
648 "When the [TT]-c[tt] option is given, an index file will be read",
649 "containing the frames in your trajectory corresponding to the clusters",
650 "(defined in another manner) that you want to analyze. For these clusters",
651 "the program will compute average violations using the third power",
652 "averaging algorithm and print them in the log file."
655 static int nlevels
= 20;
656 static real max_dr
= 0;
657 static gmx_bool bThird
= TRUE
;
663 "Number of large violations that are stored in the log file every step" },
668 "Maximum distance violation in matrix output. If less than or equal to 0 the "
669 "maximum will be determined by the data." },
670 { "-nlevels", FALSE
, etINT
, { &nlevels
}, "Number of levels in the matrix output" },
675 "Use inverse third power averaging or linear for matrix output" }
678 FILE *out
= nullptr, *aver
= nullptr, *numv
= nullptr, *maxxv
= nullptr, *xvg
= nullptr;
685 rvec
* x
, *xav
= nullptr;
690 int * index
= nullptr, *ind_fit
= nullptr;
692 t_cluster_ndx
* clust
= nullptr;
693 t_dr_result dr
, *dr_clust
= nullptr;
695 real
* vvindex
= nullptr, *w_rls
= nullptr;
696 t_pbc pbc
, *pbc_null
;
699 gmx_output_env_t
* oenv
;
700 gmx_rmpbc_t gpbc
= nullptr;
702 t_filenm fnm
[] = { { efTPR
, nullptr, nullptr, ffREAD
}, { efTRX
, "-f", nullptr, ffREAD
},
703 { efXVG
, "-ds", "drsum", ffWRITE
}, { efXVG
, "-da", "draver", ffWRITE
},
704 { efXVG
, "-dn", "drnum", ffWRITE
}, { efXVG
, "-dm", "drmax", ffWRITE
},
705 { efXVG
, "-dr", "restr", ffWRITE
}, { efLOG
, "-l", "disres", ffWRITE
},
706 { efNDX
, nullptr, "viol", ffOPTRD
}, { efPDB
, "-q", "viol", ffOPTWR
},
707 { efNDX
, "-c", "clust", ffOPTRD
}, { efXPM
, "-x", "matrix", ffOPTWR
} };
708 #define NFILE asize(fnm)
710 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_CAN_VIEW
, NFILE
, fnm
, asize(pa
), pa
,
711 asize(desc
), desc
, 0, nullptr, &oenv
))
716 fplog
= ftp2FILE(efLOG
, NFILE
, fnm
, "w");
723 t_inputrec irInstance
;
724 t_inputrec
* ir
= &irInstance
;
726 gmx::TopologyInformation topInfo
;
727 topInfo
.fillFromInputFile(ftp2fn(efTPR
, NFILE
, fnm
));
728 int ntopatoms
= topInfo
.mtop()->natoms
;
730 bPDB
= opt2bSet("-q", NFILE
, fnm
);
733 snew(xav
, ntopatoms
);
734 snew(ind_fit
, ntopatoms
);
735 snew(w_rls
, ntopatoms
);
736 for (kkk
= 0; (kkk
< ntopatoms
); kkk
++)
742 atoms
= topInfo
.copyAtoms();
744 if (atoms
->pdbinfo
== nullptr)
746 snew(atoms
->pdbinfo
, atoms
->nr
);
748 atoms
->havePdbInfo
= TRUE
;
751 gmx_mtop_generate_local_top(*topInfo
.mtop(), &top
, ir
->efep
!= efepNO
);
755 if (ir
->ePBC
!= epbcNONE
)
757 if (ir
->bPeriodicMols
)
763 g
= mk_graph(fplog
, &top
.idef
, 0, ntopatoms
, FALSE
, FALSE
);
767 if (ftp2bSet(efNDX
, NFILE
, fnm
))
769 /* TODO: Nothing is written to this file if -c is provided, but it is
771 rd_index(ftp2fn(efNDX
, NFILE
, fnm
), 1, &isize
, &index
, &grpname
);
772 xvg
= xvgropen(opt2fn("-dr", NFILE
, fnm
), "Individual Restraints", "Time (ps)", "nm", oenv
);
773 snew(vvindex
, isize
);
775 for (i
= 0; (i
< isize
); i
++)
779 sprintf(leg
[i
], "index %d", index
[i
]);
781 xvgr_legend(xvg
, isize
, leg
, oenv
);
789 init_disres(fplog
, topInfo
.mtop(), ir
, nullptr, nullptr, &fcd
, nullptr, FALSE
);
791 int natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
794 init_dr_res(&dr
, fcd
.disres
.nres
);
795 if (opt2bSet("-c", NFILE
, fnm
))
797 clust
= cluster_index(fplog
, opt2fn("-c", NFILE
, fnm
));
798 snew(dr_clust
, clust
->clust
->nr
+ 1);
799 for (i
= 0; (i
<= clust
->clust
->nr
); i
++)
801 init_dr_res(&dr_clust
[i
], fcd
.disres
.nres
);
806 out
= xvgropen(opt2fn("-ds", NFILE
, fnm
), "Sum of Violations", "Time (ps)", "nm", oenv
);
807 aver
= xvgropen(opt2fn("-da", NFILE
, fnm
), "Average Violation", "Time (ps)", "nm", oenv
);
808 numv
= xvgropen(opt2fn("-dn", NFILE
, fnm
), "# Violations", "Time (ps)", "#", oenv
);
809 maxxv
= xvgropen(opt2fn("-dm", NFILE
, fnm
), "Largest Violation", "Time (ps)", "nm", oenv
);
812 auto mdAtoms
= gmx::makeMDAtoms(fplog
, *topInfo
.mtop(), *ir
, false);
813 atoms2md(topInfo
.mtop(), ir
, -1, nullptr, ntopatoms
, mdAtoms
.get());
814 update_mdatoms(mdAtoms
->mdatoms(), ir
->fepvals
->init_lambda
);
815 if (ir
->ePBC
!= epbcNONE
)
817 gpbc
= gmx_rmpbc_init(&top
.idef
, ir
->ePBC
, natoms
);
823 if (ir
->ePBC
!= epbcNONE
)
825 if (ir
->bPeriodicMols
)
827 set_pbc(&pbc
, ir
->ePBC
, box
);
831 gmx_rmpbc(gpbc
, natoms
, box
, x
);
837 if (j
> clust
->maxframe
)
840 "There are more frames in the trajectory than in the cluster index file. "
844 my_clust
= clust
->inv_clust
[j
];
845 range_check(my_clust
, 0, clust
->clust
->nr
);
846 check_viol(fplog
, &(top
.idef
.il
[F_DISRES
]), top
.idef
.iparams
, x
, f
, pbc_null
, g
,
847 dr_clust
, my_clust
, isize
, index
, vvindex
, &fcd
);
851 check_viol(fplog
, &(top
.idef
.il
[F_DISRES
]), top
.idef
.iparams
, x
, f
, pbc_null
, g
, &dr
, 0,
852 isize
, index
, vvindex
, &fcd
);
856 reset_x(atoms
->nr
, ind_fit
, atoms
->nr
, nullptr, x
, w_rls
);
857 do_fit(atoms
->nr
, w_rls
, x
, x
);
860 /* Store the first frame of the trajectory as 'characteristic'
861 * for colouring with violations.
863 for (kkk
= 0; (kkk
< atoms
->nr
); kkk
++)
865 copy_rvec(x
[kkk
], xav
[kkk
]);
873 fprintf(xvg
, "%10g", t
);
874 for (i
= 0; (i
< isize
); i
++)
876 fprintf(xvg
, " %10g", vvindex
[i
]);
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
);
886 } while (read_next_x(oenv
, status
, &t
, x
, box
));
888 if (ir
->ePBC
!= epbcNONE
)
890 gmx_rmpbc_done(gpbc
);
895 dump_clust_stats(fplog
, fcd
.disres
, &(top
.idef
.il
[F_DISRES
]), top
.idef
.iparams
,
896 clust
->clust
, dr_clust
, clust
->grpname
, isize
, index
);
900 dump_stats(fplog
, j
, fcd
.disres
, &(top
.idef
.il
[F_DISRES
]), top
.idef
.iparams
, &dr
, isize
,
901 index
, bPDB
? atoms
.get() : nullptr);
904 write_sto_conf(opt2fn("-q", NFILE
, fnm
), "Coloured by average violation in Angstrom",
905 atoms
.get(), xav
, nullptr, ir
->ePBC
, box
);
907 dump_disre_matrix(opt2fn_null("-x", NFILE
, fnm
), &dr
, fcd
.disres
.nres
, j
, &top
.idef
,
908 topInfo
.mtop(), max_dr
, nlevels
, bThird
);
913 do_view(oenv
, opt2fn("-dn", NFILE
, fnm
), "-nxy");
914 do_view(oenv
, opt2fn("-da", NFILE
, fnm
), "-nxy");
915 do_view(oenv
, opt2fn("-ds", NFILE
, fnm
), "-nxy");
916 do_view(oenv
, opt2fn("-dm", NFILE
, fnm
), "-nxy");
923 do_view(oenv
, opt2fn("-dr", NFILE
, fnm
), "-nxy");