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