?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Made distance restraints work with threads and DD
[gromacs.git] / src / gromacs / gmxana / gmx_disre.cpp
blob94dc678cda99126a2d2c54e5054c361899d2be17
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, 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 t_toppop *top = NULL;
84 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 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, 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(NULL, n, &forceatoms[i],
208 (const rvec*)x, pbc, fcd, NULL);
209
210 if (fcd->disres.Rt_6[0] <= 0)
211 {
212 gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
213 }
214
215 rt = gmx::invsixthroot(fcd->disres.Rt_6[0]);
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[0];
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, NULL, fcd, NULL);
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[0]);
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 (bCore = FALSE; (bCore <= TRUE); bCore++)
307 {
308 for (kkk = 0; (kkk < 3); kkk++)
309 {
310 viol_tot = 0;
311 viol_max = 0;
312 nviol = 0;
313 nrestr = 0;
314 for (i = 0; (i < ndr); i++)
315 {
316 if (!bCore || drs[i].bCore)
317 {
318 switch (kkk)
319 {
320 case 0:
321 viol = drs[i].viol;
322 break;
323 case 1:
324 viol = drs[i].violT3;
325 break;
326 case 2:
327 viol = drs[i].violT6;
328 break;
329 default:
330 gmx_incons("Dumping violations");
331 }
332 viol_max = std::max(viol_max, viol);
333 if (viol > 0)
334 {
335 nviol++;
336 }
337 viol_tot += viol;
338 nrestr++;
339 }
340 }
341 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
342 {
343 fprintf(log, "\n");
344 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
345 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
346 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
347 if (nrestr > 0)
348 {
349 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
350 }
351 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
352 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
353 }
354 }
355 }
356 }
357
358 static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
359 {
360 int i;
361
362 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
363 for (i = 0; (i < ndr); i++)
364 {
365 if (bLinear && (drs[i].viol == 0))
366 {
367 break;
368 }
369 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
370 drs[i].index, yesno_names[drs[i].bCore],
371 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
372 drs[i].viol, drs[i].violT3, drs[i].violT6);
373 }
374 }
375
376 static gmx_bool is_core(int i, int isize, int index[])
377 {
378 int kk;
379 gmx_bool bIC = FALSE;
380
381 for (kk = 0; !bIC && (kk < isize); kk++)
382 {
383 bIC = (index[kk] == i);
384 }
385
386 return bIC;
387 }
388
389 static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
390 t_iparams ip[], t_dr_result *dr,
391 int isize, int index[], t_atoms *atoms)
392 {
393 int i, j, nra;
394 t_dr_stats *drs;
395
396 fprintf(log, "\n");
397 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
398 snew(drs, ndr);
399 j = 0;
400 nra = interaction_function[F_DISRES].nratoms+1;
401 for (i = 0; (i < ndr); i++)
402 {
403 /* Search for the appropriate restraint */
404 for (; (j < disres->nr); j += nra)
405 {
406 if (ip[disres->iatoms[j]].disres.label == i)
407 {
408 break;
409 }
410 }
411 drs[i].index = i;
412 drs[i].bCore = is_core(i, isize, index);
413 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
414 drs[i].r = dr->aver1[i]/nsteps;
415 drs[i].rT3 = gmx::invcbrt(dr->aver_3[i]/nsteps);
416 drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i]/nsteps);
417 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
418 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
419 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
420 if (atoms)
421 {
422 int j1 = disres->iatoms[j+1];
423 int j2 = disres->iatoms[j+2];
424 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
425 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
426 }
427 }
428 dump_viol(log, ndr, drs, FALSE);
429
430 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
431 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
432 dump_viol(log, ndr, drs, TRUE);
433
434 dump_dump(log, ndr, drs);
435
436 sfree(drs);
437 }
438
439 static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
440 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
441 char *clust_name[], int isize, int index[])
442 {
443 int i, j, k, nra, mmm = 0;
444 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
445 t_dr_stats *drs;
446
447 fprintf(fp, "\n");
448 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
449 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
450
451 snew(drs, ndr);
452
453 for (k = 0; (k < clust->nr); k++)
454 {
455 if (dr[k].nframes == 0)
456 {
457 continue;
458 }
459 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
460 {
461 gmx_fatal(FARGS, "Inconsistency in cluster %s.\n"
462 "Found %d frames in trajectory rather than the expected %d\n",
463 clust_name[k], dr[k].nframes,
464 clust->index[k+1]-clust->index[k]);
465 }
466 if (!clust_name[k])
467 {
468 gmx_fatal(FARGS, "Inconsistency with cluster %d. Invalid name", k);
469 }
470 j = 0;
471 nra = interaction_function[F_DISRES].nratoms+1;
472 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
473 for (i = 0; (i < ndr); i++)
474 {
475 /* Search for the appropriate restraint */
476 for (; (j < disres->nr); j += nra)
477 {
478 if (ip[disres->iatoms[j]].disres.label == i)
479 {
480 break;
481 }
482 }
483 drs[i].index = i;
484 drs[i].bCore = is_core(i, isize, index);
485 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
486 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
487 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
488 {
489 gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
490 }
491 drs[i].rT3 = gmx::invcbrt(dr[k].aver_3[i]/dr[k].nframes);
492 drs[i].rT6 = gmx::invsixthroot(dr[k].aver_6[i]/dr[k].nframes);
493 drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
494 drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
495 drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
496 sumV += drs[i].viol;
497 sumVT3 += drs[i].violT3;
498 sumVT6 += drs[i].violT6;
499 maxV = std::max(maxV, static_cast<double>(drs[i].viol));
500 maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
501 maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
502 }
503 if (std::strcmp(clust_name[k], "1000") == 0)
504 {
505 mmm++;
506 }
507 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
508 clust_name[k],
509 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
510
511 }
512 fflush(fp);
513 sfree(drs);
514 }
515
516 static void init_dr_res(t_dr_result *dr, int ndr)
517 {
518 snew(dr->aver1, ndr+1);
519 snew(dr->aver2, ndr+1);
520 snew(dr->aver_3, ndr+1);
521 snew(dr->aver_6, ndr+1);
522 dr->nv = 0;
523 dr->nframes = 0;
524 dr->sumv = 0;
525 dr->maxv = 0;
526 dr->averv = 0;
527 }
528
529 static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
530 int nsteps, t_idef *idef, const gmx_mtop_t *mtop,
531 real max_dr, int nlevels, gmx_bool bThird)
532 {
533 FILE *fp;
534 int *resnr;
535 int n_res, a_offset, mb, mol, a;
536 t_atoms *atoms;
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 == NULL)
543 {
544 return;
545 }
546
547 snew(resnr, mtop->natoms);
548 n_res = 0;
549 a_offset = 0;
550 for (mb = 0; mb < mtop->nmolblock; mb++)
551 {
552 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
553 for (mol = 0; mol < mtop->molblock[mb].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 = NULL, *aver = NULL, *numv = NULL, *maxxv = NULL, *xvg = NULL;
692 t_tpxheader header;
693 t_inputrec ir;
694 gmx_mtop_t mtop;
695 rvec *xtop;
696 gmx_localtop_t *top;
697 t_atoms *atoms = NULL;
698 t_fcdata fcd;
699 t_nrnb nrnb;
700 t_graph *g;
701 int ntopatoms, natoms, i, j, kkk;
702 t_trxstatus *status;
703 real t;
704 rvec *x, *xav = NULL;
705 rvec4 *f;
706 matrix box;
707 gmx_bool bPDB;
708 int isize;
709 int *index = NULL, *ind_fit = NULL;
710 char *grpname;
711 t_cluster_ndx *clust = NULL;
712 t_dr_result dr, *dr_clust = NULL;
713 char **leg;
714 real *vvindex = NULL, *w_rls = NULL;
715 t_mdatoms *mdatoms;
716 t_pbc pbc, *pbc_null;
717 int my_clust;
718 FILE *fplog;
719 gmx_output_env_t *oenv;
720 gmx_rmpbc_t gpbc = NULL;
721
722 t_filenm fnm[] = {
723 { efTPR, NULL, NULL, ffREAD },
724 { efTRX, "-f", NULL, ffREAD },
725 { efXVG, "-ds", "drsum", ffWRITE },
726 { efXVG, "-da", "draver", ffWRITE },
727 { efXVG, "-dn", "drnum", ffWRITE },
728 { efXVG, "-dm", "drmax", ffWRITE },
729 { efXVG, "-dr", "restr", ffWRITE },
730 { efLOG, "-l", "disres", ffWRITE },
731 { efNDX, NULL, "viol", ffOPTRD },
732 { efPDB, "-q", "viol", ffOPTWR },
733 { efNDX, "-c", "clust", ffOPTRD },
734 { efXPM, "-x", "matrix", ffOPTWR }
735 };
736 #define NFILE asize(fnm)
737
738 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
739 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
740 {
741 return 0;
742 }
743
744 fplog = ftp2FILE(efLOG, NFILE, fnm, "w");
745
746 if (ntop)
747 {
748 init5(ntop);
749 }
750
751 read_tpxheader(ftp2fn(efTPR, NFILE, fnm), &header, FALSE);
752 snew(xtop, header.natoms);
753 read_tpx(ftp2fn(efTPR, NFILE, fnm), &ir, box, &ntopatoms, xtop, NULL, &mtop);
754 bPDB = opt2bSet("-q", NFILE, fnm);
755 if (bPDB)
756 {
757 snew(xav, ntopatoms);
758 snew(ind_fit, ntopatoms);
759 snew(w_rls, ntopatoms);
760 for (kkk = 0; (kkk < ntopatoms); kkk++)
761 {
762 w_rls[kkk] = 1;
763 ind_fit[kkk] = kkk;
764 }
765
766 snew(atoms, 1);
767 *atoms = gmx_mtop_global_atoms(&mtop);
768
769 if (atoms->pdbinfo == NULL)
770 {
771 snew(atoms->pdbinfo, atoms->nr);
772 }
773 }
774
775 top = gmx_mtop_generate_local_top(&mtop, ir.efep != efepNO);
776
777 g = NULL;
778 pbc_null = NULL;
779 if (ir.ePBC != epbcNONE)
780 {
781 if (ir.bPeriodicMols)
782 {
783 pbc_null = &pbc;
784 }
785 else
786 {
787 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE, FALSE);
788 }
789 }
790
791 if (ftp2bSet(efNDX, NFILE, fnm))
792 {
793 /* TODO: Nothing is written to this file if -c is provided, but it is
794 * still opened... */
795 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
796 xvg = xvgropen(opt2fn("-dr", NFILE, fnm), "Individual Restraints", "Time (ps)",
797 "nm", oenv);
798 snew(vvindex, isize);
799 snew(leg, isize);
800 for (i = 0; (i < isize); i++)
801 {
802 index[i]++;
803 snew(leg[i], 12);
804 sprintf(leg[i], "index %d", index[i]);
805 }
806 xvgr_legend(xvg, isize, (const char**)leg, oenv);
807 }
808 else
809 {
810 isize = 0;
811 }
812
813 ir.dr_tau = 0.0;
814 init_disres(fplog, &mtop, &ir, NULL, &fcd, NULL, FALSE);
815
816 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
817 snew(f, 5*natoms);
818
819 init_dr_res(&dr, fcd.disres.nres);
820 if (opt2bSet("-c", NFILE, fnm))
821 {
822 clust = cluster_index(fplog, opt2fn("-c", NFILE, fnm));
823 snew(dr_clust, clust->clust->nr+1);
824 for (i = 0; (i <= clust->clust->nr); i++)
825 {
826 init_dr_res(&dr_clust[i], fcd.disres.nres);
827 }
828 }
829 else
830 {
831 out = xvgropen(opt2fn("-ds", NFILE, fnm),
832 "Sum of Violations", "Time (ps)", "nm", oenv);
833 aver = xvgropen(opt2fn("-da", NFILE, fnm),
834 "Average Violation", "Time (ps)", "nm", oenv);
835 numv = xvgropen(opt2fn("-dn", NFILE, fnm),
836 "# Violations", "Time (ps)", "#", oenv);
837 maxxv = xvgropen(opt2fn("-dm", NFILE, fnm),
838 "Largest Violation", "Time (ps)", "nm", oenv);
839 }
840
841 mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
842 atoms2md(&mtop, &ir, 0, NULL, mtop.natoms, mdatoms);
843 update_mdatoms(mdatoms, ir.fepvals->init_lambda);
844 init_nrnb(&nrnb);
845 if (ir.ePBC != epbcNONE)
846 {
847 gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
848 }
849
850 j = 0;
851 do
852 {
853 if (ir.ePBC != epbcNONE)
854 {
855 if (ir.bPeriodicMols)
856 {
857 set_pbc(&pbc, ir.ePBC, box);
858 }
859 else
860 {
861 gmx_rmpbc(gpbc, natoms, box, x);
862 }
863 }
864
865 if (clust)
866 {
867 if (j > clust->maxframe)
868 {
869 gmx_fatal(FARGS, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
870 }
871 my_clust = clust->inv_clust[j];
872 range_check(my_clust, 0, clust->clust->nr);
873 check_viol(fplog, &(top->idef.il[F_DISRES]),
874 top->idef.iparams,
875 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
876 }
877 else
878 {
879 check_viol(fplog, &(top->idef.il[F_DISRES]),
880 top->idef.iparams,
881 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
882 }
883 if (bPDB)
884 {
885 reset_x(atoms->nr, ind_fit, atoms->nr, NULL, x, w_rls);
886 do_fit(atoms->nr, w_rls, x, x);
887 if (j == 0)
888 {
889 /* Store the first frame of the trajectory as 'characteristic'
890 * for colouring with violations.
891 */
892 for (kkk = 0; (kkk < atoms->nr); kkk++)
893 {
894 copy_rvec(x[kkk], xav[kkk]);
895 }
896 }
897 }
898 if (!clust)
899 {
900 if (isize > 0)
901 {
902 fprintf(xvg, "%10g", t);
903 for (i = 0; (i < isize); i++)
904 {
905 fprintf(xvg, " %10g", vvindex[i]);
906 }
907 fprintf(xvg, "\n");
908 }
909 fprintf(out, "%10g %10g\n", t, dr.sumv);
910 fprintf(aver, "%10g %10g\n", t, dr.averv);
911 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
912 fprintf(numv, "%10g %10d\n", t, dr.nv);
913 }
914 j++;
915 }
916 while (read_next_x(oenv, status, &t, x, box));
917 close_trj(status);
918 if (ir.ePBC != epbcNONE)
919 {
920 gmx_rmpbc_done(gpbc);
921 }
922
923 if (clust)
924 {
925 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
926 top->idef.iparams, clust->clust, dr_clust,
927 clust->grpname, isize, index);
928 }
929 else
930 {
931 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
932 top->idef.iparams, &dr, isize, index,
933 bPDB ? atoms : NULL);
934 if (bPDB)
935 {
936 write_sto_conf(opt2fn("-q", NFILE, fnm),
937 "Coloured by average violation in Angstrom",
938 atoms, xav, NULL, ir.ePBC, box);
939 }
940 dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres,
941 j, &top->idef, &mtop, max_dr, nlevels, bThird);
942 xvgrclose(out);
943 xvgrclose(aver);
944 xvgrclose(numv);
945 xvgrclose(maxxv);
946 do_view(oenv, opt2fn("-dn", NFILE, fnm), "-nxy");
947 do_view(oenv, opt2fn("-da", NFILE, fnm), "-nxy");
948 do_view(oenv, opt2fn("-ds", NFILE, fnm), "-nxy");
949 do_view(oenv, opt2fn("-dm", NFILE, fnm), "-nxy");
950 }
951 if (isize > 0)
952 {
953 xvgrclose(xvg);
954 if (!clust)
955 {
956 do_view(oenv, opt2fn("-dr", NFILE, fnm), "-nxy");
957 }
958 }
959
960 gmx_ffclose(fplog);
961
962 return 0;
963 }
The ultimate molecular dynamics simulation package