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, 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.
45 #include "gromacs/fileio/enxio.h"
46 #include "gromacs/fileio/gmxfio.h"
47 #include "gromacs/fileio/xvgr.h"
48 #include "gromacs/gmxlib/network.h"
49 #include "gromacs/listed-forces/disre.h"
50 #include "gromacs/listed-forces/orires.h"
51 #include "gromacs/math/functions.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdlib/constr.h"
55 #include "gromacs/mdlib/mdebin_bar.h"
56 #include "gromacs/mdlib/mdrun.h"
57 #include "gromacs/mdtypes/energyhistory.h"
58 #include "gromacs/mdtypes/fcdata.h"
59 #include "gromacs/mdtypes/group.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/pulling/pull.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/utility/arraysize.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/smalloc.h"
69 static const char *conrmsd_nm
[] = { "Constr. rmsd", "Constr.2 rmsd" };
71 static const char *boxs_nm
[] = { "Box-X", "Box-Y", "Box-Z" };
73 static const char *tricl_boxs_nm
[] = {
74 "Box-XX", "Box-YY", "Box-ZZ",
75 "Box-YX", "Box-ZX", "Box-ZY"
78 static const char *vol_nm
[] = { "Volume" };
80 static const char *dens_nm
[] = {"Density" };
82 static const char *pv_nm
[] = {"pV" };
84 static const char *enthalpy_nm
[] = {"Enthalpy" };
86 static const char *boxvel_nm
[] = {
87 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
88 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
91 #define NBOXS asize(boxs_nm)
92 #define NTRICLBOXS asize(tricl_boxs_nm)
94 t_mdebin
*init_mdebin(ener_file_t fp_ene
,
95 const gmx_mtop_t
*mtop
,
99 const char *ener_nm
[F_NRE
];
100 static const char *vir_nm
[] = {
101 "Vir-XX", "Vir-XY", "Vir-XZ",
102 "Vir-YX", "Vir-YY", "Vir-YZ",
103 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
105 static const char *sv_nm
[] = {
106 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
107 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
108 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
110 static const char *fv_nm
[] = {
111 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
112 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
113 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
115 static const char *pres_nm
[] = {
116 "Pres-XX", "Pres-XY", "Pres-XZ",
117 "Pres-YX", "Pres-YY", "Pres-YZ",
118 "Pres-ZX", "Pres-ZY", "Pres-ZZ"
120 static const char *surft_nm
[] = {
123 static const char *mu_nm
[] = {
124 "Mu-X", "Mu-Y", "Mu-Z"
126 static const char *vcos_nm
[] = {
129 static const char *visc_nm
[] = {
132 static const char *baro_nm
[] = {
137 const gmx_groups_t
*groups
;
142 int i
, j
, ni
, nj
, n
, k
, kk
, ncon
, nset
;
147 if (EI_DYNAMICS(ir
->eI
))
149 md
->delta_t
= ir
->delta_t
;
156 groups
= &mtop
->groups
;
158 bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
159 b14
= (gmx_mtop_ftype_count(mtop
, F_LJ14
) > 0 ||
160 gmx_mtop_ftype_count(mtop
, F_LJC14_Q
) > 0);
162 ncon
= gmx_mtop_ftype_count(mtop
, F_CONSTR
);
163 nset
= gmx_mtop_ftype_count(mtop
, F_SETTLE
);
164 md
->bConstr
= (ncon
> 0 || nset
> 0);
165 md
->bConstrVir
= FALSE
;
168 if (ncon
> 0 && ir
->eConstrAlg
== econtLINCS
)
172 md
->bConstrVir
= (getenv("GMX_CONSTRAINTVIR") != nullptr);
179 /* Energy monitoring */
180 for (i
= 0; i
< egNR
; i
++)
182 md
->bEInd
[i
] = FALSE
;
185 for (i
= 0; i
< F_NRE
; i
++)
187 md
->bEner
[i
] = FALSE
;
190 md
->bEner
[i
] = !bBHAM
;
192 else if (i
== F_BHAM
)
194 md
->bEner
[i
] = bBHAM
;
198 md
->bEner
[i
] = ir
->bQMMM
;
200 else if (i
== F_RF_EXCL
)
202 md
->bEner
[i
] = (EEL_RF(ir
->coulombtype
) && ir
->cutoff_scheme
== ecutsGROUP
);
204 else if (i
== F_COUL_RECIP
)
206 md
->bEner
[i
] = EEL_FULL(ir
->coulombtype
);
208 else if (i
== F_LJ_RECIP
)
210 md
->bEner
[i
] = EVDW_PME(ir
->vdwtype
);
212 else if (i
== F_LJ14
)
216 else if (i
== F_COUL14
)
220 else if (i
== F_LJC14_Q
|| i
== F_LJC_PAIRS_NB
)
222 md
->bEner
[i
] = FALSE
;
224 else if ((i
== F_DVDL_COUL
&& ir
->fepvals
->separate_dvdl
[efptCOUL
]) ||
225 (i
== F_DVDL_VDW
&& ir
->fepvals
->separate_dvdl
[efptVDW
]) ||
226 (i
== F_DVDL_BONDED
&& ir
->fepvals
->separate_dvdl
[efptBONDED
]) ||
227 (i
== F_DVDL_RESTRAINT
&& ir
->fepvals
->separate_dvdl
[efptRESTRAINT
]) ||
228 (i
== F_DKDL
&& ir
->fepvals
->separate_dvdl
[efptMASS
]) ||
229 (i
== F_DVDL
&& ir
->fepvals
->separate_dvdl
[efptFEP
]))
231 md
->bEner
[i
] = (ir
->efep
!= efepNO
);
233 else if ((interaction_function
[i
].flags
& IF_VSITE
) ||
234 (i
== F_CONSTR
) || (i
== F_CONSTRNC
) || (i
== F_SETTLE
))
236 md
->bEner
[i
] = FALSE
;
238 else if ((i
== F_COUL_SR
) || (i
== F_EPOT
) || (i
== F_PRES
) || (i
== F_EQM
))
242 else if ((i
== F_GBPOL
) && ir
->implicit_solvent
== eisGBSA
)
246 else if ((i
== F_NPSOLVATION
) && ir
->implicit_solvent
== eisGBSA
&& (ir
->sa_algorithm
!= esaNO
))
250 else if ((i
== F_GB12
) || (i
== F_GB13
) || (i
== F_GB14
))
252 md
->bEner
[i
] = FALSE
;
254 else if ((i
== F_ETOT
) || (i
== F_EKIN
) || (i
== F_TEMP
))
256 md
->bEner
[i
] = EI_DYNAMICS(ir
->eI
);
258 else if (i
== F_DISPCORR
|| i
== F_PDISPCORR
)
260 md
->bEner
[i
] = (ir
->eDispCorr
!= edispcNO
);
262 else if (i
== F_DISRESVIOL
)
264 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_DISRES
) > 0);
266 else if (i
== F_ORIRESDEV
)
268 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_ORIRES
) > 0);
270 else if (i
== F_CONNBONDS
)
272 md
->bEner
[i
] = FALSE
;
274 else if (i
== F_COM_PULL
)
276 md
->bEner
[i
] = (ir
->bPull
&& pull_have_potential(ir
->pull_work
));
278 else if (i
== F_ECONSERVED
)
280 md
->bEner
[i
] = ((ir
->etc
== etcNOSEHOOVER
|| ir
->etc
== etcVRESCALE
) &&
281 (ir
->epc
== epcNO
|| ir
->epc
== epcMTTK
));
285 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, i
) > 0);
290 for (i
= 0; i
< F_NRE
; i
++)
294 ener_nm
[md
->f_nre
] = interaction_function
[i
].longname
;
300 md
->bDiagPres
= !TRICLINIC(ir
->ref_p
);
301 md
->ref_p
= (ir
->ref_p
[XX
][XX
]+ir
->ref_p
[YY
][YY
]+ir
->ref_p
[ZZ
][ZZ
])/DIM
;
302 md
->bTricl
= TRICLINIC(ir
->compress
) || TRICLINIC(ir
->deform
);
303 md
->bDynBox
= inputrecDynamicBox(ir
);
305 md
->bNHC_trotter
= inputrecNvtTrotter(ir
);
306 md
->bPrintNHChains
= ir
->bPrintNHChains
;
307 md
->bMTTK
= (inputrecNptTrotter(ir
) || inputrecNphTrotter(ir
));
308 md
->bMu
= inputrecNeedMutot(ir
);
310 md
->ebin
= mk_ebin();
311 /* Pass NULL for unit to let get_ebin_space determine the units
312 * for interaction_function[i].longname
314 md
->ie
= get_ebin_space(md
->ebin
, md
->f_nre
, ener_nm
, nullptr);
317 /* This should be called directly after the call for md->ie,
318 * such that md->iconrmsd follows directly in the list.
320 md
->iconrmsd
= get_ebin_space(md
->ebin
, md
->nCrmsd
, conrmsd_nm
, "");
324 md
->ib
= get_ebin_space(md
->ebin
,
325 md
->bTricl
? NTRICLBOXS
: NBOXS
,
326 md
->bTricl
? tricl_boxs_nm
: boxs_nm
,
328 md
->ivol
= get_ebin_space(md
->ebin
, 1, vol_nm
, unit_volume
);
329 md
->idens
= get_ebin_space(md
->ebin
, 1, dens_nm
, unit_density_SI
);
332 md
->ipv
= get_ebin_space(md
->ebin
, 1, pv_nm
, unit_energy
);
333 md
->ienthalpy
= get_ebin_space(md
->ebin
, 1, enthalpy_nm
, unit_energy
);
338 md
->isvir
= get_ebin_space(md
->ebin
, asize(sv_nm
), sv_nm
, unit_energy
);
339 md
->ifvir
= get_ebin_space(md
->ebin
, asize(fv_nm
), fv_nm
, unit_energy
);
341 md
->ivir
= get_ebin_space(md
->ebin
, asize(vir_nm
), vir_nm
, unit_energy
);
342 md
->ipres
= get_ebin_space(md
->ebin
, asize(pres_nm
), pres_nm
, unit_pres_bar
);
343 md
->isurft
= get_ebin_space(md
->ebin
, asize(surft_nm
), surft_nm
,
345 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
347 md
->ipc
= get_ebin_space(md
->ebin
, md
->bTricl
? 6 : 3,
348 boxvel_nm
, unit_vel
);
352 md
->imu
= get_ebin_space(md
->ebin
, asize(mu_nm
), mu_nm
, unit_dipole_D
);
354 if (ir
->cos_accel
!= 0)
356 md
->ivcos
= get_ebin_space(md
->ebin
, asize(vcos_nm
), vcos_nm
, unit_vel
);
357 md
->ivisc
= get_ebin_space(md
->ebin
, asize(visc_nm
), visc_nm
,
361 /* Energy monitoring */
362 for (i
= 0; i
< egNR
; i
++)
364 md
->bEInd
[i
] = FALSE
;
366 md
->bEInd
[egCOULSR
] = TRUE
;
367 md
->bEInd
[egLJSR
] = TRUE
;
371 md
->bEInd
[egLJSR
] = FALSE
;
372 md
->bEInd
[egBHAMSR
] = TRUE
;
376 md
->bEInd
[egLJ14
] = TRUE
;
377 md
->bEInd
[egCOUL14
] = TRUE
;
380 for (i
= 0; (i
< egNR
); i
++)
388 n
= groups
->grps
[egcENER
].nr
;
390 md
->nE
= (n
*(n
+1))/2;
392 snew(md
->igrp
, md
->nE
);
397 for (k
= 0; (k
< md
->nEc
); k
++)
399 snew(gnm
[k
], STRLEN
);
401 for (i
= 0; (i
< groups
->grps
[egcENER
].nr
); i
++)
403 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
404 for (j
= i
; (j
< groups
->grps
[egcENER
].nr
); j
++)
406 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
407 for (k
= kk
= 0; (k
< egNR
); k
++)
411 sprintf(gnm
[kk
], "%s:%s-%s", egrp_nm
[k
],
412 *(groups
->grpname
[ni
]), *(groups
->grpname
[nj
]));
416 md
->igrp
[n
] = get_ebin_space(md
->ebin
, md
->nEc
,
417 (const char **)gnm
, unit_energy
);
421 for (k
= 0; (k
< md
->nEc
); k
++)
429 gmx_incons("Number of energy terms wrong");
433 md
->nTC
= groups
->grps
[egcTC
].nr
;
434 md
->nNHC
= ir
->opts
.nhchainlength
; /* shorthand for number of NH chains */
437 md
->nTCP
= 1; /* assume only one possible coupling system for barostat
444 if (md
->etc
== etcNOSEHOOVER
)
446 if (md
->bNHC_trotter
)
448 md
->mde_n
= 2*md
->nNHC
*md
->nTC
;
452 md
->mde_n
= 2*md
->nTC
;
454 if (md
->epc
== epcMTTK
)
456 md
->mdeb_n
= 2*md
->nNHC
*md
->nTCP
;
465 snew(md
->tmp_r
, md
->mde_n
);
466 snew(md
->tmp_v
, md
->mde_n
);
467 snew(md
->grpnms
, md
->mde_n
);
470 for (i
= 0; (i
< md
->nTC
); i
++)
472 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
473 sprintf(buf
, "T-%s", *(groups
->grpname
[ni
]));
474 grpnms
[i
] = gmx_strdup(buf
);
476 md
->itemp
= get_ebin_space(md
->ebin
, md
->nTC
, (const char **)grpnms
,
479 if (md
->etc
== etcNOSEHOOVER
)
481 if (md
->bPrintNHChains
)
483 if (md
->bNHC_trotter
)
485 for (i
= 0; (i
< md
->nTC
); i
++)
487 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
488 bufi
= *(groups
->grpname
[ni
]);
489 for (j
= 0; (j
< md
->nNHC
); j
++)
491 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
492 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
493 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
494 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
497 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
498 (const char **)grpnms
, unit_invtime
);
501 for (i
= 0; (i
< md
->nTCP
); i
++)
503 bufi
= baro_nm
[0]; /* All barostat DOF's together for now. */
504 for (j
= 0; (j
< md
->nNHC
); j
++)
506 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
507 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
508 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
509 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
512 md
->itcb
= get_ebin_space(md
->ebin
, md
->mdeb_n
,
513 (const char **)grpnms
, unit_invtime
);
518 for (i
= 0; (i
< md
->nTC
); i
++)
520 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
521 bufi
= *(groups
->grpname
[ni
]);
522 sprintf(buf
, "Xi-%s", bufi
);
523 grpnms
[2*i
] = gmx_strdup(buf
);
524 sprintf(buf
, "vXi-%s", bufi
);
525 grpnms
[2*i
+1] = gmx_strdup(buf
);
527 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
528 (const char **)grpnms
, unit_invtime
);
532 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
533 md
->etc
== etcVRESCALE
)
535 for (i
= 0; (i
< md
->nTC
); i
++)
537 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
538 sprintf(buf
, "Lamb-%s", *(groups
->grpname
[ni
]));
539 grpnms
[i
] = gmx_strdup(buf
);
541 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
, (const char **)grpnms
, "");
547 md
->nU
= groups
->grps
[egcACC
].nr
;
550 snew(grpnms
, 3*md
->nU
);
551 for (i
= 0; (i
< md
->nU
); i
++)
553 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
554 sprintf(buf
, "Ux-%s", *(groups
->grpname
[ni
]));
555 grpnms
[3*i
+XX
] = gmx_strdup(buf
);
556 sprintf(buf
, "Uy-%s", *(groups
->grpname
[ni
]));
557 grpnms
[3*i
+YY
] = gmx_strdup(buf
);
558 sprintf(buf
, "Uz-%s", *(groups
->grpname
[ni
]));
559 grpnms
[3*i
+ZZ
] = gmx_strdup(buf
);
561 md
->iu
= get_ebin_space(md
->ebin
, 3*md
->nU
, (const char **)grpnms
, unit_vel
);
567 do_enxnms(fp_ene
, &md
->ebin
->nener
, &md
->ebin
->enm
);
570 md
->print_grpnms
= nullptr;
572 /* check whether we're going to write dh histograms */
574 if (ir
->fepvals
->separate_dhdl_file
== esepdhdlfileNO
)
576 /* Currently dh histograms are only written with dynamics */
577 if (EI_DYNAMICS(ir
->eI
))
581 mde_delta_h_coll_init(md
->dhc
, ir
);
583 md
->fp_dhdl
= nullptr;
584 snew(md
->dE
, ir
->fepvals
->n_lambda
);
588 md
->fp_dhdl
= fp_dhdl
;
589 snew(md
->dE
, ir
->fepvals
->n_lambda
);
594 snew(md
->temperatures
, ir
->fepvals
->n_lambda
);
595 for (i
= 0; i
< ir
->fepvals
->n_lambda
; i
++)
597 md
->temperatures
[i
] = ir
->simtempvals
->temperatures
[i
];
603 /* print a lambda vector to a string
604 fep = the inputrec's FEP input data
605 i = the index of the lambda vector
606 get_native_lambda = whether to print the native lambda
607 get_names = whether to print the names rather than the values
608 str = the pre-allocated string buffer to print to. */
609 static void print_lambda_vector(t_lambda
*fep
, int i
,
610 gmx_bool get_native_lambda
, gmx_bool get_names
,
616 for (j
= 0; j
< efptNR
; j
++)
618 if (fep
->separate_dvdl
[j
])
623 str
[0] = 0; /* reset the string */
626 str
+= sprintf(str
, "("); /* set the opening parenthesis*/
628 for (j
= 0; j
< efptNR
; j
++)
630 if (fep
->separate_dvdl
[j
])
634 if (get_native_lambda
&& fep
->init_lambda
>= 0)
636 str
+= sprintf(str
, "%.4f", fep
->init_lambda
);
640 str
+= sprintf(str
, "%.4f", fep
->all_lambda
[j
][i
]);
645 str
+= sprintf(str
, "%s", efpt_singular_names
[j
]);
647 /* print comma for the next item */
650 str
+= sprintf(str
, ", ");
657 /* and add the closing parenthesis */
663 extern FILE *open_dhdl(const char *filename
, const t_inputrec
*ir
,
664 const gmx_output_env_t
*oenv
)
667 const char *dhdl
= "dH/d\\lambda", *deltag
= "\\DeltaH", *lambda
= "\\lambda",
668 *lambdastate
= "\\lambda state";
669 char title
[STRLEN
], label_x
[STRLEN
], label_y
[STRLEN
];
670 int i
, nps
, nsets
, nsets_de
, nsetsbegin
;
671 int n_lambda_terms
= 0;
672 t_lambda
*fep
= ir
->fepvals
; /* for simplicity */
673 t_expanded
*expand
= ir
->expandedvals
;
675 char buf
[STRLEN
], lambda_vec_str
[STRLEN
], lambda_name_str
[STRLEN
];
681 gmx_bool write_pV
= FALSE
;
683 /* count the number of different lambda terms */
684 for (i
= 0; i
< efptNR
; i
++)
686 if (fep
->separate_dvdl
[i
])
692 if (fep
->n_lambda
== 0)
694 sprintf(title
, "%s", dhdl
);
695 sprintf(label_x
, "Time (ps)");
696 sprintf(label_y
, "%s (%s %s)",
697 dhdl
, unit_energy
, "[\\lambda]\\S-1\\N");
701 sprintf(title
, "%s and %s", dhdl
, deltag
);
702 sprintf(label_x
, "Time (ps)");
703 sprintf(label_y
, "%s and %s (%s %s)",
704 dhdl
, deltag
, unit_energy
, "[\\8l\\4]\\S-1\\N");
706 fp
= gmx_fio_fopen(filename
, "w+");
707 xvgr_header(fp
, title
, label_x
, label_y
, exvggtXNY
, oenv
);
711 bufplace
= sprintf(buf
, "T = %g (K) ",
714 if ((ir
->efep
!= efepSLOWGROWTH
) && (ir
->efep
!= efepEXPANDED
))
716 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
718 /* compatibility output */
719 sprintf(&(buf
[bufplace
]), "%s = %.4f", lambda
, fep
->init_lambda
);
723 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, FALSE
,
725 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, TRUE
,
727 sprintf(&(buf
[bufplace
]), "%s %d: %s = %s",
728 lambdastate
, fep
->init_fep_state
,
729 lambda_name_str
, lambda_vec_str
);
732 xvgr_subtitle(fp
, buf
, oenv
);
736 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
738 nsets_dhdl
= n_lambda_terms
;
740 /* count the number of delta_g states */
741 nsets_de
= fep
->lambda_stop_n
- fep
->lambda_start_n
;
743 nsets
= nsets_dhdl
+ nsets_de
; /* dhdl + fep differences */
745 if (fep
->n_lambda
> 0 && (expand
->elmcmove
> elmcmoveNO
))
747 nsets
+= 1; /*add fep state for expanded ensemble */
750 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
752 nsets
+= 1; /* add energy to the dhdl as well */
756 if ((ir
->epc
!= epcNO
) && (fep
->n_lambda
> 0) && (fep
->init_lambda
< 0))
758 nsetsextend
+= 1; /* for PV term, other terms possible if required for
759 the reduced potential (only needed with foreign
760 lambda, and only output when init_lambda is not
761 set in order to maintain compatibility of the
765 snew(setname
, nsetsextend
);
767 if (expand
->elmcmove
> elmcmoveNO
)
769 /* state for the fep_vals, if we have alchemical sampling */
770 sprintf(buf
, "%s", "Thermodynamic state");
771 setname
[s
] = gmx_strdup(buf
);
775 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
777 switch (fep
->edHdLPrintEnergy
)
779 case edHdLPrintEnergyPOTENTIAL
:
780 sprintf(buf
, "%s (%s)", "Potential Energy", unit_energy
);
782 case edHdLPrintEnergyTOTAL
:
783 case edHdLPrintEnergyYES
:
785 sprintf(buf
, "%s (%s)", "Total Energy", unit_energy
);
787 setname
[s
] = gmx_strdup(buf
);
791 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
793 for (i
= 0; i
< efptNR
; i
++)
795 if (fep
->separate_dvdl
[i
])
798 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
800 /* compatibility output */
801 sprintf(buf
, "%s %s %.4f", dhdl
, lambda
, fep
->init_lambda
);
805 double lam
= fep
->init_lambda
;
806 if (fep
->init_lambda
< 0)
808 lam
= fep
->all_lambda
[i
][fep
->init_fep_state
];
810 sprintf(buf
, "%s %s = %.4f", dhdl
, efpt_singular_names
[i
],
813 setname
[s
] = gmx_strdup(buf
);
819 if (fep
->n_lambda
> 0)
821 /* g_bar has to determine the lambda values used in this simulation
822 * from this xvg legend.
825 if (expand
->elmcmove
> elmcmoveNO
)
827 nsetsbegin
= 1; /* for including the expanded ensemble */
834 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
838 nsetsbegin
+= nsets_dhdl
;
840 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
842 print_lambda_vector(fep
, i
, FALSE
, FALSE
, lambda_vec_str
);
843 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
845 /* for compatible dhdl.xvg files */
846 nps
= sprintf(buf
, "%s %s %s", deltag
, lambda
, lambda_vec_str
);
850 nps
= sprintf(buf
, "%s %s to %s", deltag
, lambda
, lambda_vec_str
);
855 /* print the temperature for this state if doing simulated annealing */
856 sprintf(&buf
[nps
], "T = %g (%s)",
857 ir
->simtempvals
->temperatures
[s
-(nsetsbegin
)],
860 setname
[s
] = gmx_strdup(buf
);
865 sprintf(buf
, "pV (%s)", unit_energy
);
866 setname
[nsetsextend
-1] = gmx_strdup(buf
); /* the first entry after
870 xvgr_legend(fp
, nsetsextend
, (const char **)setname
, oenv
);
872 for (s
= 0; s
< nsetsextend
; s
++)
882 static void copy_energy(t_mdebin
*md
, real e
[], real ecpy
[])
886 for (i
= j
= 0; (i
< F_NRE
); i
++)
895 gmx_incons("Number of energy terms wrong");
899 void upd_mdebin(t_mdebin
*md
,
904 gmx_enerdata_t
*enerd
,
913 gmx_ekindata_t
*ekind
,
917 int i
, j
, k
, kk
, n
, gid
;
918 real crmsd
[2], tmp6
[6];
919 real bs
[NTRICLBOXS
], vol
, dens
, pv
, enthalpy
;
922 double store_dhdl
[efptNR
];
923 real store_energy
= 0;
926 /* Do NOT use the box in the state variable, but the separate box provided
927 * as an argument. This is because we sometimes need to write the box from
928 * the last timestep to match the trajectory frames.
930 copy_energy(md
, enerd
->term
, ecopy
);
931 add_ebin(md
->ebin
, md
->ie
, md
->f_nre
, ecopy
, bSum
);
934 crmsd
[0] = constr_rmsd(constr
);
935 add_ebin(md
->ebin
, md
->iconrmsd
, md
->nCrmsd
, crmsd
, FALSE
);
957 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
958 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
959 add_ebin(md
->ebin
, md
->ib
, nboxs
, bs
, bSum
);
960 add_ebin(md
->ebin
, md
->ivol
, 1, &vol
, bSum
);
961 add_ebin(md
->ebin
, md
->idens
, 1, &dens
, bSum
);
965 /* This is pV (in kJ/mol). The pressure is the reference pressure,
966 not the instantaneous pressure */
967 pv
= vol
*md
->ref_p
/PRESFAC
;
969 add_ebin(md
->ebin
, md
->ipv
, 1, &pv
, bSum
);
970 enthalpy
= pv
+ enerd
->term
[F_ETOT
];
971 add_ebin(md
->ebin
, md
->ienthalpy
, 1, &enthalpy
, bSum
);
976 add_ebin(md
->ebin
, md
->isvir
, 9, svir
[0], bSum
);
977 add_ebin(md
->ebin
, md
->ifvir
, 9, fvir
[0], bSum
);
979 add_ebin(md
->ebin
, md
->ivir
, 9, vir
[0], bSum
);
980 add_ebin(md
->ebin
, md
->ipres
, 9, pres
[0], bSum
);
981 tmp
= (pres
[ZZ
][ZZ
]-(pres
[XX
][XX
]+pres
[YY
][YY
])*0.5)*box
[ZZ
][ZZ
];
982 add_ebin(md
->ebin
, md
->isurft
, 1, &tmp
, bSum
);
983 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
985 tmp6
[0] = state
->boxv
[XX
][XX
];
986 tmp6
[1] = state
->boxv
[YY
][YY
];
987 tmp6
[2] = state
->boxv
[ZZ
][ZZ
];
988 tmp6
[3] = state
->boxv
[YY
][XX
];
989 tmp6
[4] = state
->boxv
[ZZ
][XX
];
990 tmp6
[5] = state
->boxv
[ZZ
][YY
];
991 add_ebin(md
->ebin
, md
->ipc
, md
->bTricl
? 6 : 3, tmp6
, bSum
);
995 add_ebin(md
->ebin
, md
->imu
, 3, mu_tot
, bSum
);
997 if (ekind
&& ekind
->cosacc
.cos_accel
!= 0)
999 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
1000 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
1001 add_ebin(md
->ebin
, md
->ivcos
, 1, &(ekind
->cosacc
.vcos
), bSum
);
1002 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
1003 tmp
= 1/(ekind
->cosacc
.cos_accel
/(ekind
->cosacc
.vcos
*PICO
)
1004 *dens
*gmx::square(box
[ZZ
][ZZ
]*NANO
/(2*M_PI
)));
1005 add_ebin(md
->ebin
, md
->ivisc
, 1, &tmp
, bSum
);
1010 for (i
= 0; (i
< md
->nEg
); i
++)
1012 for (j
= i
; (j
< md
->nEg
); j
++)
1014 gid
= GID(i
, j
, md
->nEg
);
1015 for (k
= kk
= 0; (k
< egNR
); k
++)
1019 eee
[kk
++] = enerd
->grpp
.ener
[k
][gid
];
1022 add_ebin(md
->ebin
, md
->igrp
[n
], md
->nEc
, eee
, bSum
);
1030 for (i
= 0; (i
< md
->nTC
); i
++)
1032 md
->tmp_r
[i
] = ekind
->tcstat
[i
].T
;
1034 add_ebin(md
->ebin
, md
->itemp
, md
->nTC
, md
->tmp_r
, bSum
);
1036 if (md
->etc
== etcNOSEHOOVER
)
1038 /* whether to print Nose-Hoover chains: */
1039 if (md
->bPrintNHChains
)
1041 if (md
->bNHC_trotter
)
1043 for (i
= 0; (i
< md
->nTC
); i
++)
1045 for (j
= 0; j
< md
->nNHC
; j
++)
1048 md
->tmp_r
[2*k
] = state
->nosehoover_xi
[k
];
1049 md
->tmp_r
[2*k
+1] = state
->nosehoover_vxi
[k
];
1052 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1056 for (i
= 0; (i
< md
->nTCP
); i
++)
1058 for (j
= 0; j
< md
->nNHC
; j
++)
1061 md
->tmp_r
[2*k
] = state
->nhpres_xi
[k
];
1062 md
->tmp_r
[2*k
+1] = state
->nhpres_vxi
[k
];
1065 add_ebin(md
->ebin
, md
->itcb
, md
->mdeb_n
, md
->tmp_r
, bSum
);
1070 for (i
= 0; (i
< md
->nTC
); i
++)
1072 md
->tmp_r
[2*i
] = state
->nosehoover_xi
[i
];
1073 md
->tmp_r
[2*i
+1] = state
->nosehoover_vxi
[i
];
1075 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1079 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
1080 md
->etc
== etcVRESCALE
)
1082 for (i
= 0; (i
< md
->nTC
); i
++)
1084 md
->tmp_r
[i
] = ekind
->tcstat
[i
].lambda
;
1086 add_ebin(md
->ebin
, md
->itc
, md
->nTC
, md
->tmp_r
, bSum
);
1090 if (ekind
&& md
->nU
> 1)
1092 for (i
= 0; (i
< md
->nU
); i
++)
1094 copy_rvec(ekind
->grpstat
[i
].u
, md
->tmp_v
[i
]);
1096 add_ebin(md
->ebin
, md
->iu
, 3*md
->nU
, md
->tmp_v
[0], bSum
);
1099 ebin_increase_count(md
->ebin
, bSum
);
1101 /* BAR + thermodynamic integration values */
1102 if ((md
->fp_dhdl
|| md
->dhc
) && bDoDHDL
)
1104 for (i
= 0; i
< enerd
->n_lambda
-1; i
++)
1106 /* zero for simulated tempering */
1107 md
->dE
[i
] = enerd
->enerpart_lambda
[i
+1]-enerd
->enerpart_lambda
[0];
1108 if (md
->temperatures
!= nullptr)
1110 /* MRS: is this right, given the way we have defined the exchange probabilities? */
1111 /* is this even useful to have at all? */
1112 md
->dE
[i
] += (md
->temperatures
[i
]/
1113 md
->temperatures
[state
->fep_state
]-1.0)*
1114 enerd
->term
[F_EKIN
];
1120 fprintf(md
->fp_dhdl
, "%.4f", time
);
1121 /* the current free energy state */
1123 /* print the current state if we are doing expanded ensemble */
1124 if (expand
->elmcmove
> elmcmoveNO
)
1126 fprintf(md
->fp_dhdl
, " %4d", state
->fep_state
);
1128 /* total energy (for if the temperature changes */
1130 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
1132 switch (fep
->edHdLPrintEnergy
)
1134 case edHdLPrintEnergyPOTENTIAL
:
1135 store_energy
= enerd
->term
[F_EPOT
];
1137 case edHdLPrintEnergyTOTAL
:
1138 case edHdLPrintEnergyYES
:
1140 store_energy
= enerd
->term
[F_ETOT
];
1142 fprintf(md
->fp_dhdl
, " %#.8g", store_energy
);
1145 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
1147 for (i
= 0; i
< efptNR
; i
++)
1149 if (fep
->separate_dvdl
[i
])
1151 /* assumes F_DVDL is first */
1152 fprintf(md
->fp_dhdl
, " %#.8g", enerd
->term
[F_DVDL
+i
]);
1156 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
1158 fprintf(md
->fp_dhdl
, " %#.8g", md
->dE
[i
]);
1162 (md
->epc
!= epcNO
) &&
1163 (enerd
->n_lambda
> 0) &&
1164 (fep
->init_lambda
< 0))
1166 fprintf(md
->fp_dhdl
, " %#.8g", pv
); /* PV term only needed when
1167 there are alternate state
1168 lambda and we're not in
1169 compatibility mode */
1171 fprintf(md
->fp_dhdl
, "\n");
1172 /* and the binary free energy output */
1174 if (md
->dhc
&& bDoDHDL
)
1177 for (i
= 0; i
< efptNR
; i
++)
1179 if (fep
->separate_dvdl
[i
])
1181 /* assumes F_DVDL is first */
1182 store_dhdl
[idhdl
] = enerd
->term
[F_DVDL
+i
];
1186 store_energy
= enerd
->term
[F_ETOT
];
1187 /* store_dh is dE */
1188 mde_delta_h_coll_add_dh(md
->dhc
,
1189 (double)state
->fep_state
,
1193 md
->dE
+ fep
->lambda_start_n
,
1200 void upd_mdebin_step(t_mdebin
*md
)
1202 ebin_increase_count(md
->ebin
, FALSE
);
1205 static void npr(FILE *log
, int n
, char c
)
1207 for (; (n
> 0); n
--)
1209 fprintf(log
, "%c", c
);
1213 static void pprint(FILE *log
, const char *s
, t_mdebin
*md
)
1217 char buf1
[22], buf2
[22];
1220 fprintf(log
, "\t<====== ");
1221 npr(log
, slen
, CHAR
);
1222 fprintf(log
, " ==>\n");
1223 fprintf(log
, "\t<==== %s ====>\n", s
);
1224 fprintf(log
, "\t<== ");
1225 npr(log
, slen
, CHAR
);
1226 fprintf(log
, " ======>\n\n");
1228 fprintf(log
, "\tStatistics over %s steps using %s frames\n",
1229 gmx_step_str(md
->ebin
->nsteps_sim
, buf1
),
1230 gmx_step_str(md
->ebin
->nsum_sim
, buf2
));
1234 void print_ebin_header(FILE *log
, gmx_int64_t steps
, double time
)
1238 fprintf(log
, " %12s %12s\n"
1240 "Step", "Time", gmx_step_str(steps
, buf
), time
);
1243 void print_ebin(ener_file_t fp_ene
, gmx_bool bEne
, gmx_bool bDR
, gmx_bool bOR
,
1245 gmx_int64_t step
, double time
,
1247 t_mdebin
*md
, t_fcdata
*fcd
,
1248 gmx_groups_t
*groups
, t_grpopts
*opts
)
1250 /*static char **grpnms=NULL;*/
1252 int i
, j
, n
, ni
, nj
, b
;
1254 real
*disre_rm3tav
, *disre_rt
;
1256 /* these are for the old-style blocks (1 subblock, only reals), because
1257 there can be only one per ID for these */
1270 fr
.nsteps
= md
->ebin
->nsteps
;
1271 fr
.dt
= md
->delta_t
;
1272 fr
.nsum
= md
->ebin
->nsum
;
1273 fr
.nre
= (bEne
) ? md
->ebin
->nener
: 0;
1274 fr
.ener
= md
->ebin
->e
;
1275 ndisre
= bDR
? fcd
->disres
.npair
: 0;
1276 disre_rm3tav
= fcd
->disres
.rm3tav
;
1277 disre_rt
= fcd
->disres
.rt
;
1278 /* Optional additional old-style (real-only) blocks. */
1279 for (i
= 0; i
< enxNR
; i
++)
1283 if (fcd
->orires
.nr
> 0 && bOR
)
1285 diagonalize_orires_tensors(&(fcd
->orires
));
1286 nr
[enxOR
] = fcd
->orires
.nr
;
1287 block
[enxOR
] = fcd
->orires
.otav
;
1289 nr
[enxORI
] = (fcd
->orires
.oinsl
!= fcd
->orires
.otav
) ?
1291 block
[enxORI
] = fcd
->orires
.oinsl
;
1292 id
[enxORI
] = enxORI
;
1293 nr
[enxORT
] = fcd
->orires
.nex
*12;
1294 block
[enxORT
] = fcd
->orires
.eig
;
1295 id
[enxORT
] = enxORT
;
1298 /* whether we are going to wrte anything out: */
1299 if (fr
.nre
|| ndisre
|| nr
[enxOR
] || nr
[enxORI
])
1302 /* the old-style blocks go first */
1304 for (i
= 0; i
< enxNR
; i
++)
1311 add_blocks_enxframe(&fr
, fr
.nblock
);
1312 for (b
= 0; b
< fr
.nblock
; b
++)
1314 add_subblocks_enxblock(&(fr
.block
[b
]), 1);
1315 fr
.block
[b
].id
= id
[b
];
1316 fr
.block
[b
].sub
[0].nr
= nr
[b
];
1318 fr
.block
[b
].sub
[0].type
= xdr_datatype_float
;
1319 fr
.block
[b
].sub
[0].fval
= block
[b
];
1321 fr
.block
[b
].sub
[0].type
= xdr_datatype_double
;
1322 fr
.block
[b
].sub
[0].dval
= block
[b
];
1326 /* check for disre block & fill it. */
1331 add_blocks_enxframe(&fr
, fr
.nblock
);
1333 add_subblocks_enxblock(&(fr
.block
[db
]), 2);
1334 fr
.block
[db
].id
= enxDISRE
;
1335 fr
.block
[db
].sub
[0].nr
= ndisre
;
1336 fr
.block
[db
].sub
[1].nr
= ndisre
;
1338 fr
.block
[db
].sub
[0].type
= xdr_datatype_float
;
1339 fr
.block
[db
].sub
[1].type
= xdr_datatype_float
;
1340 fr
.block
[db
].sub
[0].fval
= disre_rt
;
1341 fr
.block
[db
].sub
[1].fval
= disre_rm3tav
;
1343 fr
.block
[db
].sub
[0].type
= xdr_datatype_double
;
1344 fr
.block
[db
].sub
[1].type
= xdr_datatype_double
;
1345 fr
.block
[db
].sub
[0].dval
= disre_rt
;
1346 fr
.block
[db
].sub
[1].dval
= disre_rm3tav
;
1349 /* here we can put new-style blocks */
1351 /* Free energy perturbation blocks */
1354 mde_delta_h_coll_handle_block(md
->dhc
, &fr
, fr
.nblock
);
1357 /* we can now free & reset the data in the blocks */
1360 mde_delta_h_coll_reset(md
->dhc
);
1363 /* do the actual I/O */
1364 do_enx(fp_ene
, &fr
);
1367 /* We have stored the sums, so reset the sum history */
1368 reset_ebin_sums(md
->ebin
);
1376 pprint(log
, "A V E R A G E S", md
);
1382 pprint(log
, "R M S - F L U C T U A T I O N S", md
);
1386 gmx_fatal(FARGS
, "Invalid print mode (%d)", mode
);
1391 for (i
= 0; i
< opts
->ngtc
; i
++)
1393 if (opts
->annealing
[i
] != eannNO
)
1395 fprintf(log
, "Current ref_t for group %s: %8.1f\n",
1396 *(groups
->grpname
[groups
->grps
[egcTC
].nm_ind
[i
]]),
1400 if (mode
== eprNORMAL
&& fcd
->orires
.nr
> 0)
1402 print_orires_log(log
, &(fcd
->orires
));
1404 fprintf(log
, " Energies (%s)\n", unit_energy
);
1405 pr_ebin(log
, md
->ebin
, md
->ie
, md
->f_nre
+md
->nCrmsd
, 5, mode
, TRUE
);
1408 if (mode
== eprAVER
)
1412 pr_ebin(log
, md
->ebin
, md
->ib
, md
->bTricl
? NTRICLBOXS
: NBOXS
, 5,
1418 fprintf(log
, " Constraint Virial (%s)\n", unit_energy
);
1419 pr_ebin(log
, md
->ebin
, md
->isvir
, 9, 3, mode
, FALSE
);
1421 fprintf(log
, " Force Virial (%s)\n", unit_energy
);
1422 pr_ebin(log
, md
->ebin
, md
->ifvir
, 9, 3, mode
, FALSE
);
1425 fprintf(log
, " Total Virial (%s)\n", unit_energy
);
1426 pr_ebin(log
, md
->ebin
, md
->ivir
, 9, 3, mode
, FALSE
);
1428 fprintf(log
, " Pressure (%s)\n", unit_pres_bar
);
1429 pr_ebin(log
, md
->ebin
, md
->ipres
, 9, 3, mode
, FALSE
);
1433 fprintf(log
, " Total Dipole (%s)\n", unit_dipole_D
);
1434 pr_ebin(log
, md
->ebin
, md
->imu
, 3, 3, mode
, FALSE
);
1440 if (md
->print_grpnms
== nullptr)
1442 snew(md
->print_grpnms
, md
->nE
);
1444 for (i
= 0; (i
< md
->nEg
); i
++)
1446 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
1447 for (j
= i
; (j
< md
->nEg
); j
++)
1449 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
1450 sprintf(buf
, "%s-%s", *(groups
->grpname
[ni
]),
1451 *(groups
->grpname
[nj
]));
1452 md
->print_grpnms
[n
++] = gmx_strdup(buf
);
1456 sprintf(buf
, "Epot (%s)", unit_energy
);
1457 fprintf(log
, "%15s ", buf
);
1458 for (i
= 0; (i
< egNR
); i
++)
1462 fprintf(log
, "%12s ", egrp_nm
[i
]);
1466 for (i
= 0; (i
< md
->nE
); i
++)
1468 fprintf(log
, "%15s", md
->print_grpnms
[i
]);
1469 pr_ebin(log
, md
->ebin
, md
->igrp
[i
], md
->nEc
, md
->nEc
, mode
,
1476 pr_ebin(log
, md
->ebin
, md
->itemp
, md
->nTC
, 4, mode
, TRUE
);
1481 fprintf(log
, "%15s %12s %12s %12s\n",
1482 "Group", "Ux", "Uy", "Uz");
1483 for (i
= 0; (i
< md
->nU
); i
++)
1485 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
1486 fprintf(log
, "%15s", *groups
->grpname
[ni
]);
1487 pr_ebin(log
, md
->ebin
, md
->iu
+3*i
, 3, 3, mode
, FALSE
);
1496 void update_energyhistory(energyhistory_t
* enerhist
, t_mdebin
* mdebin
)
1498 const t_ebin
* const ebin
= mdebin
->ebin
;
1500 enerhist
->nsteps
= ebin
->nsteps
;
1501 enerhist
->nsum
= ebin
->nsum
;
1502 enerhist
->nsteps_sim
= ebin
->nsteps_sim
;
1503 enerhist
->nsum_sim
= ebin
->nsum_sim
;
1507 /* This will only actually resize the first time */
1508 enerhist
->ener_ave
.resize(ebin
->nener
);
1509 enerhist
->ener_sum
.resize(ebin
->nener
);
1511 for (int i
= 0; i
< ebin
->nener
; i
++)
1513 enerhist
->ener_ave
[i
] = ebin
->e
[i
].eav
;
1514 enerhist
->ener_sum
[i
] = ebin
->e
[i
].esum
;
1518 if (ebin
->nsum_sim
> 0)
1520 /* This will only actually resize the first time */
1521 enerhist
->ener_sum_sim
.resize(ebin
->nener
);
1523 for (int i
= 0; i
< ebin
->nener
; i
++)
1525 enerhist
->ener_sum_sim
[i
] = ebin
->e_sim
[i
].esum
;
1530 mde_delta_h_coll_update_energyhistory(mdebin
->dhc
, enerhist
);
1534 void restore_energyhistory_from_state(t_mdebin
* mdebin
,
1535 energyhistory_t
* enerhist
)
1537 unsigned int nener
= static_cast<unsigned int>(mdebin
->ebin
->nener
);
1539 if ((enerhist
->nsum
> 0 && nener
!= enerhist
->ener_sum
.size()) ||
1540 (enerhist
->nsum_sim
> 0 && nener
!= enerhist
->ener_sum_sim
.size()))
1542 gmx_fatal(FARGS
, "Mismatch between number of energies in run input (%d) and checkpoint file (%u or %u).",
1543 nener
, enerhist
->ener_sum
.size(), enerhist
->ener_sum_sim
.size());
1546 mdebin
->ebin
->nsteps
= enerhist
->nsteps
;
1547 mdebin
->ebin
->nsum
= enerhist
->nsum
;
1548 mdebin
->ebin
->nsteps_sim
= enerhist
->nsteps_sim
;
1549 mdebin
->ebin
->nsum_sim
= enerhist
->nsum_sim
;
1551 for (int i
= 0; i
< mdebin
->ebin
->nener
; i
++)
1553 mdebin
->ebin
->e
[i
].eav
=
1554 (enerhist
->nsum
> 0 ? enerhist
->ener_ave
[i
] : 0);
1555 mdebin
->ebin
->e
[i
].esum
=
1556 (enerhist
->nsum
> 0 ? enerhist
->ener_sum
[i
] : 0);
1557 mdebin
->ebin
->e_sim
[i
].esum
=
1558 (enerhist
->nsum_sim
> 0 ? enerhist
->ener_sum_sim
[i
] : 0);
1562 mde_delta_h_coll_restore_energyhistory(mdebin
->dhc
, enerhist
->deltaHForeignLambdas
.get());