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.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/commandline/viewit.h"
48 #include "gromacs/correlationfunctions/autocorr.h"
49 #include "gromacs/fileio/enxio.h"
50 #include "gromacs/fileio/gmxfio.h"
51 #include "gromacs/fileio/tpxio.h"
52 #include "gromacs/fileio/trxio.h"
53 #include "gromacs/fileio/xvgr.h"
54 #include "gromacs/gmxana/gmx_ana.h"
55 #include "gromacs/gmxana/gstat.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/energyoutput.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/topology/ifunc.h"
63 #include "gromacs/topology/mtop_lookup.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/topology/topology.h"
66 #include "gromacs/trajectory/energyframe.h"
67 #include "gromacs/trajectoryanalysis/topologyinformation.h"
68 #include "gromacs/utility/arraysize.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/gmxassert.h"
72 #include "gromacs/utility/pleasecite.h"
73 #include "gromacs/utility/smalloc.h"
74 #include "gromacs/utility/strconvert.h"
76 static real minthird
= -1.0 / 3.0, minsixth
= -1.0 / 6.0;
78 static double mypow(double x
, double y
)
82 return std::pow(x
, y
);
90 static real
blk_value(t_enxblock
* blk
, int sub
, int index
)
92 range_check(index
, 0, blk
->sub
[sub
].nr
);
93 if (blk
->sub
[sub
].type
== xdr_datatype_float
)
95 return blk
->sub
[sub
].fval
[index
];
97 else if (blk
->sub
[sub
].type
== xdr_datatype_double
)
99 return blk
->sub
[sub
].dval
[index
];
103 gmx_incons("Unknown datatype in t_enxblock");
107 static int* select_it(int nre
, char* nm
[], int* nset
)
112 gmx_bool bVerbose
= TRUE
;
114 if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
119 fprintf(stderr
, "Select the terms you want from the following list\n");
120 fprintf(stderr
, "End your selection with 0\n");
124 for (k
= 0; (k
< nre
);)
126 for (j
= 0; (j
< 4) && (k
< nre
); j
++, k
++)
128 fprintf(stderr
, " %3d=%14s", k
+ 1, nm
[k
]);
130 fprintf(stderr
, "\n");
137 if (1 != scanf("%d", &n
))
139 gmx_fatal(FARGS
, "Error reading user input");
141 if ((n
> 0) && (n
<= nre
))
148 for (i
= (*nset
) = 0; (i
< nre
); i
++)
161 static void get_orires_parms(const char* topnm
, t_inputrec
* ir
, int* nor
, int* nex
, int** label
, real
** obs
)
171 read_tpx(topnm
, ir
, box
, &natoms
, nullptr, nullptr, &mtop
);
172 top
= gmx_mtop_t_to_t_topology(&mtop
, FALSE
);
174 ip
= top
.idef
.iparams
;
175 iatom
= top
.idef
.il
[F_ORIRES
].iatoms
;
177 /* Count how many distance restraint there are... */
178 nb
= top
.idef
.il
[F_ORIRES
].nr
;
181 gmx_fatal(FARGS
, "No orientation restraints in topology!\n");
188 for (i
= 0; i
< nb
; i
+= 3)
190 (*label
)[i
/ 3] = ip
[iatom
[i
]].orires
.label
;
191 (*obs
)[i
/ 3] = ip
[iatom
[i
]].orires
.obs
;
192 if (ip
[iatom
[i
]].orires
.ex
>= *nex
)
194 *nex
= ip
[iatom
[i
]].orires
.ex
+ 1;
197 fprintf(stderr
, "Found %d orientation restraints with %d experiments", *nor
, *nex
);
198 done_top_mtop(&top
, &mtop
);
201 static int get_bounds(real
** bounds
, int** index
, int** dr_pair
, int* npairs
, gmx_localtop_t
* top
)
203 t_functype
* functype
;
205 int i
, j
, k
, type
, ftype
, natom
;
212 functype
= top
->idef
.functype
;
213 ip
= top
->idef
.iparams
;
215 /* Count how many distance restraint there are... */
216 nb
= top
->idef
.il
[F_DISRES
].nr
;
219 gmx_fatal(FARGS
, "No distance restraints in topology!\n");
222 /* Allocate memory */
227 /* Fill the bound array */
229 for (i
= 0; (i
< top
->idef
.ntypes
); i
++)
232 if (ftype
== F_DISRES
)
235 label1
= ip
[i
].disres
.label
;
236 b
[nb
] = ip
[i
].disres
.up1
;
243 /* Fill the index array */
245 disres
= &(top
->idef
.il
[F_DISRES
]);
246 iatom
= disres
->iatoms
;
247 for (i
= j
= k
= 0; (i
< disres
->nr
);)
250 ftype
= top
->idef
.functype
[type
];
251 natom
= interaction_function
[ftype
].nratoms
+ 1;
252 if (label1
!= top
->idef
.iparams
[type
].disres
.label
)
255 label1
= top
->idef
.iparams
[type
].disres
.label
;
265 gmx_incons("get_bounds for distance restraints");
275 calc_violations(real rt
[], real rav3
[], int nb
, const int index
[], real bounds
[], real
* viol
, double* st
, double* sa
)
277 const real sixth
= 1.0 / 6.0;
279 double rsum
, rav
, sumaver
, sumt
;
283 for (i
= 0; (i
< nb
); i
++)
287 for (j
= index
[i
]; (j
< index
[i
+ 1]); j
++)
291 viol
[j
] += mypow(rt
[j
], -3.0);
293 rav
+= gmx::square(rav3
[j
]);
294 rsum
+= mypow(rt
[j
], -6);
296 rsum
= std::max(0.0, mypow(rsum
, -sixth
) - bounds
[i
]);
297 rav
= std::max(0.0, mypow(rav
, -sixth
) - bounds
[i
]);
306 static void analyse_disre(const char* voutfn
,
313 const gmx_output_env_t
* oenv
)
316 double sum
, sumt
, sumaver
;
319 /* Subtract bounds from distances, to calculate violations */
320 calc_violations(violaver
, violaver
, nbounds
, pair
, bounds
, nullptr, &sumt
, &sumaver
);
323 fprintf(stdout
, "\nSum of violations averaged over simulation: %g nm\n", sumaver
);
324 fprintf(stdout
, "Largest violation averaged over simulation: %g nm\n\n", sumt
);
326 vout
= xvgropen(voutfn
, "r\\S-3\\N average violations", "DR Index", "nm", oenv
);
329 for (i
= 0; (i
< nbounds
); i
++)
331 /* Do ensemble averaging */
333 for (j
= pair
[i
]; (j
< pair
[i
+ 1]); j
++)
335 sumaver
+= gmx::square(violaver
[j
] / real(nframes
));
337 sumaver
= std::max(0.0, mypow(sumaver
, minsixth
) - bounds
[i
]);
340 sum
= std::max(sum
, sumaver
);
341 fprintf(vout
, "%10d %10.5e\n", index
[i
], sumaver
);
344 for (j
= 0; (j
< dr
.ndr
); j
++)
346 fprintf(vout
, "%10d %10.5e\n", j
, mypow(violaver
[j
] / real(nframes
), minthird
));
351 fprintf(stdout
, "\nSum of violations averaged over simulation: %g nm\n", sumt
);
352 fprintf(stdout
, "Largest violation averaged over simulation: %g nm\n\n", sum
);
354 do_view(oenv
, voutfn
, "-graphtype bar");
357 static void print_time(FILE* fp
, double t
)
359 fprintf(fp
, "%12.6f", t
);
362 int gmx_nmr(int argc
, char* argv
[])
364 const char* desc
[] = {
365 "[THISMODULE] extracts distance or orientation restraint",
366 "data from an energy file. The user is prompted to interactively",
367 "select the desired terms.[PAR]",
369 "When the [TT]-viol[tt] option is set, the time averaged",
370 "violations are plotted and the running time-averaged and",
371 "instantaneous sum of violations are recalculated. Additionally",
372 "running time-averaged and instantaneous distances between",
373 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
375 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
376 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
377 "The first two options plot the orientation, the last three the",
378 "deviations of the orientations from the experimental values.",
379 "The options that end on an 'a' plot the average over time",
380 "as a function of restraint. The options that end on a 't'",
381 "prompt the user for restraint label numbers and plot the data",
382 "as a function of time. Option [TT]-odr[tt] plots the RMS",
383 "deviation as a function of restraint.",
384 "When the run used time or ensemble averaged orientation restraints,",
385 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
386 "not ensemble-averaged orientations and deviations instead of",
387 "the time and ensemble averages.[PAR]",
389 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
390 "tensor for each orientation restraint experiment. With option",
391 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
395 static gmx_bool bPrAll
= FALSE
;
396 static gmx_bool bDp
= FALSE
, bOrinst
= FALSE
, bOvec
= FALSE
;
399 { "-dp", FALSE
, etBOOL
, { &bDp
}, "Print energies in high precision" },
400 { "-skip", FALSE
, etINT
, { &skip
}, "Skip number of frames between data points" },
405 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is "
407 { "-orinst", FALSE
, etBOOL
, { &bOrinst
}, "Analyse instantaneous orientation data" },
408 { "-ovec", FALSE
, etBOOL
, { &bOvec
}, "Also plot the eigenvectors with [TT]-oten[tt]" }
410 const char* drleg
[] = { "Running average", "Instantaneous" };
412 FILE /* *out = NULL,*/ *out_disre
= nullptr, *fp_pairs
= nullptr, *fort
= nullptr,
413 *fodt
= nullptr, *foten
= nullptr;
417 gmx_enxnm_t
* enm
= nullptr;
419 int nre
, teller
, teller_disre
;
420 int nor
= 0, nex
= 0, norfr
= 0, enx_i
= 0;
421 real
*bounds
= nullptr, *violaver
= nullptr, *oobs
= nullptr, *orient
= nullptr, *odrms
= nullptr;
422 int * index
= nullptr, *pair
= nullptr, norsel
= 0, *orsel
= nullptr, *or_label
= nullptr;
423 int nbounds
= 0, npairs
;
424 gmx_bool bDisRe
, bDRAll
, bORA
, bORT
, bODA
, bODR
, bODT
, bORIRE
, bOTEN
;
426 double sumaver
, sumt
;
427 int * set
= nullptr, i
, j
, k
, nset
, sss
;
428 char ** pairleg
, **odtleg
, **otenleg
;
429 char** leg
= nullptr;
430 const char *anm_j
, *anm_k
, *resnm_j
, *resnm_k
;
431 int resnr_j
, resnr_k
;
432 const char* orinst_sub
= "@ subtitle \"instantaneous\"\n";
434 gmx_output_env_t
* oenv
;
435 t_enxblock
* blk_disre
= nullptr;
438 t_filenm fnm
[] = { { efEDR
, "-f", nullptr, ffREAD
},
439 { efEDR
, "-f2", nullptr, ffOPTRD
},
440 { efTPR
, "-s", nullptr, ffOPTRD
},
441 // { efXVG, "-o", "energy", ffWRITE },
442 { efXVG
, "-viol", "violaver", ffOPTWR
},
443 { efXVG
, "-pairs", "pairs", ffOPTWR
},
444 { efXVG
, "-ora", "orienta", ffOPTWR
},
445 { efXVG
, "-ort", "orientt", ffOPTWR
},
446 { efXVG
, "-oda", "orideva", ffOPTWR
},
447 { efXVG
, "-odr", "oridevr", ffOPTWR
},
448 { efXVG
, "-odt", "oridevt", ffOPTWR
},
449 { efXVG
, "-oten", "oriten", ffOPTWR
} };
450 #define NFILE asize(fnm)
454 if (!parse_common_args(&argc
, argv
, PCA_CAN_VIEW
| PCA_CAN_BEGIN
| PCA_CAN_END
, NFILE
, fnm
,
455 npargs
, pa
, asize(desc
), desc
, 0, nullptr, &oenv
))
460 bDRAll
= opt2bSet("-pairs", NFILE
, fnm
);
461 bDisRe
= opt2bSet("-viol", NFILE
, fnm
) || bDRAll
;
462 bORA
= opt2bSet("-ora", NFILE
, fnm
);
463 bORT
= opt2bSet("-ort", NFILE
, fnm
);
464 bODA
= opt2bSet("-oda", NFILE
, fnm
);
465 bODR
= opt2bSet("-odr", NFILE
, fnm
);
466 bODT
= opt2bSet("-odt", NFILE
, fnm
);
467 bORIRE
= bORA
|| bORT
|| bODA
|| bODR
|| bODT
;
468 bOTEN
= opt2bSet("-oten", NFILE
, fnm
);
469 if (!(bDRAll
|| bDisRe
|| bORA
|| bORT
|| bODA
|| bODR
|| bODT
|| bORIRE
|| bOTEN
))
471 printf("No output selected. Run with -h to see options. Terminating program.\n");
476 fp
= open_enx(ftp2fn(efEDR
, NFILE
, fnm
), "r");
477 do_enxnms(fp
, &nre
, &enm
);
478 free_enxnms(nre
, enm
);
480 t_inputrec irInstance
;
481 t_inputrec
* ir
= &irInstance
;
483 gmx::TopologyInformation topInfo
;
488 get_orires_parms(ftp2fn(efTPR
, NFILE
, fnm
), ir
, &nor
, &nex
, &or_label
, &oobs
);
512 fprintf(stderr
, "Select the orientation restraint labels you want (-1 is all)\n");
513 fprintf(stderr
, "End your selection with 0\n");
519 srenew(orsel
, j
+ 1);
520 if (1 != scanf("%d", &(orsel
[j
])))
522 gmx_fatal(FARGS
, "Error reading user input");
524 } while (orsel
[j
] > 0);
527 fprintf(stderr
, "Selecting all %d orientation restraints\n", nor
);
530 for (i
= 0; i
< nor
; i
++)
537 /* Build the selection */
539 for (i
= 0; i
< j
; i
++)
541 for (k
= 0; k
< nor
; k
++)
543 if (or_label
[k
] == orsel
[i
])
552 fprintf(stderr
, "Orientation restraint label %d not found\n", orsel
[i
]);
556 snew(odtleg
, norsel
);
557 for (i
= 0; i
< norsel
; i
++)
559 snew(odtleg
[i
], 256);
560 sprintf(odtleg
[i
], "%d", or_label
[orsel
[i
]]);
564 fort
= xvgropen(opt2fn("-ort", NFILE
, fnm
), "Calculated orientations",
565 "Time (ps)", "", oenv
);
566 if (bOrinst
&& output_env_get_print_xvgr_codes(oenv
))
568 fprintf(fort
, "%s", orinst_sub
);
570 xvgr_legend(fort
, norsel
, odtleg
, oenv
);
574 fodt
= xvgropen(opt2fn("-odt", NFILE
, fnm
), "Orientation restraint deviation",
575 "Time (ps)", "", oenv
);
576 if (bOrinst
&& output_env_get_print_xvgr_codes(oenv
))
578 fprintf(fodt
, "%s", orinst_sub
);
580 xvgr_legend(fodt
, norsel
, odtleg
, oenv
);
582 for (i
= 0; i
< norsel
; i
++)
591 foten
= xvgropen(opt2fn("-oten", NFILE
, fnm
), "Order tensor", "Time (ps)", "", oenv
);
592 snew(otenleg
, bOvec
? nex
* 12 : nex
* 3);
593 for (i
= 0; i
< nex
; i
++)
595 for (j
= 0; j
< 3; j
++)
597 sprintf(buf
, "eig%d", j
+ 1);
598 otenleg
[(bOvec
? 12 : 3) * i
+ j
] = gmx_strdup(buf
);
602 for (j
= 0; j
< 9; j
++)
604 sprintf(buf
, "vec%d%s", j
/ 3 + 1, j
% 3 == 0 ? "x" : (j
% 3 == 1 ? "y" : "z"));
605 otenleg
[12 * i
+ 3 + j
] = gmx_strdup(buf
);
609 xvgr_legend(foten
, bOvec
? nex
* 12 : nex
* 3, otenleg
, oenv
);
610 for (j
= 0; j
< 3; j
++)
620 topInfo
.fillFromInputFile(ftp2fn(efTPR
, NFILE
, fnm
));
621 gmx_mtop_generate_local_top(*topInfo
.mtop(), &top
, ir
->efep
!= efepNO
);
623 nbounds
= get_bounds(&bounds
, &index
, &pair
, &npairs
, &top
);
624 snew(violaver
, npairs
);
625 out_disre
= xvgropen(opt2fn("-o", NFILE
, fnm
), "Sum of Violations", "Time (ps)", "nm", oenv
);
626 xvgr_legend(out_disre
, 2, drleg
, oenv
);
629 fp_pairs
= xvgropen(opt2fn("-pairs", NFILE
, fnm
), "Pair Distances", "Time (ps)",
630 "Distance (nm)", oenv
);
631 if (output_env_get_print_xvgr_codes(oenv
))
633 fprintf(fp_pairs
, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n", ir
->dr_tau
);
638 /* Initiate counters */
643 /* This loop searches for the first frame (when -b option is given),
644 * or when this has been found it reads just one energy frame
648 bCont
= do_enx(fp
, &fr
);
651 timecheck
= check_times(fr
.t
);
653 } while (bCont
&& (timecheck
< 0));
655 if ((timecheck
== 0) && bCont
)
658 * Define distance restraint legends. Can only be done after
659 * the first frame has been read... (Then we know how many there are)
661 blk_disre
= find_block_id_enxframe(&fr
, enxDISRE
, nullptr);
662 if (bDisRe
&& bDRAll
&& !leg
&& blk_disre
)
667 fa
= top
.idef
.il
[F_DISRES
].iatoms
;
668 ip
= top
.idef
.iparams
;
669 if (blk_disre
->nsub
!= 2 || (blk_disre
->sub
[0].nr
!= blk_disre
->sub
[1].nr
))
671 gmx_incons("Number of disre sub-blocks not equal to 2");
674 ndisre
= blk_disre
->sub
[0].nr
;
675 if (ndisre
!= top
.idef
.il
[F_DISRES
].nr
/ 3)
678 "Number of disre pairs in the energy file (%d) does not match the "
679 "number in the run input file (%d)\n",
680 ndisre
, top
.idef
.il
[F_DISRES
].nr
/ 3);
682 snew(pairleg
, ndisre
);
684 for (i
= 0; i
< ndisre
; i
++)
686 snew(pairleg
[i
], 30);
689 mtopGetAtomAndResidueName(topInfo
.mtop(), j
, &molb
, &anm_j
, &resnr_j
, &resnm_j
, nullptr);
690 mtopGetAtomAndResidueName(topInfo
.mtop(), k
, &molb
, &anm_k
, &resnr_k
, &resnm_k
, nullptr);
691 sprintf(pairleg
[i
], "%d %s %d %s (%d)", resnr_j
, anm_j
, resnr_k
, anm_k
,
692 ip
[fa
[3 * i
]].disres
.label
);
694 set
= select_it(ndisre
, pairleg
, &nset
);
696 for (i
= 0; (i
< nset
); i
++)
698 snew(leg
[2 * i
], 32);
699 sprintf(leg
[2 * i
], "a %s", pairleg
[set
[i
]]);
700 snew(leg
[2 * i
+ 1], 32);
701 sprintf(leg
[2 * i
+ 1], "i %s", pairleg
[set
[i
]]);
703 xvgr_legend(fp_pairs
, 2 * nset
, leg
, oenv
);
707 * Printing time, only when we do not want to skip frames
709 if (!skip
|| teller
% skip
== 0)
713 /*******************************************
714 * D I S T A N C E R E S T R A I N T S
715 *******************************************/
718 GMX_RELEASE_ASSERT(blk_disre
!= nullptr,
719 "Trying to dereference NULL blk_disre pointer");
721 float* disre_rt
= blk_disre
->sub
[0].fval
;
722 float* disre_rm3tav
= blk_disre
->sub
[1].fval
;
724 double* disre_rt
= blk_disre
->sub
[0].dval
;
725 double* disre_rm3tav
= blk_disre
->sub
[1].dval
;
728 print_time(out_disre
, fr
.t
);
729 if (violaver
== nullptr)
731 snew(violaver
, ndisre
);
734 /* Subtract bounds from distances, to calculate violations */
735 calc_violations(disre_rt
, disre_rm3tav
, nbounds
, pair
, bounds
, violaver
,
738 fprintf(out_disre
, " %8.4f %8.4f\n", sumaver
, sumt
);
741 print_time(fp_pairs
, fr
.t
);
742 for (i
= 0; (i
< nset
); i
++)
745 fprintf(fp_pairs
, " %8.4f", mypow(disre_rm3tav
[sss
], minthird
));
746 fprintf(fp_pairs
, " %8.4f", disre_rt
[sss
]);
748 fprintf(fp_pairs
, "\n");
754 /*******************************************
756 *******************************************/
759 t_enxblock
* blk
= find_block_id_enxframe(&fr
, enx_i
, nullptr);
764 gmx_fatal(FARGS
, "Orientational restraints read in incorrectly.");
767 if (blk
->sub
[0].nr
!= nor
)
770 "Number of orientation restraints in energy file (%d) does "
771 "not match with the topology (%d)",
772 blk
->sub
[0].nr
, nor
);
776 for (i
= 0; i
< nor
; i
++)
778 orient
[i
] += blk_value(blk
, 0, i
);
783 for (i
= 0; i
< nor
; i
++)
785 real v
= blk_value(blk
, 0, i
);
786 odrms
[i
] += gmx::square(v
- oobs
[i
]);
791 fprintf(fort
, " %10f", fr
.t
);
792 for (i
= 0; i
< norsel
; i
++)
794 fprintf(fort
, " %g", blk_value(blk
, 0, orsel
[i
]));
800 fprintf(fodt
, " %10f", fr
.t
);
801 for (i
= 0; i
< norsel
; i
++)
803 fprintf(fodt
, " %g", blk_value(blk
, 0, orsel
[i
]) - oobs
[orsel
[i
]]);
809 blk
= find_block_id_enxframe(&fr
, enxORT
, nullptr);
814 gmx_fatal(FARGS
, "Orientational restraints read in incorrectly");
817 if (blk
->sub
[0].nr
!= nex
* 12)
820 "Number of orientation experiments in energy file (%d) does "
821 "not match with the topology (%d)",
822 blk
->sub
[0].nr
/ 12, nex
);
824 fprintf(foten
, " %10f", fr
.t
);
825 for (i
= 0; i
< nex
; i
++)
827 for (j
= 0; j
< (bOvec
? 12 : 3); j
++)
829 fprintf(foten
, " %g", blk_value(blk
, 0, i
* 12 + j
));
832 fprintf(foten
, "\n");
838 } while (bCont
&& (timecheck
== 0));
841 fprintf(stderr
, "\n");
845 xvgrclose(out_disre
);
863 FILE* out
= xvgropen(opt2fn("-ora", NFILE
, fnm
), "Average calculated orientations",
864 "Restraint label", "", oenv
);
865 if (bOrinst
&& output_env_get_print_xvgr_codes(oenv
))
867 fprintf(out
, "%s", orinst_sub
);
869 for (i
= 0; i
< nor
; i
++)
871 fprintf(out
, "%5d %g\n", or_label
[i
], orient
[i
] / real(norfr
));
877 FILE* out
= xvgropen(opt2fn("-oda", NFILE
, fnm
), "Average restraint deviation",
878 "Restraint label", "", oenv
);
879 if (bOrinst
&& output_env_get_print_xvgr_codes(oenv
))
881 fprintf(out
, "%s", orinst_sub
);
883 for (i
= 0; i
< nor
; i
++)
885 fprintf(out
, "%5d %g\n", or_label
[i
], orient
[i
] / real(norfr
) - oobs
[i
]);
891 FILE* out
= xvgropen(opt2fn("-odr", NFILE
, fnm
), "RMS orientation restraint deviations",
892 "Restraint label", "", oenv
);
893 if (bOrinst
&& output_env_get_print_xvgr_codes(oenv
))
895 fprintf(out
, "%s", orinst_sub
);
897 for (i
= 0; i
< nor
; i
++)
899 fprintf(out
, "%5d %g\n", or_label
[i
], std::sqrt(odrms
[i
] / real(norfr
)));
903 // Clean up orires variables.
916 analyse_disre(opt2fn("-viol", NFILE
, fnm
), teller_disre
, violaver
, bounds
, index
, pair
,
920 const char* nxy
= "-nxy";
922 do_view(oenv
, opt2fn_null("-ora", NFILE
, fnm
), nxy
);
923 do_view(oenv
, opt2fn_null("-ort", NFILE
, fnm
), nxy
);
924 do_view(oenv
, opt2fn_null("-oda", NFILE
, fnm
), nxy
);
925 do_view(oenv
, opt2fn_null("-odr", NFILE
, fnm
), nxy
);
926 do_view(oenv
, opt2fn_null("-odt", NFILE
, fnm
), nxy
);
927 do_view(oenv
, opt2fn_null("-oten", NFILE
, fnm
), nxy
);
929 output_env_done(oenv
);