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,2018, 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.
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.
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.
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.
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.
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.
47 #include "gromacs/awh/awh.h"
48 #include "gromacs/fileio/enxio.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/listed-forces/disre.h"
53 #include "gromacs/listed-forces/orires.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdlib/constr.h"
58 #include "gromacs/mdlib/mdebin_bar.h"
59 #include "gromacs/mdlib/mdrun.h"
60 #include "gromacs/mdtypes/energyhistory.h"
61 #include "gromacs/mdtypes/fcdata.h"
62 #include "gromacs/mdtypes/group.h"
63 #include "gromacs/mdtypes/inputrec.h"
64 #include "gromacs/mdtypes/md_enums.h"
65 #include "gromacs/mdtypes/state.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/pulling/pull.h"
68 #include "gromacs/topology/mtop_util.h"
69 #include "gromacs/trajectory/energyframe.h"
70 #include "gromacs/utility/arraysize.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/smalloc.h"
73 #include "gromacs/utility/stringutil.h"
75 static const char *conrmsd_nm
[] = { "Constr. rmsd", "Constr.2 rmsd" };
77 static const char *boxs_nm
[] = { "Box-X", "Box-Y", "Box-Z" };
79 static const char *tricl_boxs_nm
[] = {
80 "Box-XX", "Box-YY", "Box-ZZ",
81 "Box-YX", "Box-ZX", "Box-ZY"
84 static const char *vol_nm
[] = { "Volume" };
86 static const char *dens_nm
[] = {"Density" };
88 static const char *pv_nm
[] = {"pV" };
90 static const char *enthalpy_nm
[] = {"Enthalpy" };
92 static const char *boxvel_nm
[] = {
93 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
94 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
97 #define NBOXS asize(boxs_nm)
98 #define NTRICLBOXS asize(tricl_boxs_nm)
100 const char *egrp_nm
[egNR
+1] = {
101 "Coul-SR", "LJ-SR", "Buck-SR",
102 "Coul-14", "LJ-14", nullptr
105 t_mdebin
*init_mdebin(ener_file_t fp_ene
,
106 const gmx_mtop_t
*mtop
,
107 const t_inputrec
*ir
,
111 const char *ener_nm
[F_NRE
];
112 static const char *vir_nm
[] = {
113 "Vir-XX", "Vir-XY", "Vir-XZ",
114 "Vir-YX", "Vir-YY", "Vir-YZ",
115 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
117 static const char *sv_nm
[] = {
118 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
119 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
120 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
122 static const char *fv_nm
[] = {
123 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
124 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
125 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
127 static const char *pres_nm
[] = {
128 "Pres-XX", "Pres-XY", "Pres-XZ",
129 "Pres-YX", "Pres-YY", "Pres-YZ",
130 "Pres-ZX", "Pres-ZY", "Pres-ZZ"
132 static const char *surft_nm
[] = {
135 static const char *mu_nm
[] = {
136 "Mu-X", "Mu-Y", "Mu-Z"
138 static const char *vcos_nm
[] = {
141 static const char *visc_nm
[] = {
144 static const char *baro_nm
[] = {
148 const gmx_groups_t
*groups
;
153 int i
, j
, ni
, nj
, n
, k
, kk
, ncon
, nset
;
158 if (EI_DYNAMICS(ir
->eI
))
160 md
->delta_t
= ir
->delta_t
;
167 groups
= &mtop
->groups
;
169 bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
170 b14
= (gmx_mtop_ftype_count(mtop
, F_LJ14
) > 0 ||
171 gmx_mtop_ftype_count(mtop
, F_LJC14_Q
) > 0);
173 ncon
= gmx_mtop_ftype_count(mtop
, F_CONSTR
);
174 nset
= gmx_mtop_ftype_count(mtop
, F_SETTLE
);
175 md
->bConstr
= (ncon
> 0 || nset
> 0) && !isRerun
;
176 md
->bConstrVir
= FALSE
;
179 if (ncon
> 0 && ir
->eConstrAlg
== econtLINCS
)
183 md
->bConstrVir
= (getenv("GMX_CONSTRAINTVIR") != nullptr);
190 /* Energy monitoring */
191 for (i
= 0; i
< egNR
; i
++)
193 md
->bEInd
[i
] = FALSE
;
196 for (i
= 0; i
< F_NRE
; i
++)
198 md
->bEner
[i
] = FALSE
;
200 (i
== F_EKIN
|| i
== F_ETOT
|| i
== F_ECONSERVED
||
201 i
== F_TEMP
|| i
== F_PDISPCORR
|| i
== F_PRES
))
207 md
->bEner
[i
] = !bBHAM
;
209 else if (i
== F_BHAM
)
211 md
->bEner
[i
] = bBHAM
;
215 md
->bEner
[i
] = ir
->bQMMM
;
217 else if (i
== F_RF_EXCL
)
219 md
->bEner
[i
] = (EEL_RF(ir
->coulombtype
) && ir
->cutoff_scheme
== ecutsGROUP
);
221 else if (i
== F_COUL_RECIP
)
223 md
->bEner
[i
] = EEL_FULL(ir
->coulombtype
);
225 else if (i
== F_LJ_RECIP
)
227 md
->bEner
[i
] = EVDW_PME(ir
->vdwtype
);
229 else if (i
== F_LJ14
)
233 else if (i
== F_COUL14
)
237 else if (i
== F_LJC14_Q
|| i
== F_LJC_PAIRS_NB
)
239 md
->bEner
[i
] = FALSE
;
241 else if ((i
== F_DVDL_COUL
&& ir
->fepvals
->separate_dvdl
[efptCOUL
]) ||
242 (i
== F_DVDL_VDW
&& ir
->fepvals
->separate_dvdl
[efptVDW
]) ||
243 (i
== F_DVDL_BONDED
&& ir
->fepvals
->separate_dvdl
[efptBONDED
]) ||
244 (i
== F_DVDL_RESTRAINT
&& ir
->fepvals
->separate_dvdl
[efptRESTRAINT
]) ||
245 (i
== F_DKDL
&& ir
->fepvals
->separate_dvdl
[efptMASS
]) ||
246 (i
== F_DVDL
&& ir
->fepvals
->separate_dvdl
[efptFEP
]))
248 md
->bEner
[i
] = (ir
->efep
!= efepNO
);
250 else if ((interaction_function
[i
].flags
& IF_VSITE
) ||
251 (i
== F_CONSTR
) || (i
== F_CONSTRNC
) || (i
== F_SETTLE
))
253 md
->bEner
[i
] = FALSE
;
255 else if ((i
== F_COUL_SR
) || (i
== F_EPOT
) || (i
== F_PRES
) || (i
== F_EQM
))
259 else if ((i
== F_ETOT
) || (i
== F_EKIN
) || (i
== F_TEMP
))
261 md
->bEner
[i
] = EI_DYNAMICS(ir
->eI
);
263 else if (i
== F_DISPCORR
|| i
== F_PDISPCORR
)
265 md
->bEner
[i
] = (ir
->eDispCorr
!= edispcNO
);
267 else if (i
== F_DISRESVIOL
)
269 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_DISRES
) > 0);
271 else if (i
== F_ORIRESDEV
)
273 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_ORIRES
) > 0);
275 else if (i
== F_CONNBONDS
)
277 md
->bEner
[i
] = FALSE
;
279 else if (i
== F_COM_PULL
)
281 md
->bEner
[i
] = ((ir
->bPull
&& pull_have_potential(ir
->pull_work
)) ||
284 else if (i
== F_ECONSERVED
)
286 md
->bEner
[i
] = (integratorHasConservedEnergyQuantity(ir
));
290 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, i
) > 0);
295 for (i
= 0; i
< F_NRE
; i
++)
299 ener_nm
[md
->f_nre
] = interaction_function
[i
].longname
;
304 md
->epc
= isRerun
? epcNO
: ir
->epc
;
305 md
->bDiagPres
= !TRICLINIC(ir
->ref_p
) && !isRerun
;
306 md
->ref_p
= (ir
->ref_p
[XX
][XX
]+ir
->ref_p
[YY
][YY
]+ir
->ref_p
[ZZ
][ZZ
])/DIM
;
307 md
->bTricl
= TRICLINIC(ir
->compress
) || TRICLINIC(ir
->deform
);
308 md
->bDynBox
= inputrecDynamicBox(ir
);
309 md
->etc
= isRerun
? etcNO
: ir
->etc
;
310 md
->bNHC_trotter
= inputrecNvtTrotter(ir
) && !isRerun
;
311 md
->bPrintNHChains
= ir
->bPrintNHChains
&& !isRerun
;
312 md
->bMTTK
= (inputrecNptTrotter(ir
) || inputrecNphTrotter(ir
)) && !isRerun
;
313 md
->bMu
= inputrecNeedMutot(ir
);
314 md
->bPres
= !isRerun
;
316 md
->ebin
= mk_ebin();
317 /* Pass NULL for unit to let get_ebin_space determine the units
318 * for interaction_function[i].longname
320 md
->ie
= get_ebin_space(md
->ebin
, md
->f_nre
, ener_nm
, nullptr);
323 /* This should be called directly after the call for md->ie,
324 * such that md->iconrmsd follows directly in the list.
326 md
->iconrmsd
= get_ebin_space(md
->ebin
, md
->nCrmsd
, conrmsd_nm
, "");
330 md
->ib
= get_ebin_space(md
->ebin
,
331 md
->bTricl
? NTRICLBOXS
: NBOXS
,
332 md
->bTricl
? tricl_boxs_nm
: boxs_nm
,
334 md
->ivol
= get_ebin_space(md
->ebin
, 1, vol_nm
, unit_volume
);
335 md
->idens
= get_ebin_space(md
->ebin
, 1, dens_nm
, unit_density_SI
);
338 md
->ipv
= get_ebin_space(md
->ebin
, 1, pv_nm
, unit_energy
);
339 md
->ienthalpy
= get_ebin_space(md
->ebin
, 1, enthalpy_nm
, unit_energy
);
344 md
->isvir
= get_ebin_space(md
->ebin
, asize(sv_nm
), sv_nm
, unit_energy
);
345 md
->ifvir
= get_ebin_space(md
->ebin
, asize(fv_nm
), fv_nm
, unit_energy
);
349 md
->ivir
= get_ebin_space(md
->ebin
, asize(vir_nm
), vir_nm
, unit_energy
);
350 md
->ipres
= get_ebin_space(md
->ebin
, asize(pres_nm
), pres_nm
, unit_pres_bar
);
351 md
->isurft
= get_ebin_space(md
->ebin
, asize(surft_nm
), surft_nm
,
354 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
356 md
->ipc
= get_ebin_space(md
->ebin
, md
->bTricl
? 6 : 3,
357 boxvel_nm
, unit_vel
);
361 md
->imu
= get_ebin_space(md
->ebin
, asize(mu_nm
), mu_nm
, unit_dipole_D
);
363 if (ir
->cos_accel
!= 0)
365 md
->ivcos
= get_ebin_space(md
->ebin
, asize(vcos_nm
), vcos_nm
, unit_vel
);
366 md
->ivisc
= get_ebin_space(md
->ebin
, asize(visc_nm
), visc_nm
,
370 /* Energy monitoring */
371 for (i
= 0; i
< egNR
; i
++)
373 md
->bEInd
[i
] = FALSE
;
375 md
->bEInd
[egCOULSR
] = TRUE
;
376 md
->bEInd
[egLJSR
] = TRUE
;
380 md
->bEInd
[egLJSR
] = FALSE
;
381 md
->bEInd
[egBHAMSR
] = TRUE
;
385 md
->bEInd
[egLJ14
] = TRUE
;
386 md
->bEInd
[egCOUL14
] = TRUE
;
389 for (i
= 0; (i
< egNR
); i
++)
397 n
= groups
->grps
[egcENER
].nr
;
399 md
->nE
= (n
*(n
+1))/2;
401 snew(md
->igrp
, md
->nE
);
406 for (k
= 0; (k
< md
->nEc
); k
++)
408 snew(gnm
[k
], STRLEN
);
410 for (i
= 0; (i
< groups
->grps
[egcENER
].nr
); i
++)
412 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
413 for (j
= i
; (j
< groups
->grps
[egcENER
].nr
); j
++)
415 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
416 for (k
= kk
= 0; (k
< egNR
); k
++)
420 sprintf(gnm
[kk
], "%s:%s-%s", egrp_nm
[k
],
421 *(groups
->grpname
[ni
]), *(groups
->grpname
[nj
]));
425 md
->igrp
[n
] = get_ebin_space(md
->ebin
, md
->nEc
,
430 for (k
= 0; (k
< md
->nEc
); k
++)
438 gmx_incons("Number of energy terms wrong");
442 md
->nTC
= isRerun
? 0 : groups
->grps
[egcTC
].nr
;
443 md
->nNHC
= ir
->opts
.nhchainlength
; /* shorthand for number of NH chains */
446 md
->nTCP
= 1; /* assume only one possible coupling system for barostat
453 if (md
->etc
== etcNOSEHOOVER
)
455 if (md
->bNHC_trotter
)
457 md
->mde_n
= 2*md
->nNHC
*md
->nTC
;
461 md
->mde_n
= 2*md
->nTC
;
463 if (md
->epc
== epcMTTK
)
465 md
->mdeb_n
= 2*md
->nNHC
*md
->nTCP
;
474 snew(md
->tmp_r
, md
->mde_n
);
475 snew(md
->tmp_v
, md
->mde_n
);
477 snew(grpnms
, md
->mde_n
);
479 for (i
= 0; (i
< md
->nTC
); i
++)
481 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
482 sprintf(buf
, "T-%s", *(groups
->grpname
[ni
]));
483 grpnms
[i
] = gmx_strdup(buf
);
485 md
->itemp
= get_ebin_space(md
->ebin
, md
->nTC
, grpnms
,
488 if (md
->etc
== etcNOSEHOOVER
)
490 if (md
->bPrintNHChains
)
492 if (md
->bNHC_trotter
)
494 for (i
= 0; (i
< md
->nTC
); i
++)
496 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
497 bufi
= *(groups
->grpname
[ni
]);
498 for (j
= 0; (j
< md
->nNHC
); j
++)
500 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
501 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
502 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
503 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
506 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
507 grpnms
, unit_invtime
);
510 for (i
= 0; (i
< md
->nTCP
); i
++)
512 bufi
= baro_nm
[0]; /* All barostat DOF's together for now. */
513 for (j
= 0; (j
< md
->nNHC
); j
++)
515 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
516 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
517 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
518 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
521 md
->itcb
= get_ebin_space(md
->ebin
, md
->mdeb_n
,
522 grpnms
, unit_invtime
);
527 for (i
= 0; (i
< md
->nTC
); i
++)
529 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
530 bufi
= *(groups
->grpname
[ni
]);
531 sprintf(buf
, "Xi-%s", bufi
);
532 grpnms
[2*i
] = gmx_strdup(buf
);
533 sprintf(buf
, "vXi-%s", bufi
);
534 grpnms
[2*i
+1] = gmx_strdup(buf
);
536 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
537 grpnms
, unit_invtime
);
541 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
542 md
->etc
== etcVRESCALE
)
544 for (i
= 0; (i
< md
->nTC
); i
++)
546 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
547 sprintf(buf
, "Lamb-%s", *(groups
->grpname
[ni
]));
548 grpnms
[i
] = gmx_strdup(buf
);
550 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
, grpnms
, "");
553 for (i
= 0; i
< md
->mde_n
; i
++)
559 md
->nU
= groups
->grps
[egcACC
].nr
;
562 snew(grpnms
, 3*md
->nU
);
563 for (i
= 0; (i
< md
->nU
); i
++)
565 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
566 sprintf(buf
, "Ux-%s", *(groups
->grpname
[ni
]));
567 grpnms
[3*i
+XX
] = gmx_strdup(buf
);
568 sprintf(buf
, "Uy-%s", *(groups
->grpname
[ni
]));
569 grpnms
[3*i
+YY
] = gmx_strdup(buf
);
570 sprintf(buf
, "Uz-%s", *(groups
->grpname
[ni
]));
571 grpnms
[3*i
+ZZ
] = gmx_strdup(buf
);
573 md
->iu
= get_ebin_space(md
->ebin
, 3*md
->nU
, grpnms
, unit_vel
);
579 do_enxnms(fp_ene
, &md
->ebin
->nener
, &md
->ebin
->enm
);
582 md
->print_grpnms
= nullptr;
584 /* check whether we're going to write dh histograms */
586 if (ir
->fepvals
->separate_dhdl_file
== esepdhdlfileNO
)
588 /* Currently dh histograms are only written with dynamics */
589 if (EI_DYNAMICS(ir
->eI
))
593 mde_delta_h_coll_init(md
->dhc
, ir
);
595 md
->fp_dhdl
= nullptr;
596 snew(md
->dE
, ir
->fepvals
->n_lambda
);
600 md
->fp_dhdl
= fp_dhdl
;
601 snew(md
->dE
, ir
->fepvals
->n_lambda
);
606 snew(md
->temperatures
, ir
->fepvals
->n_lambda
);
607 for (i
= 0; i
< ir
->fepvals
->n_lambda
; i
++)
609 md
->temperatures
[i
] = ir
->simtempvals
->temperatures
[i
];
615 void done_mdebin(t_mdebin
*mdebin
)
618 sfree(mdebin
->tmp_r
);
619 sfree(mdebin
->tmp_v
);
620 done_ebin(mdebin
->ebin
);
621 done_mde_delta_h_coll(mdebin
->dhc
);
623 sfree(mdebin
->temperatures
);
628 /* print a lambda vector to a string
629 fep = the inputrec's FEP input data
630 i = the index of the lambda vector
631 get_native_lambda = whether to print the native lambda
632 get_names = whether to print the names rather than the values
633 str = the pre-allocated string buffer to print to. */
634 static void print_lambda_vector(t_lambda
*fep
, int i
,
635 gmx_bool get_native_lambda
, gmx_bool get_names
,
641 for (j
= 0; j
< efptNR
; j
++)
643 if (fep
->separate_dvdl
[j
])
648 str
[0] = 0; /* reset the string */
651 str
+= sprintf(str
, "("); /* set the opening parenthesis*/
653 for (j
= 0; j
< efptNR
; j
++)
655 if (fep
->separate_dvdl
[j
])
659 if (get_native_lambda
&& fep
->init_lambda
>= 0)
661 str
+= sprintf(str
, "%.4f", fep
->init_lambda
);
665 str
+= sprintf(str
, "%.4f", fep
->all_lambda
[j
][i
]);
670 str
+= sprintf(str
, "%s", efpt_singular_names
[j
]);
672 /* print comma for the next item */
675 str
+= sprintf(str
, ", ");
682 /* and add the closing parenthesis */
688 extern FILE *open_dhdl(const char *filename
, const t_inputrec
*ir
,
689 const gmx_output_env_t
*oenv
)
692 const char *dhdl
= "dH/d\\lambda", *deltag
= "\\DeltaH", *lambda
= "\\lambda",
693 *lambdastate
= "\\lambda state";
694 int i
, nsets
, nsets_de
, nsetsbegin
;
695 int n_lambda_terms
= 0;
696 t_lambda
*fep
= ir
->fepvals
; /* for simplicity */
697 t_expanded
*expand
= ir
->expandedvals
;
698 char lambda_vec_str
[STRLEN
], lambda_name_str
[STRLEN
];
703 gmx_bool write_pV
= FALSE
;
705 /* count the number of different lambda terms */
706 for (i
= 0; i
< efptNR
; i
++)
708 if (fep
->separate_dvdl
[i
])
714 std::string title
, label_x
, label_y
;
715 if (fep
->n_lambda
== 0)
717 title
= gmx::formatString("%s", dhdl
);
718 label_x
= gmx::formatString("Time (ps)");
719 label_y
= gmx::formatString("%s (%s %s)",
720 dhdl
, unit_energy
, "[\\lambda]\\S-1\\N");
724 title
= gmx::formatString("%s and %s", dhdl
, deltag
);
725 label_x
= gmx::formatString("Time (ps)");
726 label_y
= gmx::formatString("%s and %s (%s %s)",
727 dhdl
, deltag
, unit_energy
, "[\\8l\\4]\\S-1\\N");
729 fp
= gmx_fio_fopen(filename
, "w+");
730 xvgr_header(fp
, title
.c_str(), label_x
, label_y
, exvggtXNY
, oenv
);
735 buf
= gmx::formatString("T = %g (K) ", ir
->opts
.ref_t
[0]);
737 if ((ir
->efep
!= efepSLOWGROWTH
) && (ir
->efep
!= efepEXPANDED
))
739 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
741 /* compatibility output */
742 buf
+= gmx::formatString("%s = %.4f", lambda
, fep
->init_lambda
);
746 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, FALSE
,
748 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, TRUE
,
750 buf
+= gmx::formatString("%s %d: %s = %s",
751 lambdastate
, fep
->init_fep_state
,
752 lambda_name_str
, lambda_vec_str
);
755 xvgr_subtitle(fp
, buf
.c_str(), oenv
);
759 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
761 nsets_dhdl
= n_lambda_terms
;
763 /* count the number of delta_g states */
764 nsets_de
= fep
->lambda_stop_n
- fep
->lambda_start_n
;
766 nsets
= nsets_dhdl
+ nsets_de
; /* dhdl + fep differences */
768 if (fep
->n_lambda
> 0 && (expand
->elmcmove
> elmcmoveNO
))
770 nsets
+= 1; /*add fep state for expanded ensemble */
773 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
775 nsets
+= 1; /* add energy to the dhdl as well */
779 if ((ir
->epc
!= epcNO
) && (fep
->n_lambda
> 0) && (fep
->init_lambda
< 0))
781 nsetsextend
+= 1; /* for PV term, other terms possible if required for
782 the reduced potential (only needed with foreign
783 lambda, and only output when init_lambda is not
784 set in order to maintain compatibility of the
788 std::vector
<std::string
> setname(nsetsextend
);
790 if (expand
->elmcmove
> elmcmoveNO
)
792 /* state for the fep_vals, if we have alchemical sampling */
793 setname
[s
++] = "Thermodynamic state";
796 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
799 switch (fep
->edHdLPrintEnergy
)
801 case edHdLPrintEnergyPOTENTIAL
:
802 energy
= gmx::formatString("%s (%s)", "Potential Energy", unit_energy
);
804 case edHdLPrintEnergyTOTAL
:
805 case edHdLPrintEnergyYES
:
807 energy
= gmx::formatString("%s (%s)", "Total Energy", unit_energy
);
809 setname
[s
++] = energy
;
812 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
814 for (i
= 0; i
< efptNR
; i
++)
816 if (fep
->separate_dvdl
[i
])
818 std::string derivative
;
819 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
821 /* compatibility output */
822 derivative
= gmx::formatString("%s %s %.4f", dhdl
, lambda
, fep
->init_lambda
);
826 double lam
= fep
->init_lambda
;
827 if (fep
->init_lambda
< 0)
829 lam
= fep
->all_lambda
[i
][fep
->init_fep_state
];
831 derivative
= gmx::formatString("%s %s = %.4f", dhdl
, efpt_singular_names
[i
],
834 setname
[s
++] = derivative
;
839 if (fep
->n_lambda
> 0)
841 /* g_bar has to determine the lambda values used in this simulation
842 * from this xvg legend.
845 if (expand
->elmcmove
> elmcmoveNO
)
847 nsetsbegin
= 1; /* for including the expanded ensemble */
854 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
858 nsetsbegin
+= nsets_dhdl
;
860 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
862 print_lambda_vector(fep
, i
, FALSE
, FALSE
, lambda_vec_str
);
864 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
866 /* for compatible dhdl.xvg files */
867 buf
= gmx::formatString("%s %s %s", deltag
, lambda
, lambda_vec_str
);
871 buf
= gmx::formatString("%s %s to %s", deltag
, lambda
, lambda_vec_str
);
876 /* print the temperature for this state if doing simulated annealing */
877 buf
+= gmx::formatString("T = %g (%s)",
878 ir
->simtempvals
->temperatures
[s
-(nsetsbegin
)],
885 setname
[s
++] = gmx::formatString("pV (%s)", unit_energy
);
888 xvgrLegend(fp
, setname
, oenv
);
894 static void copy_energy(t_mdebin
*md
, const real e
[], real ecpy
[])
898 for (i
= j
= 0; (i
< F_NRE
); i
++)
907 gmx_incons("Number of energy terms wrong");
911 void upd_mdebin(t_mdebin
*md
,
916 gmx_enerdata_t
*enerd
,
925 gmx_ekindata_t
*ekind
,
927 const gmx::Constraints
*constr
)
929 int i
, j
, k
, kk
, n
, gid
;
930 real crmsd
[2], tmp6
[6];
931 real bs
[NTRICLBOXS
], vol
, dens
, pv
, enthalpy
;
934 double store_dhdl
[efptNR
];
935 real store_energy
= 0;
938 /* Do NOT use the box in the state variable, but the separate box provided
939 * as an argument. This is because we sometimes need to write the box from
940 * the last timestep to match the trajectory frames.
942 copy_energy(md
, enerd
->term
, ecopy
);
943 add_ebin(md
->ebin
, md
->ie
, md
->f_nre
, ecopy
, bSum
);
946 crmsd
[0] = constr
->rmsd();
947 add_ebin(md
->ebin
, md
->iconrmsd
, md
->nCrmsd
, crmsd
, FALSE
);
969 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
970 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
971 add_ebin(md
->ebin
, md
->ib
, nboxs
, bs
, bSum
);
972 add_ebin(md
->ebin
, md
->ivol
, 1, &vol
, bSum
);
973 add_ebin(md
->ebin
, md
->idens
, 1, &dens
, bSum
);
977 /* This is pV (in kJ/mol). The pressure is the reference pressure,
978 not the instantaneous pressure */
979 pv
= vol
*md
->ref_p
/PRESFAC
;
981 add_ebin(md
->ebin
, md
->ipv
, 1, &pv
, bSum
);
982 enthalpy
= pv
+ enerd
->term
[F_ETOT
];
983 add_ebin(md
->ebin
, md
->ienthalpy
, 1, &enthalpy
, bSum
);
988 add_ebin(md
->ebin
, md
->isvir
, 9, svir
[0], bSum
);
989 add_ebin(md
->ebin
, md
->ifvir
, 9, fvir
[0], bSum
);
993 add_ebin(md
->ebin
, md
->ivir
, 9, vir
[0], bSum
);
994 add_ebin(md
->ebin
, md
->ipres
, 9, pres
[0], bSum
);
995 tmp
= (pres
[ZZ
][ZZ
]-(pres
[XX
][XX
]+pres
[YY
][YY
])*0.5)*box
[ZZ
][ZZ
];
996 add_ebin(md
->ebin
, md
->isurft
, 1, &tmp
, bSum
);
998 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
1000 tmp6
[0] = state
->boxv
[XX
][XX
];
1001 tmp6
[1] = state
->boxv
[YY
][YY
];
1002 tmp6
[2] = state
->boxv
[ZZ
][ZZ
];
1003 tmp6
[3] = state
->boxv
[YY
][XX
];
1004 tmp6
[4] = state
->boxv
[ZZ
][XX
];
1005 tmp6
[5] = state
->boxv
[ZZ
][YY
];
1006 add_ebin(md
->ebin
, md
->ipc
, md
->bTricl
? 6 : 3, tmp6
, bSum
);
1010 add_ebin(md
->ebin
, md
->imu
, 3, mu_tot
, bSum
);
1012 if (ekind
&& ekind
->cosacc
.cos_accel
!= 0)
1014 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
1015 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
1016 add_ebin(md
->ebin
, md
->ivcos
, 1, &(ekind
->cosacc
.vcos
), bSum
);
1017 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
1018 tmp
= 1/(ekind
->cosacc
.cos_accel
/(ekind
->cosacc
.vcos
*PICO
)
1019 *dens
*gmx::square(box
[ZZ
][ZZ
]*NANO
/(2*M_PI
)));
1020 add_ebin(md
->ebin
, md
->ivisc
, 1, &tmp
, bSum
);
1025 for (i
= 0; (i
< md
->nEg
); i
++)
1027 for (j
= i
; (j
< md
->nEg
); j
++)
1029 gid
= GID(i
, j
, md
->nEg
);
1030 for (k
= kk
= 0; (k
< egNR
); k
++)
1034 eee
[kk
++] = enerd
->grpp
.ener
[k
][gid
];
1037 add_ebin(md
->ebin
, md
->igrp
[n
], md
->nEc
, eee
, bSum
);
1045 for (i
= 0; (i
< md
->nTC
); i
++)
1047 md
->tmp_r
[i
] = ekind
->tcstat
[i
].T
;
1049 add_ebin(md
->ebin
, md
->itemp
, md
->nTC
, md
->tmp_r
, bSum
);
1051 if (md
->etc
== etcNOSEHOOVER
)
1053 /* whether to print Nose-Hoover chains: */
1054 if (md
->bPrintNHChains
)
1056 if (md
->bNHC_trotter
)
1058 for (i
= 0; (i
< md
->nTC
); i
++)
1060 for (j
= 0; j
< md
->nNHC
; j
++)
1063 md
->tmp_r
[2*k
] = state
->nosehoover_xi
[k
];
1064 md
->tmp_r
[2*k
+1] = state
->nosehoover_vxi
[k
];
1067 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1071 for (i
= 0; (i
< md
->nTCP
); i
++)
1073 for (j
= 0; j
< md
->nNHC
; j
++)
1076 md
->tmp_r
[2*k
] = state
->nhpres_xi
[k
];
1077 md
->tmp_r
[2*k
+1] = state
->nhpres_vxi
[k
];
1080 add_ebin(md
->ebin
, md
->itcb
, md
->mdeb_n
, md
->tmp_r
, bSum
);
1085 for (i
= 0; (i
< md
->nTC
); i
++)
1087 md
->tmp_r
[2*i
] = state
->nosehoover_xi
[i
];
1088 md
->tmp_r
[2*i
+1] = state
->nosehoover_vxi
[i
];
1090 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1094 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
1095 md
->etc
== etcVRESCALE
)
1097 for (i
= 0; (i
< md
->nTC
); i
++)
1099 md
->tmp_r
[i
] = ekind
->tcstat
[i
].lambda
;
1101 add_ebin(md
->ebin
, md
->itc
, md
->nTC
, md
->tmp_r
, bSum
);
1105 if (ekind
&& md
->nU
> 1)
1107 for (i
= 0; (i
< md
->nU
); i
++)
1109 copy_rvec(ekind
->grpstat
[i
].u
, md
->tmp_v
[i
]);
1111 add_ebin(md
->ebin
, md
->iu
, 3*md
->nU
, md
->tmp_v
[0], bSum
);
1114 ebin_increase_count(md
->ebin
, bSum
);
1116 /* BAR + thermodynamic integration values */
1117 if ((md
->fp_dhdl
|| md
->dhc
) && bDoDHDL
)
1119 for (i
= 0; i
< enerd
->n_lambda
-1; i
++)
1121 /* zero for simulated tempering */
1122 md
->dE
[i
] = enerd
->enerpart_lambda
[i
+1]-enerd
->enerpart_lambda
[0];
1123 if (md
->temperatures
!= nullptr)
1125 /* MRS: is this right, given the way we have defined the exchange probabilities? */
1126 /* is this even useful to have at all? */
1127 md
->dE
[i
] += (md
->temperatures
[i
]/
1128 md
->temperatures
[state
->fep_state
]-1.0)*
1129 enerd
->term
[F_EKIN
];
1135 fprintf(md
->fp_dhdl
, "%.4f", time
);
1136 /* the current free energy state */
1138 /* print the current state if we are doing expanded ensemble */
1139 if (expand
->elmcmove
> elmcmoveNO
)
1141 fprintf(md
->fp_dhdl
, " %4d", state
->fep_state
);
1143 /* total energy (for if the temperature changes */
1145 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
1147 switch (fep
->edHdLPrintEnergy
)
1149 case edHdLPrintEnergyPOTENTIAL
:
1150 store_energy
= enerd
->term
[F_EPOT
];
1152 case edHdLPrintEnergyTOTAL
:
1153 case edHdLPrintEnergyYES
:
1155 store_energy
= enerd
->term
[F_ETOT
];
1157 fprintf(md
->fp_dhdl
, " %#.8g", store_energy
);
1160 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
1162 for (i
= 0; i
< efptNR
; i
++)
1164 if (fep
->separate_dvdl
[i
])
1166 /* assumes F_DVDL is first */
1167 fprintf(md
->fp_dhdl
, " %#.8g", enerd
->term
[F_DVDL
+i
]);
1171 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
1173 fprintf(md
->fp_dhdl
, " %#.8g", md
->dE
[i
]);
1177 (md
->epc
!= epcNO
) &&
1178 (enerd
->n_lambda
> 0) &&
1179 (fep
->init_lambda
< 0))
1181 fprintf(md
->fp_dhdl
, " %#.8g", pv
); /* PV term only needed when
1182 there are alternate state
1183 lambda and we're not in
1184 compatibility mode */
1186 fprintf(md
->fp_dhdl
, "\n");
1187 /* and the binary free energy output */
1189 if (md
->dhc
&& bDoDHDL
)
1192 for (i
= 0; i
< efptNR
; i
++)
1194 if (fep
->separate_dvdl
[i
])
1196 /* assumes F_DVDL is first */
1197 store_dhdl
[idhdl
] = enerd
->term
[F_DVDL
+i
];
1201 store_energy
= enerd
->term
[F_ETOT
];
1202 /* store_dh is dE */
1203 mde_delta_h_coll_add_dh(md
->dhc
,
1204 static_cast<double>(state
->fep_state
),
1208 md
->dE
+ fep
->lambda_start_n
,
1215 void upd_mdebin_step(t_mdebin
*md
)
1217 ebin_increase_count(md
->ebin
, FALSE
);
1220 static void npr(FILE *log
, int n
, char c
)
1222 for (; (n
> 0); n
--)
1224 fprintf(log
, "%c", c
);
1228 static void pprint(FILE *log
, const char *s
, t_mdebin
*md
)
1232 char buf1
[22], buf2
[22];
1235 fprintf(log
, "\t<====== ");
1236 npr(log
, slen
, CHAR
);
1237 fprintf(log
, " ==>\n");
1238 fprintf(log
, "\t<==== %s ====>\n", s
);
1239 fprintf(log
, "\t<== ");
1240 npr(log
, slen
, CHAR
);
1241 fprintf(log
, " ======>\n\n");
1243 fprintf(log
, "\tStatistics over %s steps using %s frames\n",
1244 gmx_step_str(md
->ebin
->nsteps_sim
, buf1
),
1245 gmx_step_str(md
->ebin
->nsum_sim
, buf2
));
1249 void print_ebin_header(FILE *log
, int64_t steps
, double time
)
1253 fprintf(log
, " %12s %12s\n"
1255 "Step", "Time", gmx_step_str(steps
, buf
), time
);
1258 // TODO It is too many responsibilities for this function to handle
1259 // both .edr and .log output for both per-time and time-average data.
1260 void print_ebin(ener_file_t fp_ene
, gmx_bool bEne
, gmx_bool bDR
, gmx_bool bOR
,
1262 int64_t step
, double time
,
1264 t_mdebin
*md
, t_fcdata
*fcd
,
1265 gmx_groups_t
*groups
, t_grpopts
*opts
,
1268 /*static char **grpnms=NULL;*/
1270 int i
, j
, n
, ni
, nj
, b
;
1272 real
*disre_rm3tav
, *disre_rt
;
1274 /* these are for the old-style blocks (1 subblock, only reals), because
1275 there can be only one per ID for these */
1282 if (mode
== eprAVER
&& md
->ebin
->nsum_sim
<= 0)
1286 fprintf(log
, "Not enough data recorded to report energy averages\n");
1297 fr
.nsteps
= md
->ebin
->nsteps
;
1298 fr
.dt
= md
->delta_t
;
1299 fr
.nsum
= md
->ebin
->nsum
;
1300 fr
.nre
= (bEne
) ? md
->ebin
->nener
: 0;
1301 fr
.ener
= md
->ebin
->e
;
1302 ndisre
= bDR
? fcd
->disres
.npair
: 0;
1303 disre_rm3tav
= fcd
->disres
.rm3tav
;
1304 disre_rt
= fcd
->disres
.rt
;
1305 /* Optional additional old-style (real-only) blocks. */
1306 for (i
= 0; i
< enxNR
; i
++)
1310 if (fcd
->orires
.nr
> 0 && bOR
)
1312 diagonalize_orires_tensors(&(fcd
->orires
));
1313 nr
[enxOR
] = fcd
->orires
.nr
;
1314 block
[enxOR
] = fcd
->orires
.otav
;
1316 nr
[enxORI
] = (fcd
->orires
.oinsl
!= fcd
->orires
.otav
) ?
1318 block
[enxORI
] = fcd
->orires
.oinsl
;
1319 id
[enxORI
] = enxORI
;
1320 nr
[enxORT
] = fcd
->orires
.nex
*12;
1321 block
[enxORT
] = fcd
->orires
.eig
;
1322 id
[enxORT
] = enxORT
;
1325 /* whether we are going to wrte anything out: */
1326 if (fr
.nre
|| ndisre
|| nr
[enxOR
] || nr
[enxORI
])
1329 /* the old-style blocks go first */
1331 for (i
= 0; i
< enxNR
; i
++)
1338 add_blocks_enxframe(&fr
, fr
.nblock
);
1339 for (b
= 0; b
< fr
.nblock
; b
++)
1341 add_subblocks_enxblock(&(fr
.block
[b
]), 1);
1342 fr
.block
[b
].id
= id
[b
];
1343 fr
.block
[b
].sub
[0].nr
= nr
[b
];
1345 fr
.block
[b
].sub
[0].type
= xdr_datatype_float
;
1346 fr
.block
[b
].sub
[0].fval
= block
[b
];
1348 fr
.block
[b
].sub
[0].type
= xdr_datatype_double
;
1349 fr
.block
[b
].sub
[0].dval
= block
[b
];
1353 /* check for disre block & fill it. */
1358 add_blocks_enxframe(&fr
, fr
.nblock
);
1360 add_subblocks_enxblock(&(fr
.block
[db
]), 2);
1361 fr
.block
[db
].id
= enxDISRE
;
1362 fr
.block
[db
].sub
[0].nr
= ndisre
;
1363 fr
.block
[db
].sub
[1].nr
= ndisre
;
1365 fr
.block
[db
].sub
[0].type
= xdr_datatype_float
;
1366 fr
.block
[db
].sub
[1].type
= xdr_datatype_float
;
1367 fr
.block
[db
].sub
[0].fval
= disre_rt
;
1368 fr
.block
[db
].sub
[1].fval
= disre_rm3tav
;
1370 fr
.block
[db
].sub
[0].type
= xdr_datatype_double
;
1371 fr
.block
[db
].sub
[1].type
= xdr_datatype_double
;
1372 fr
.block
[db
].sub
[0].dval
= disre_rt
;
1373 fr
.block
[db
].sub
[1].dval
= disre_rm3tav
;
1376 /* here we can put new-style blocks */
1378 /* Free energy perturbation blocks */
1381 mde_delta_h_coll_handle_block(md
->dhc
, &fr
, fr
.nblock
);
1384 /* we can now free & reset the data in the blocks */
1387 mde_delta_h_coll_reset(md
->dhc
);
1390 /* AWH bias blocks. */
1391 if (awh
!= nullptr) // TODO: add boolean in t_mdebin. See in mdebin.h.
1393 awh
->writeToEnergyFrame(step
, &fr
);
1396 /* do the actual I/O */
1397 do_enx(fp_ene
, &fr
);
1400 /* We have stored the sums, so reset the sum history */
1401 reset_ebin_sums(md
->ebin
);
1409 pprint(log
, "A V E R A G E S", md
);
1415 pprint(log
, "R M S - F L U C T U A T I O N S", md
);
1419 gmx_fatal(FARGS
, "Invalid print mode (%d)", mode
);
1424 for (i
= 0; i
< opts
->ngtc
; i
++)
1426 if (opts
->annealing
[i
] != eannNO
)
1428 fprintf(log
, "Current ref_t for group %s: %8.1f\n",
1429 *(groups
->grpname
[groups
->grps
[egcTC
].nm_ind
[i
]]),
1433 if (mode
== eprNORMAL
&& fcd
->orires
.nr
> 0)
1435 print_orires_log(log
, &(fcd
->orires
));
1437 fprintf(log
, " Energies (%s)\n", unit_energy
);
1438 pr_ebin(log
, md
->ebin
, md
->ie
, md
->f_nre
+md
->nCrmsd
, 5, mode
, TRUE
);
1441 if (mode
== eprAVER
)
1445 pr_ebin(log
, md
->ebin
, md
->ib
, md
->bTricl
? NTRICLBOXS
: NBOXS
, 5,
1451 fprintf(log
, " Constraint Virial (%s)\n", unit_energy
);
1452 pr_ebin(log
, md
->ebin
, md
->isvir
, 9, 3, mode
, FALSE
);
1454 fprintf(log
, " Force Virial (%s)\n", unit_energy
);
1455 pr_ebin(log
, md
->ebin
, md
->ifvir
, 9, 3, mode
, FALSE
);
1460 fprintf(log
, " Total Virial (%s)\n", unit_energy
);
1461 pr_ebin(log
, md
->ebin
, md
->ivir
, 9, 3, mode
, FALSE
);
1463 fprintf(log
, " Pressure (%s)\n", unit_pres_bar
);
1464 pr_ebin(log
, md
->ebin
, md
->ipres
, 9, 3, mode
, FALSE
);
1469 fprintf(log
, " Total Dipole (%s)\n", unit_dipole_D
);
1470 pr_ebin(log
, md
->ebin
, md
->imu
, 3, 3, mode
, FALSE
);
1476 if (md
->print_grpnms
== nullptr)
1478 snew(md
->print_grpnms
, md
->nE
);
1480 for (i
= 0; (i
< md
->nEg
); i
++)
1482 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
1483 for (j
= i
; (j
< md
->nEg
); j
++)
1485 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
1486 sprintf(buf
, "%s-%s", *(groups
->grpname
[ni
]),
1487 *(groups
->grpname
[nj
]));
1488 md
->print_grpnms
[n
++] = gmx_strdup(buf
);
1492 sprintf(buf
, "Epot (%s)", unit_energy
);
1493 fprintf(log
, "%15s ", buf
);
1494 for (i
= 0; (i
< egNR
); i
++)
1498 fprintf(log
, "%12s ", egrp_nm
[i
]);
1502 for (i
= 0; (i
< md
->nE
); i
++)
1504 fprintf(log
, "%15s", md
->print_grpnms
[i
]);
1505 pr_ebin(log
, md
->ebin
, md
->igrp
[i
], md
->nEc
, md
->nEc
, mode
,
1512 pr_ebin(log
, md
->ebin
, md
->itemp
, md
->nTC
, 4, mode
, TRUE
);
1517 fprintf(log
, "%15s %12s %12s %12s\n",
1518 "Group", "Ux", "Uy", "Uz");
1519 for (i
= 0; (i
< md
->nU
); i
++)
1521 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
1522 fprintf(log
, "%15s", *groups
->grpname
[ni
]);
1523 pr_ebin(log
, md
->ebin
, md
->iu
+3*i
, 3, 3, mode
, FALSE
);
1532 void update_energyhistory(energyhistory_t
* enerhist
, const t_mdebin
* mdebin
)
1534 const t_ebin
* const ebin
= mdebin
->ebin
;
1536 enerhist
->nsteps
= ebin
->nsteps
;
1537 enerhist
->nsum
= ebin
->nsum
;
1538 enerhist
->nsteps_sim
= ebin
->nsteps_sim
;
1539 enerhist
->nsum_sim
= ebin
->nsum_sim
;
1543 /* This will only actually resize the first time */
1544 enerhist
->ener_ave
.resize(ebin
->nener
);
1545 enerhist
->ener_sum
.resize(ebin
->nener
);
1547 for (int i
= 0; i
< ebin
->nener
; i
++)
1549 enerhist
->ener_ave
[i
] = ebin
->e
[i
].eav
;
1550 enerhist
->ener_sum
[i
] = ebin
->e
[i
].esum
;
1554 if (ebin
->nsum_sim
> 0)
1556 /* This will only actually resize the first time */
1557 enerhist
->ener_sum_sim
.resize(ebin
->nener
);
1559 for (int i
= 0; i
< ebin
->nener
; i
++)
1561 enerhist
->ener_sum_sim
[i
] = ebin
->e_sim
[i
].esum
;
1566 mde_delta_h_coll_update_energyhistory(mdebin
->dhc
, enerhist
);
1570 void restore_energyhistory_from_state(t_mdebin
* mdebin
,
1571 const energyhistory_t
* enerhist
)
1573 unsigned int nener
= static_cast<unsigned int>(mdebin
->ebin
->nener
);
1575 GMX_RELEASE_ASSERT(enerhist
, "Need valid history to restore");
1577 if ((enerhist
->nsum
> 0 && nener
!= enerhist
->ener_sum
.size()) ||
1578 (enerhist
->nsum_sim
> 0 && nener
!= enerhist
->ener_sum_sim
.size()))
1580 gmx_fatal(FARGS
, "Mismatch between number of energies in run input (%u) and checkpoint file (%zu or %zu).",
1581 nener
, enerhist
->ener_sum
.size(), enerhist
->ener_sum_sim
.size());
1584 mdebin
->ebin
->nsteps
= enerhist
->nsteps
;
1585 mdebin
->ebin
->nsum
= enerhist
->nsum
;
1586 mdebin
->ebin
->nsteps_sim
= enerhist
->nsteps_sim
;
1587 mdebin
->ebin
->nsum_sim
= enerhist
->nsum_sim
;
1589 for (int i
= 0; i
< mdebin
->ebin
->nener
; i
++)
1591 mdebin
->ebin
->e
[i
].eav
=
1592 (enerhist
->nsum
> 0 ? enerhist
->ener_ave
[i
] : 0);
1593 mdebin
->ebin
->e
[i
].esum
=
1594 (enerhist
->nsum
> 0 ? enerhist
->ener_sum
[i
] : 0);
1595 mdebin
->ebin
->e_sim
[i
].esum
=
1596 (enerhist
->nsum_sim
> 0 ? enerhist
->ener_sum_sim
[i
] : 0);
1600 mde_delta_h_coll_restore_energyhistory(mdebin
->dhc
, enerhist
->deltaHForeignLambdas
.get());