?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Introduce SimulatorBuilder
[gromacs.git] / src / gromacs / gmxana / gmx_disre.cpp
blob9c6fe6564fbaf32d275c7f4d015389c5466a302e
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
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/listed_forces/disre.h"
56 #include "gromacs/math/do_fit.h"
57 #include "gromacs/math/functions.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/force.h"
60 #include "gromacs/mdlib/mdatoms.h"
61 #include "gromacs/mdtypes/fcdata.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/pbcutil/ishift.h"
65 #include "gromacs/pbcutil/mshift.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/pbcutil/rmpbc.h"
68 #include "gromacs/topology/index.h"
69 #include "gromacs/topology/mtop_util.h"
70 #include "gromacs/topology/topology.h"
71 #include "gromacs/trajectoryanalysis/topologyinformation.h"
72 #include "gromacs/utility/arraysize.h"
73 #include "gromacs/utility/fatalerror.h"
74 #include "gromacs/utility/futil.h"
75 #include "gromacs/utility/smalloc.h"
76
77 typedef struct {
78 int n;
79 real v;
80 } t_toppop;
81
82 static t_toppop *top = nullptr;
83 static int ntop = 0;
84
85 typedef struct {
86 int nv, nframes;
87 real sumv, averv, maxv;
88 real *aver1, *aver2, *aver_3, *aver_6;
89 } t_dr_result;
90
91 static void init5(int n)
92 {
93 ntop = n;
94 snew(top, ntop);
95 }
96
97 static void reset5()
98 {
99 int i;
100
101 for (i = 0; (i < ntop); i++)
102 {
103 top[i].n = -1;
104 top[i].v = 0;
105 }
106 }
107
108 static void add5(int ndr, real viol)
109 {
110 int i, mini;
111
112 mini = 0;
113 for (i = 1; (i < ntop); i++)
114 {
115 if (top[i].v < top[mini].v)
116 {
117 mini = i;
118 }
119 }
120 if (viol > top[mini].v)
121 {
122 top[mini].v = viol;
123 top[mini].n = ndr;
124 }
125 }
126
127 static void print5(FILE *fp)
128 {
129 int i;
130
131 std::sort(top, top+ntop, [](const t_toppop &a, const t_toppop &b) {return a.v > b.v; }); //reverse sort
132 fprintf(fp, "Index:");
133 for (i = 0; (i < ntop); i++)
134 {
135 fprintf(fp, " %6d", top[i].n);
136 }
137 fprintf(fp, "\nViol: ");
138 for (i = 0; (i < ntop); i++)
139 {
140 fprintf(fp, " %6.3f", top[i].v);
141 }
142 fprintf(fp, "\n");
143 }
144
145 static void check_viol(FILE *log,
146 t_ilist *disres, t_iparams forceparams[],
147 rvec x[], rvec4 f[],
148 t_pbc *pbc, t_graph *g, t_dr_result dr[],
149 int clust_id, int isize, const int index[], real vvindex[],
150 t_fcdata *fcd)
151 {
152 t_iatom *forceatoms;
153 int i, j, nat, n, type, nviol, ndr, label;
154 real rt, mviol, tviol, viol, lam, dvdl, drt;
155 rvec *fshift;
156 static gmx_bool bFirst = TRUE;
157
158 lam = 0;
159 dvdl = 0;
160 tviol = 0;
161 nviol = 0;
162 mviol = 0;
163 ndr = 0;
164 if (ntop)
165 {
166 reset5();
167 }
168 forceatoms = disres->iatoms;
169 for (j = 0; (j < isize); j++)
170 {
171 vvindex[j] = 0;
172 }
173 nat = interaction_function[F_DISRES].nratoms+1;
174 for (i = 0; (i < disres->nr); )
175 {
176 type = forceatoms[i];
177 n = 0;
178 label = forceparams[type].disres.label;
179 if (debug)
180 {
181 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
182 ndr, label, i, n);
183 }
184 if (ndr != label)
185 {
186 gmx_fatal(FARGS, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
187 }
188 do
189 {
190 n += nat;
191 }
192 while (((i+n) < disres->nr) &&
193 (forceparams[forceatoms[i+n]].disres.label == label));
194
195 calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i],
196 x, pbc, fcd, nullptr);
197
198 if (fcd->disres.Rt_6[label] <= 0)
199 {
200 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[label]);
201 }
202
203 rt = gmx::invsixthroot(fcd->disres.Rt_6[label]);
204 dr[clust_id].aver1[ndr] += rt;
205 dr[clust_id].aver2[ndr] += gmx::square(rt);
206 drt = 1.0/gmx::power3(rt);
207 dr[clust_id].aver_3[ndr] += drt;
208 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label];
209
210 snew(fshift, SHIFTS);
211 ta_disres(n, &forceatoms[i], forceparams,
212 x, f, fshift,
213 pbc, g, lam, &dvdl, nullptr, fcd, nullptr);
214 sfree(fshift);
215 viol = fcd->disres.sumviol;
216
217 if (viol > 0)
218 {
219 nviol++;
220 if (ntop)
221 {
222 add5(forceparams[type].disres.label, viol);
223 }
224 if (viol > mviol)
225 {
226 mviol = viol;
227 }
228 tviol += viol;
229 for (j = 0; (j < isize); j++)
230 {
231 if (index[j] == forceparams[type].disres.label)
232 {
233 vvindex[j] = gmx::invsixthroot(fcd->disres.Rt_6[label]);
234 }
235 }
236 }
237 ndr++;
238 i += n;
239 }
240 dr[clust_id].nv = nviol;
241 dr[clust_id].maxv = mviol;
242 dr[clust_id].sumv = tviol;
243 dr[clust_id].averv = tviol/ndr;
244 dr[clust_id].nframes++;
245
246 if (bFirst)
247 {
248 fprintf(stderr, "\nThere are %d restraints and %d pairs\n",
249 ndr, disres->nr/nat);
250 bFirst = FALSE;
251 }
252 if (ntop)
253 {
254 print5(log);
255 }
256 }
257
258 typedef struct {
259 int index;
260 gmx_bool bCore;
261 real up1, r, rT3, rT6, viol, violT3, violT6;
262 } t_dr_stats;
263
264 static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
265 {
266 static const char *core[] = { "All restraints", "Core restraints" };
267 static const char *tp[] = { "linear", "third power", "sixth power" };
268 real viol_tot, viol_max, viol = 0;
269 gmx_bool bCore;
270 int nviol, nrestr;
271 int i, kkk;
272
273 for (int iCore = 0; iCore < 2; iCore++)
274 {
275 bCore = (iCore == 1);
276 for (kkk = 0; (kkk < 3); kkk++)
277 {
278 viol_tot = 0;
279 viol_max = 0;
280 nviol = 0;
281 nrestr = 0;
282 for (i = 0; (i < ndr); i++)
283 {
284 if (!bCore || drs[i].bCore)
285 {
286 switch (kkk)
287 {
288 case 0:
289 viol = drs[i].viol;
290 break;
291 case 1:
292 viol = drs[i].violT3;
293 break;
294 case 2:
295 viol = drs[i].violT6;
296 break;
297 default:
298 gmx_incons("Dumping violations");
299 }
300 viol_max = std::max(viol_max, viol);
301 if (viol > 0)
302 {
303 nviol++;
304 }
305 viol_tot += viol;
306 nrestr++;
307 }
308 }
309 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
310 {
311 fprintf(log, "\n");
312 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
313 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
314 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
315 if (nrestr > 0)
316 {
317 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
318 }
319 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
320 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
321 }
322 }
323 }
324 }
325
326 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
327 {
328 int i;
329
330 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
331 for (i = 0; (i < ndr); i++)
332 {
333 if (bLinear && (drs[i].viol == 0))
334 {
335 break;
336 }
337 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
338 drs[i].index, yesno_names[drs[i].bCore],
339 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
340 drs[i].viol, drs[i].violT3, drs[i].violT6);
341 }
342 }
343
344 static gmx_bool is_core(int i, int isize, const int index[])
345 {
346 int kk;
347 gmx_bool bIC = FALSE;
348
349 for (kk = 0; !bIC && (kk < isize); kk++)
350 {
351 bIC = (index[kk] == i);
352 }
353
354 return bIC;
355 }
356
357 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
358 t_iparams ip[], t_dr_result *dr,
359 int isize, int index[], t_atoms *atoms)
360 {
361 int i, j, nra;
362 t_dr_stats *drs;
363
364 fprintf(log, "\n");
365 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
366 snew(drs, ndr);
367 j = 0;
368 nra = interaction_function[F_DISRES].nratoms+1;
369 for (i = 0; (i < ndr); i++)
370 {
371 /* Search for the appropriate restraint */
372 for (; (j < disres->nr); j += nra)
373 {
374 if (ip[disres->iatoms[j]].disres.label == i)
375 {
376 break;
377 }
378 }
379 drs[i].index = i;
380 drs[i].bCore = is_core(i, isize, index);
381 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
382 drs[i].r = dr->aver1[i]/nsteps;
383 drs[i].rT3 = gmx::invcbrt(dr->aver_3[i]/nsteps);
384 drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i]/nsteps);
385 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
386 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
387 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
388 if (atoms)
389 {
390 int j1 = disres->iatoms[j+1];
391 int j2 = disres->iatoms[j+2];
392 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
393 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
394 }
395 }
396 dump_viol(log, ndr, drs, FALSE);
397
398 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
399 std::sort(drs, drs+ndr, [](const t_dr_stats &a, const t_dr_stats &b)
400 {return a.viol > b.viol; }); //Reverse sort
401 dump_viol(log, ndr, drs, TRUE);
402
403 dump_dump(log, ndr, drs);
404
405 sfree(drs);
406 }
407
408 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
409 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
410 char *clust_name[], int isize, int index[])
411 {
412 int i, j, k, nra, mmm = 0;
413 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
414 t_dr_stats *drs;
415
416 fprintf(fp, "\n");
417 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
418 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
419
420 snew(drs, ndr);
421
422 for (k = 0; (k < clust->nr); k++)
423 {
424 if (dr[k].nframes == 0)
425 {
426 continue;
427 }
428 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
429 {
430 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
431 "Found %d frames in trajectory rather than the expected %d\n",
432 clust_name[k], dr[k].nframes,
433 clust->index[k+1]-clust->index[k]);
434 }
435 if (!clust_name[k])
436 {
437 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
438 }
439 j = 0;
440 nra = interaction_function[F_DISRES].nratoms+1;
441 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
442 for (i = 0; (i < ndr); i++)
443 {
444 /* Search for the appropriate restraint */
445 for (; (j < disres->nr); j += nra)
446 {
447 if (ip[disres->iatoms[j]].disres.label == i)
448 {
449 break;
450 }
451 }
452 drs[i].index = i;
453 drs[i].bCore = is_core(i, isize, index);
454 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
455 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
456 if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
457 {
458 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
459 }
460 drs[i].rT3 = gmx::invcbrt(dr[k].aver_3[i]/dr[k].nframes);
461 drs[i].rT6 = gmx::invsixthroot(dr[k].aver_6[i]/dr[k].nframes);
462 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
463 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
464 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
465 sumV += drs[i].viol;
466 sumVT3 += drs[i].violT3;
467 sumVT6 += drs[i].violT6;
468 maxV = std::max(maxV, static_cast<double>(drs[i].viol));
469 maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
470 maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
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 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.nres, &(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.nres, &(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