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, 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.
39 #include "gromacs/legacyheaders/mdebin.h"
45 #include "gromacs/fileio/enxio.h"
46 #include "gromacs/fileio/gmxfio.h"
47 #include "gromacs/fileio/xvgr.h"
48 #include "gromacs/legacyheaders/disre.h"
49 #include "gromacs/legacyheaders/mdrun.h"
50 #include "gromacs/legacyheaders/names.h"
51 #include "gromacs/legacyheaders/network.h"
52 #include "gromacs/legacyheaders/orires.h"
53 #include "gromacs/legacyheaders/typedefs.h"
54 #include "gromacs/legacyheaders/types/group.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/pbcutil/pbc.h"
60 #include "gromacs/pulling/pull.h"
61 #include "gromacs/topology/mtop_util.h"
62 #include "gromacs/utility/arraysize.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/smalloc.h"
66 static const char *conrmsd_nm
[] = { "Constr. rmsd", "Constr.2 rmsd" };
68 static const char *boxs_nm
[] = { "Box-X", "Box-Y", "Box-Z" };
70 static const char *tricl_boxs_nm
[] = {
71 "Box-XX", "Box-YY", "Box-ZZ",
72 "Box-YX", "Box-ZX", "Box-ZY"
75 static const char *vol_nm
[] = { "Volume" };
77 static const char *dens_nm
[] = {"Density" };
79 static const char *pv_nm
[] = {"pV" };
81 static const char *enthalpy_nm
[] = {"Enthalpy" };
83 static const char *boxvel_nm
[] = {
84 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
85 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
88 #define NBOXS asize(boxs_nm)
89 #define NTRICLBOXS asize(tricl_boxs_nm)
91 t_mdebin
*init_mdebin(ener_file_t fp_ene
,
92 const gmx_mtop_t
*mtop
,
96 const char *ener_nm
[F_NRE
];
97 static const char *vir_nm
[] = {
98 "Vir-XX", "Vir-XY", "Vir-XZ",
99 "Vir-YX", "Vir-YY", "Vir-YZ",
100 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
102 static const char *sv_nm
[] = {
103 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
104 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
105 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
107 static const char *fv_nm
[] = {
108 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
109 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
110 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
112 static const char *pres_nm
[] = {
113 "Pres-XX", "Pres-XY", "Pres-XZ",
114 "Pres-YX", "Pres-YY", "Pres-YZ",
115 "Pres-ZX", "Pres-ZY", "Pres-ZZ"
117 static const char *surft_nm
[] = {
120 static const char *mu_nm
[] = {
121 "Mu-X", "Mu-Y", "Mu-Z"
123 static const char *vcos_nm
[] = {
126 static const char *visc_nm
[] = {
129 static const char *baro_nm
[] = {
134 const gmx_groups_t
*groups
;
139 int i
, j
, ni
, nj
, n
, k
, kk
, ncon
, nset
;
149 if (EI_DYNAMICS(ir
->eI
))
151 md
->delta_t
= ir
->delta_t
;
158 groups
= &mtop
->groups
;
160 bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
161 b14
= (gmx_mtop_ftype_count(mtop
, F_LJ14
) > 0 ||
162 gmx_mtop_ftype_count(mtop
, F_LJC14_Q
) > 0);
164 ncon
= gmx_mtop_ftype_count(mtop
, F_CONSTR
);
165 nset
= gmx_mtop_ftype_count(mtop
, F_SETTLE
);
166 md
->bConstr
= (ncon
> 0 || nset
> 0);
167 md
->bConstrVir
= FALSE
;
170 if (ncon
> 0 && ir
->eConstrAlg
== econtLINCS
)
181 md
->bConstrVir
= (getenv("GMX_CONSTRAINTVIR") != NULL
);
188 /* Energy monitoring */
189 for (i
= 0; i
< egNR
; i
++)
191 md
->bEInd
[i
] = FALSE
;
194 for (i
= 0; i
< F_NRE
; i
++)
196 md
->bEner
[i
] = FALSE
;
199 md
->bEner
[i
] = !bBHAM
;
201 else if (i
== F_BHAM
)
203 md
->bEner
[i
] = bBHAM
;
207 md
->bEner
[i
] = ir
->bQMMM
;
209 else if (i
== F_COUL_LR
)
211 md
->bEner
[i
] = (ir
->rcoulomb
> ir
->rlist
);
213 else if (i
== F_LJ_LR
)
215 md
->bEner
[i
] = (!bBHAM
&& ir
->rvdw
> ir
->rlist
);
217 else if (i
== F_BHAM_LR
)
219 md
->bEner
[i
] = (bBHAM
&& ir
->rvdw
> ir
->rlist
);
221 else if (i
== F_RF_EXCL
)
223 md
->bEner
[i
] = (EEL_RF(ir
->coulombtype
) && ir
->coulombtype
!= eelRF_NEC
&& ir
->cutoff_scheme
== ecutsGROUP
);
225 else if (i
== F_COUL_RECIP
)
227 md
->bEner
[i
] = EEL_FULL(ir
->coulombtype
);
229 else if (i
== F_LJ_RECIP
)
231 md
->bEner
[i
] = EVDW_PME(ir
->vdwtype
);
233 else if (i
== F_LJ14
)
237 else if (i
== F_COUL14
)
241 else if (i
== F_LJC14_Q
|| i
== F_LJC_PAIRS_NB
)
243 md
->bEner
[i
] = FALSE
;
245 else if ((i
== F_DVDL_COUL
&& ir
->fepvals
->separate_dvdl
[efptCOUL
]) ||
246 (i
== F_DVDL_VDW
&& ir
->fepvals
->separate_dvdl
[efptVDW
]) ||
247 (i
== F_DVDL_BONDED
&& ir
->fepvals
->separate_dvdl
[efptBONDED
]) ||
248 (i
== F_DVDL_RESTRAINT
&& ir
->fepvals
->separate_dvdl
[efptRESTRAINT
]) ||
249 (i
== F_DKDL
&& ir
->fepvals
->separate_dvdl
[efptMASS
]) ||
250 (i
== F_DVDL
&& ir
->fepvals
->separate_dvdl
[efptFEP
]))
252 md
->bEner
[i
] = (ir
->efep
!= efepNO
);
254 else if ((interaction_function
[i
].flags
& IF_VSITE
) ||
255 (i
== F_CONSTR
) || (i
== F_CONSTRNC
) || (i
== F_SETTLE
))
257 md
->bEner
[i
] = FALSE
;
259 else if ((i
== F_COUL_SR
) || (i
== F_EPOT
) || (i
== F_PRES
) || (i
== F_EQM
))
263 else if ((i
== F_GBPOL
) && ir
->implicit_solvent
== eisGBSA
)
267 else if ((i
== F_NPSOLVATION
) && ir
->implicit_solvent
== eisGBSA
&& (ir
->sa_algorithm
!= esaNO
))
271 else if ((i
== F_GB12
) || (i
== F_GB13
) || (i
== F_GB14
))
273 md
->bEner
[i
] = FALSE
;
275 else if ((i
== F_ETOT
) || (i
== F_EKIN
) || (i
== F_TEMP
))
277 md
->bEner
[i
] = EI_DYNAMICS(ir
->eI
);
279 else if (i
== F_DISPCORR
|| i
== F_PDISPCORR
)
281 md
->bEner
[i
] = (ir
->eDispCorr
!= edispcNO
);
283 else if (i
== F_DISRESVIOL
)
285 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_DISRES
) > 0);
287 else if (i
== F_ORIRESDEV
)
289 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_ORIRES
) > 0);
291 else if (i
== F_CONNBONDS
)
293 md
->bEner
[i
] = FALSE
;
295 else if (i
== F_COM_PULL
)
297 md
->bEner
[i
] = (ir
->bPull
&& pull_have_potential(ir
->pull_work
));
299 else if (i
== F_ECONSERVED
)
301 md
->bEner
[i
] = ((ir
->etc
== etcNOSEHOOVER
|| ir
->etc
== etcVRESCALE
) &&
302 (ir
->epc
== epcNO
|| ir
->epc
== epcMTTK
));
306 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, i
) > 0);
310 /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
311 if (ir
->bAdress
&& !debug
)
313 for (i
= 0; i
< F_NRE
; i
++)
315 md
->bEner
[i
] = FALSE
;
332 for (i
= 0; i
< F_NRE
; i
++)
336 ener_nm
[md
->f_nre
] = interaction_function
[i
].longname
;
342 md
->bDiagPres
= !TRICLINIC(ir
->ref_p
);
343 md
->ref_p
= (ir
->ref_p
[XX
][XX
]+ir
->ref_p
[YY
][YY
]+ir
->ref_p
[ZZ
][ZZ
])/DIM
;
344 md
->bTricl
= TRICLINIC(ir
->compress
) || TRICLINIC(ir
->deform
);
345 md
->bDynBox
= DYNAMIC_BOX(*ir
);
347 md
->bNHC_trotter
= IR_NVT_TROTTER(ir
);
348 md
->bPrintNHChains
= ir
->bPrintNHChains
;
349 md
->bMTTK
= (IR_NPT_TROTTER(ir
) || IR_NPH_TROTTER(ir
));
350 md
->bMu
= NEED_MUTOT(*ir
);
352 md
->ebin
= mk_ebin();
353 /* Pass NULL for unit to let get_ebin_space determine the units
354 * for interaction_function[i].longname
356 md
->ie
= get_ebin_space(md
->ebin
, md
->f_nre
, ener_nm
, NULL
);
359 /* This should be called directly after the call for md->ie,
360 * such that md->iconrmsd follows directly in the list.
362 md
->iconrmsd
= get_ebin_space(md
->ebin
, md
->nCrmsd
, conrmsd_nm
, "");
366 md
->ib
= get_ebin_space(md
->ebin
,
367 md
->bTricl
? NTRICLBOXS
: NBOXS
,
368 md
->bTricl
? tricl_boxs_nm
: boxs_nm
,
370 md
->ivol
= get_ebin_space(md
->ebin
, 1, vol_nm
, unit_volume
);
371 md
->idens
= get_ebin_space(md
->ebin
, 1, dens_nm
, unit_density_SI
);
374 md
->ipv
= get_ebin_space(md
->ebin
, 1, pv_nm
, unit_energy
);
375 md
->ienthalpy
= get_ebin_space(md
->ebin
, 1, enthalpy_nm
, unit_energy
);
380 md
->isvir
= get_ebin_space(md
->ebin
, asize(sv_nm
), sv_nm
, unit_energy
);
381 md
->ifvir
= get_ebin_space(md
->ebin
, asize(fv_nm
), fv_nm
, unit_energy
);
385 md
->ivir
= get_ebin_space(md
->ebin
, asize(vir_nm
), vir_nm
, unit_energy
);
389 md
->ipres
= get_ebin_space(md
->ebin
, asize(pres_nm
), pres_nm
, unit_pres_bar
);
393 md
->isurft
= get_ebin_space(md
->ebin
, asize(surft_nm
), surft_nm
,
396 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
398 md
->ipc
= get_ebin_space(md
->ebin
, md
->bTricl
? 6 : 3,
399 boxvel_nm
, unit_vel
);
403 md
->imu
= get_ebin_space(md
->ebin
, asize(mu_nm
), mu_nm
, unit_dipole_D
);
405 if (ir
->cos_accel
!= 0)
407 md
->ivcos
= get_ebin_space(md
->ebin
, asize(vcos_nm
), vcos_nm
, unit_vel
);
408 md
->ivisc
= get_ebin_space(md
->ebin
, asize(visc_nm
), visc_nm
,
412 /* Energy monitoring */
413 for (i
= 0; i
< egNR
; i
++)
415 md
->bEInd
[i
] = FALSE
;
417 md
->bEInd
[egCOULSR
] = TRUE
;
418 md
->bEInd
[egLJSR
] = TRUE
;
420 if (ir
->rcoulomb
> ir
->rlist
)
422 md
->bEInd
[egCOULLR
] = TRUE
;
426 if (ir
->rvdw
> ir
->rlist
)
428 md
->bEInd
[egLJLR
] = TRUE
;
433 md
->bEInd
[egLJSR
] = FALSE
;
434 md
->bEInd
[egBHAMSR
] = TRUE
;
435 if (ir
->rvdw
> ir
->rlist
)
437 md
->bEInd
[egBHAMLR
] = TRUE
;
442 md
->bEInd
[egLJ14
] = TRUE
;
443 md
->bEInd
[egCOUL14
] = TRUE
;
446 for (i
= 0; (i
< egNR
); i
++)
454 n
= groups
->grps
[egcENER
].nr
;
455 /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
458 /*standard simulation*/
460 md
->nE
= (n
*(n
+1))/2;
464 /*AdResS simulation*/
471 snew(md
->igrp
, md
->nE
);
476 for (k
= 0; (k
< md
->nEc
); k
++)
478 snew(gnm
[k
], STRLEN
);
480 for (i
= 0; (i
< groups
->grps
[egcENER
].nr
); i
++)
482 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
483 for (j
= i
; (j
< groups
->grps
[egcENER
].nr
); j
++)
485 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
486 for (k
= kk
= 0; (k
< egNR
); k
++)
490 sprintf(gnm
[kk
], "%s:%s-%s", egrp_nm
[k
],
491 *(groups
->grpname
[ni
]), *(groups
->grpname
[nj
]));
495 md
->igrp
[n
] = get_ebin_space(md
->ebin
, md
->nEc
,
496 (const char **)gnm
, unit_energy
);
500 for (k
= 0; (k
< md
->nEc
); k
++)
508 gmx_incons("Number of energy terms wrong");
512 md
->nTC
= groups
->grps
[egcTC
].nr
;
513 md
->nNHC
= ir
->opts
.nhchainlength
; /* shorthand for number of NH chains */
516 md
->nTCP
= 1; /* assume only one possible coupling system for barostat
523 if (md
->etc
== etcNOSEHOOVER
)
525 if (md
->bNHC_trotter
)
527 md
->mde_n
= 2*md
->nNHC
*md
->nTC
;
531 md
->mde_n
= 2*md
->nTC
;
533 if (md
->epc
== epcMTTK
)
535 md
->mdeb_n
= 2*md
->nNHC
*md
->nTCP
;
544 snew(md
->tmp_r
, md
->mde_n
);
545 snew(md
->tmp_v
, md
->mde_n
);
546 snew(md
->grpnms
, md
->mde_n
);
549 for (i
= 0; (i
< md
->nTC
); i
++)
551 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
552 sprintf(buf
, "T-%s", *(groups
->grpname
[ni
]));
553 grpnms
[i
] = gmx_strdup(buf
);
555 md
->itemp
= get_ebin_space(md
->ebin
, md
->nTC
, (const char **)grpnms
,
558 if (md
->etc
== etcNOSEHOOVER
)
560 if (md
->bPrintNHChains
)
562 if (md
->bNHC_trotter
)
564 for (i
= 0; (i
< md
->nTC
); i
++)
566 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
567 bufi
= *(groups
->grpname
[ni
]);
568 for (j
= 0; (j
< md
->nNHC
); j
++)
570 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
571 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
572 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
573 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
576 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
577 (const char **)grpnms
, unit_invtime
);
580 for (i
= 0; (i
< md
->nTCP
); i
++)
582 bufi
= baro_nm
[0]; /* All barostat DOF's together for now. */
583 for (j
= 0; (j
< md
->nNHC
); j
++)
585 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
586 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
587 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
588 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
591 md
->itcb
= get_ebin_space(md
->ebin
, md
->mdeb_n
,
592 (const char **)grpnms
, unit_invtime
);
597 for (i
= 0; (i
< md
->nTC
); i
++)
599 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
600 bufi
= *(groups
->grpname
[ni
]);
601 sprintf(buf
, "Xi-%s", bufi
);
602 grpnms
[2*i
] = gmx_strdup(buf
);
603 sprintf(buf
, "vXi-%s", bufi
);
604 grpnms
[2*i
+1] = gmx_strdup(buf
);
606 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
607 (const char **)grpnms
, unit_invtime
);
611 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
612 md
->etc
== etcVRESCALE
)
614 for (i
= 0; (i
< md
->nTC
); i
++)
616 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
617 sprintf(buf
, "Lamb-%s", *(groups
->grpname
[ni
]));
618 grpnms
[i
] = gmx_strdup(buf
);
620 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
, (const char **)grpnms
, "");
626 md
->nU
= groups
->grps
[egcACC
].nr
;
629 snew(grpnms
, 3*md
->nU
);
630 for (i
= 0; (i
< md
->nU
); i
++)
632 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
633 sprintf(buf
, "Ux-%s", *(groups
->grpname
[ni
]));
634 grpnms
[3*i
+XX
] = gmx_strdup(buf
);
635 sprintf(buf
, "Uy-%s", *(groups
->grpname
[ni
]));
636 grpnms
[3*i
+YY
] = gmx_strdup(buf
);
637 sprintf(buf
, "Uz-%s", *(groups
->grpname
[ni
]));
638 grpnms
[3*i
+ZZ
] = gmx_strdup(buf
);
640 md
->iu
= get_ebin_space(md
->ebin
, 3*md
->nU
, (const char **)grpnms
, unit_vel
);
646 do_enxnms(fp_ene
, &md
->ebin
->nener
, &md
->ebin
->enm
);
649 md
->print_grpnms
= NULL
;
651 /* check whether we're going to write dh histograms */
653 if (ir
->fepvals
->separate_dhdl_file
== esepdhdlfileNO
)
655 /* Currently dh histograms are only written with dynamics */
656 if (EI_DYNAMICS(ir
->eI
))
660 mde_delta_h_coll_init(md
->dhc
, ir
);
663 snew(md
->dE
, ir
->fepvals
->n_lambda
);
667 md
->fp_dhdl
= fp_dhdl
;
668 snew(md
->dE
, ir
->fepvals
->n_lambda
);
673 snew(md
->temperatures
, ir
->fepvals
->n_lambda
);
674 for (i
= 0; i
< ir
->fepvals
->n_lambda
; i
++)
676 md
->temperatures
[i
] = ir
->simtempvals
->temperatures
[i
];
682 /* print a lambda vector to a string
683 fep = the inputrec's FEP input data
684 i = the index of the lambda vector
685 get_native_lambda = whether to print the native lambda
686 get_names = whether to print the names rather than the values
687 str = the pre-allocated string buffer to print to. */
688 static void print_lambda_vector(t_lambda
*fep
, int i
,
689 gmx_bool get_native_lambda
, gmx_bool get_names
,
695 for (j
= 0; j
< efptNR
; j
++)
697 if (fep
->separate_dvdl
[j
])
702 str
[0] = 0; /* reset the string */
705 str
+= sprintf(str
, "("); /* set the opening parenthesis*/
707 for (j
= 0; j
< efptNR
; j
++)
709 if (fep
->separate_dvdl
[j
])
713 if (get_native_lambda
&& fep
->init_lambda
>= 0)
715 str
+= sprintf(str
, "%.4f", fep
->init_lambda
);
719 str
+= sprintf(str
, "%.4f", fep
->all_lambda
[j
][i
]);
724 str
+= sprintf(str
, "%s", efpt_singular_names
[j
]);
726 /* print comma for the next item */
729 str
+= sprintf(str
, ", ");
736 /* and add the closing parenthesis */
742 extern FILE *open_dhdl(const char *filename
, const t_inputrec
*ir
,
743 const gmx_output_env_t
*oenv
)
746 const char *dhdl
= "dH/d\\lambda", *deltag
= "\\DeltaH", *lambda
= "\\lambda",
747 *lambdastate
= "\\lambda state";
748 char title
[STRLEN
], label_x
[STRLEN
], label_y
[STRLEN
];
749 int i
, nps
, nsets
, nsets_de
, nsetsbegin
;
750 int n_lambda_terms
= 0;
751 t_lambda
*fep
= ir
->fepvals
; /* for simplicity */
752 t_expanded
*expand
= ir
->expandedvals
;
754 char buf
[STRLEN
], lambda_vec_str
[STRLEN
], lambda_name_str
[STRLEN
];
760 gmx_bool write_pV
= FALSE
;
762 /* count the number of different lambda terms */
763 for (i
= 0; i
< efptNR
; i
++)
765 if (fep
->separate_dvdl
[i
])
771 if (fep
->n_lambda
== 0)
773 sprintf(title
, "%s", dhdl
);
774 sprintf(label_x
, "Time (ps)");
775 sprintf(label_y
, "%s (%s %s)",
776 dhdl
, unit_energy
, "[\\lambda]\\S-1\\N");
780 sprintf(title
, "%s and %s", dhdl
, deltag
);
781 sprintf(label_x
, "Time (ps)");
782 sprintf(label_y
, "%s and %s (%s %s)",
783 dhdl
, deltag
, unit_energy
, "[\\8l\\4]\\S-1\\N");
785 fp
= gmx_fio_fopen(filename
, "w+");
786 xvgr_header(fp
, title
, label_x
, label_y
, exvggtXNY
, oenv
);
790 bufplace
= sprintf(buf
, "T = %g (K) ",
793 if ((ir
->efep
!= efepSLOWGROWTH
) && (ir
->efep
!= efepEXPANDED
))
795 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
797 /* compatibility output */
798 sprintf(&(buf
[bufplace
]), "%s = %.4f", lambda
, fep
->init_lambda
);
802 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, FALSE
,
804 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, TRUE
,
806 sprintf(&(buf
[bufplace
]), "%s %d: %s = %s",
807 lambdastate
, fep
->init_fep_state
,
808 lambda_name_str
, lambda_vec_str
);
811 xvgr_subtitle(fp
, buf
, oenv
);
815 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
817 nsets_dhdl
= n_lambda_terms
;
819 /* count the number of delta_g states */
820 nsets_de
= fep
->lambda_stop_n
- fep
->lambda_start_n
;
822 nsets
= nsets_dhdl
+ nsets_de
; /* dhdl + fep differences */
824 if (fep
->n_lambda
> 0 && (expand
->elmcmove
> elmcmoveNO
))
826 nsets
+= 1; /*add fep state for expanded ensemble */
829 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
831 nsets
+= 1; /* add energy to the dhdl as well */
835 if ((ir
->epc
!= epcNO
) && (fep
->n_lambda
> 0) && (fep
->init_lambda
< 0))
837 nsetsextend
+= 1; /* for PV term, other terms possible if required for
838 the reduced potential (only needed with foreign
839 lambda, and only output when init_lambda is not
840 set in order to maintain compatibility of the
844 snew(setname
, nsetsextend
);
846 if (expand
->elmcmove
> elmcmoveNO
)
848 /* state for the fep_vals, if we have alchemical sampling */
849 sprintf(buf
, "%s", "Thermodynamic state");
850 setname
[s
] = gmx_strdup(buf
);
854 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
856 switch (fep
->edHdLPrintEnergy
)
858 case edHdLPrintEnergyPOTENTIAL
:
859 sprintf(buf
, "%s (%s)", "Potential Energy", unit_energy
);
861 case edHdLPrintEnergyTOTAL
:
862 case edHdLPrintEnergyYES
:
864 sprintf(buf
, "%s (%s)", "Total Energy", unit_energy
);
866 setname
[s
] = gmx_strdup(buf
);
870 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
872 for (i
= 0; i
< efptNR
; i
++)
874 if (fep
->separate_dvdl
[i
])
877 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
879 /* compatibility output */
880 sprintf(buf
, "%s %s %.4f", dhdl
, lambda
, fep
->init_lambda
);
884 double lam
= fep
->init_lambda
;
885 if (fep
->init_lambda
< 0)
887 lam
= fep
->all_lambda
[i
][fep
->init_fep_state
];
889 sprintf(buf
, "%s %s = %.4f", dhdl
, efpt_singular_names
[i
],
892 setname
[s
] = gmx_strdup(buf
);
898 if (fep
->n_lambda
> 0)
900 /* g_bar has to determine the lambda values used in this simulation
901 * from this xvg legend.
904 if (expand
->elmcmove
> elmcmoveNO
)
906 nsetsbegin
= 1; /* for including the expanded ensemble */
913 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
917 nsetsbegin
+= nsets_dhdl
;
919 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
921 print_lambda_vector(fep
, i
, FALSE
, FALSE
, lambda_vec_str
);
922 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
924 /* for compatible dhdl.xvg files */
925 nps
= sprintf(buf
, "%s %s %s", deltag
, lambda
, lambda_vec_str
);
929 nps
= sprintf(buf
, "%s %s to %s", deltag
, lambda
, lambda_vec_str
);
934 /* print the temperature for this state if doing simulated annealing */
935 sprintf(&buf
[nps
], "T = %g (%s)",
936 ir
->simtempvals
->temperatures
[s
-(nsetsbegin
)],
939 setname
[s
] = gmx_strdup(buf
);
944 sprintf(buf
, "pV (%s)", unit_energy
);
945 setname
[nsetsextend
-1] = gmx_strdup(buf
); /* the first entry after
949 xvgr_legend(fp
, nsetsextend
, (const char **)setname
, oenv
);
951 for (s
= 0; s
< nsetsextend
; s
++)
961 static void copy_energy(t_mdebin
*md
, real e
[], real ecpy
[])
965 for (i
= j
= 0; (i
< F_NRE
); i
++)
974 gmx_incons("Number of energy terms wrong");
978 void upd_mdebin(t_mdebin
*md
,
983 gmx_enerdata_t
*enerd
,
992 gmx_ekindata_t
*ekind
,
996 int i
, j
, k
, kk
, n
, gid
;
997 real crmsd
[2], tmp6
[6];
998 real bs
[NTRICLBOXS
], vol
, dens
, pv
, enthalpy
;
1001 double store_dhdl
[efptNR
];
1002 real store_energy
= 0;
1005 /* Do NOT use the box in the state variable, but the separate box provided
1006 * as an argument. This is because we sometimes need to write the box from
1007 * the last timestep to match the trajectory frames.
1009 copy_energy(md
, enerd
->term
, ecopy
);
1010 add_ebin(md
->ebin
, md
->ie
, md
->f_nre
, ecopy
, bSum
);
1013 crmsd
[0] = constr_rmsd(constr
, FALSE
);
1016 crmsd
[1] = constr_rmsd(constr
, TRUE
);
1018 add_ebin(md
->ebin
, md
->iconrmsd
, md
->nCrmsd
, crmsd
, FALSE
);
1025 bs
[0] = box
[XX
][XX
];
1026 bs
[1] = box
[YY
][YY
];
1027 bs
[2] = box
[ZZ
][ZZ
];
1028 bs
[3] = box
[YY
][XX
];
1029 bs
[4] = box
[ZZ
][XX
];
1030 bs
[5] = box
[ZZ
][YY
];
1035 bs
[0] = box
[XX
][XX
];
1036 bs
[1] = box
[YY
][YY
];
1037 bs
[2] = box
[ZZ
][ZZ
];
1040 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
1041 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
1042 add_ebin(md
->ebin
, md
->ib
, nboxs
, bs
, bSum
);
1043 add_ebin(md
->ebin
, md
->ivol
, 1, &vol
, bSum
);
1044 add_ebin(md
->ebin
, md
->idens
, 1, &dens
, bSum
);
1048 /* This is pV (in kJ/mol). The pressure is the reference pressure,
1049 not the instantaneous pressure */
1050 pv
= vol
*md
->ref_p
/PRESFAC
;
1052 add_ebin(md
->ebin
, md
->ipv
, 1, &pv
, bSum
);
1053 enthalpy
= pv
+ enerd
->term
[F_ETOT
];
1054 add_ebin(md
->ebin
, md
->ienthalpy
, 1, &enthalpy
, bSum
);
1059 add_ebin(md
->ebin
, md
->isvir
, 9, svir
[0], bSum
);
1060 add_ebin(md
->ebin
, md
->ifvir
, 9, fvir
[0], bSum
);
1064 add_ebin(md
->ebin
, md
->ivir
, 9, vir
[0], bSum
);
1068 add_ebin(md
->ebin
, md
->ipres
, 9, pres
[0], bSum
);
1072 tmp
= (pres
[ZZ
][ZZ
]-(pres
[XX
][XX
]+pres
[YY
][YY
])*0.5)*box
[ZZ
][ZZ
];
1073 add_ebin(md
->ebin
, md
->isurft
, 1, &tmp
, bSum
);
1075 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
1077 tmp6
[0] = state
->boxv
[XX
][XX
];
1078 tmp6
[1] = state
->boxv
[YY
][YY
];
1079 tmp6
[2] = state
->boxv
[ZZ
][ZZ
];
1080 tmp6
[3] = state
->boxv
[YY
][XX
];
1081 tmp6
[4] = state
->boxv
[ZZ
][XX
];
1082 tmp6
[5] = state
->boxv
[ZZ
][YY
];
1083 add_ebin(md
->ebin
, md
->ipc
, md
->bTricl
? 6 : 3, tmp6
, bSum
);
1087 add_ebin(md
->ebin
, md
->imu
, 3, mu_tot
, bSum
);
1089 if (ekind
&& ekind
->cosacc
.cos_accel
!= 0)
1091 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
1092 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
1093 add_ebin(md
->ebin
, md
->ivcos
, 1, &(ekind
->cosacc
.vcos
), bSum
);
1094 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
1095 tmp
= 1/(ekind
->cosacc
.cos_accel
/(ekind
->cosacc
.vcos
*PICO
)
1096 *dens
*sqr(box
[ZZ
][ZZ
]*NANO
/(2*M_PI
)));
1097 add_ebin(md
->ebin
, md
->ivisc
, 1, &tmp
, bSum
);
1102 for (i
= 0; (i
< md
->nEg
); i
++)
1104 for (j
= i
; (j
< md
->nEg
); j
++)
1106 gid
= GID(i
, j
, md
->nEg
);
1107 for (k
= kk
= 0; (k
< egNR
); k
++)
1111 eee
[kk
++] = enerd
->grpp
.ener
[k
][gid
];
1114 add_ebin(md
->ebin
, md
->igrp
[n
], md
->nEc
, eee
, bSum
);
1122 for (i
= 0; (i
< md
->nTC
); i
++)
1124 md
->tmp_r
[i
] = ekind
->tcstat
[i
].T
;
1126 add_ebin(md
->ebin
, md
->itemp
, md
->nTC
, md
->tmp_r
, bSum
);
1128 if (md
->etc
== etcNOSEHOOVER
)
1130 /* whether to print Nose-Hoover chains: */
1131 if (md
->bPrintNHChains
)
1133 if (md
->bNHC_trotter
)
1135 for (i
= 0; (i
< md
->nTC
); i
++)
1137 for (j
= 0; j
< md
->nNHC
; j
++)
1140 md
->tmp_r
[2*k
] = state
->nosehoover_xi
[k
];
1141 md
->tmp_r
[2*k
+1] = state
->nosehoover_vxi
[k
];
1144 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1148 for (i
= 0; (i
< md
->nTCP
); i
++)
1150 for (j
= 0; j
< md
->nNHC
; j
++)
1153 md
->tmp_r
[2*k
] = state
->nhpres_xi
[k
];
1154 md
->tmp_r
[2*k
+1] = state
->nhpres_vxi
[k
];
1157 add_ebin(md
->ebin
, md
->itcb
, md
->mdeb_n
, md
->tmp_r
, bSum
);
1162 for (i
= 0; (i
< md
->nTC
); i
++)
1164 md
->tmp_r
[2*i
] = state
->nosehoover_xi
[i
];
1165 md
->tmp_r
[2*i
+1] = state
->nosehoover_vxi
[i
];
1167 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1171 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
1172 md
->etc
== etcVRESCALE
)
1174 for (i
= 0; (i
< md
->nTC
); i
++)
1176 md
->tmp_r
[i
] = ekind
->tcstat
[i
].lambda
;
1178 add_ebin(md
->ebin
, md
->itc
, md
->nTC
, md
->tmp_r
, bSum
);
1182 if (ekind
&& md
->nU
> 1)
1184 for (i
= 0; (i
< md
->nU
); i
++)
1186 copy_rvec(ekind
->grpstat
[i
].u
, md
->tmp_v
[i
]);
1188 add_ebin(md
->ebin
, md
->iu
, 3*md
->nU
, md
->tmp_v
[0], bSum
);
1191 ebin_increase_count(md
->ebin
, bSum
);
1193 /* BAR + thermodynamic integration values */
1194 if ((md
->fp_dhdl
|| md
->dhc
) && bDoDHDL
)
1196 for (i
= 0; i
< enerd
->n_lambda
-1; i
++)
1198 /* zero for simulated tempering */
1199 md
->dE
[i
] = enerd
->enerpart_lambda
[i
+1]-enerd
->enerpart_lambda
[0];
1200 if (md
->temperatures
!= NULL
)
1202 /* MRS: is this right, given the way we have defined the exchange probabilities? */
1203 /* is this even useful to have at all? */
1204 md
->dE
[i
] += (md
->temperatures
[i
]/
1205 md
->temperatures
[state
->fep_state
]-1.0)*
1206 enerd
->term
[F_EKIN
];
1212 fprintf(md
->fp_dhdl
, "%.4f", time
);
1213 /* the current free energy state */
1215 /* print the current state if we are doing expanded ensemble */
1216 if (expand
->elmcmove
> elmcmoveNO
)
1218 fprintf(md
->fp_dhdl
, " %4d", state
->fep_state
);
1220 /* total energy (for if the temperature changes */
1222 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
1224 switch (fep
->edHdLPrintEnergy
)
1226 case edHdLPrintEnergyPOTENTIAL
:
1227 store_energy
= enerd
->term
[F_EPOT
];
1229 case edHdLPrintEnergyTOTAL
:
1230 case edHdLPrintEnergyYES
:
1232 store_energy
= enerd
->term
[F_ETOT
];
1234 fprintf(md
->fp_dhdl
, " %#.8g", store_energy
);
1237 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
1239 for (i
= 0; i
< efptNR
; i
++)
1241 if (fep
->separate_dvdl
[i
])
1243 /* assumes F_DVDL is first */
1244 fprintf(md
->fp_dhdl
, " %#.8g", enerd
->term
[F_DVDL
+i
]);
1248 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
1250 fprintf(md
->fp_dhdl
, " %#.8g", md
->dE
[i
]);
1254 (md
->epc
!= epcNO
) &&
1255 (enerd
->n_lambda
> 0) &&
1256 (fep
->init_lambda
< 0))
1258 fprintf(md
->fp_dhdl
, " %#.8g", pv
); /* PV term only needed when
1259 there are alternate state
1260 lambda and we're not in
1261 compatibility mode */
1263 fprintf(md
->fp_dhdl
, "\n");
1264 /* and the binary free energy output */
1266 if (md
->dhc
&& bDoDHDL
)
1269 for (i
= 0; i
< efptNR
; i
++)
1271 if (fep
->separate_dvdl
[i
])
1273 /* assumes F_DVDL is first */
1274 store_dhdl
[idhdl
] = enerd
->term
[F_DVDL
+i
];
1278 store_energy
= enerd
->term
[F_ETOT
];
1279 /* store_dh is dE */
1280 mde_delta_h_coll_add_dh(md
->dhc
,
1281 (double)state
->fep_state
,
1285 md
->dE
+ fep
->lambda_start_n
,
1292 void upd_mdebin_step(t_mdebin
*md
)
1294 ebin_increase_count(md
->ebin
, FALSE
);
1297 static void npr(FILE *log
, int n
, char c
)
1299 for (; (n
> 0); n
--)
1301 fprintf(log
, "%c", c
);
1305 static void pprint(FILE *log
, const char *s
, t_mdebin
*md
)
1309 char buf1
[22], buf2
[22];
1312 fprintf(log
, "\t<====== ");
1313 npr(log
, slen
, CHAR
);
1314 fprintf(log
, " ==>\n");
1315 fprintf(log
, "\t<==== %s ====>\n", s
);
1316 fprintf(log
, "\t<== ");
1317 npr(log
, slen
, CHAR
);
1318 fprintf(log
, " ======>\n\n");
1320 fprintf(log
, "\tStatistics over %s steps using %s frames\n",
1321 gmx_step_str(md
->ebin
->nsteps_sim
, buf1
),
1322 gmx_step_str(md
->ebin
->nsum_sim
, buf2
));
1326 void print_ebin_header(FILE *log
, gmx_int64_t steps
, double time
, real lambda
)
1330 fprintf(log
, " %12s %12s %12s\n"
1331 " %12s %12.5f %12.5f\n\n",
1332 "Step", "Time", "Lambda", gmx_step_str(steps
, buf
), time
, lambda
);
1335 void print_ebin(ener_file_t fp_ene
, gmx_bool bEne
, gmx_bool bDR
, gmx_bool bOR
,
1337 gmx_int64_t step
, double time
,
1338 int mode
, gmx_bool bCompact
,
1339 t_mdebin
*md
, t_fcdata
*fcd
,
1340 gmx_groups_t
*groups
, t_grpopts
*opts
)
1342 /*static char **grpnms=NULL;*/
1344 int i
, j
, n
, ni
, nj
, b
;
1346 real
*disre_rm3tav
, *disre_rt
;
1348 /* these are for the old-style blocks (1 subblock, only reals), because
1349 there can be only one per ID for these */
1362 fr
.nsteps
= md
->ebin
->nsteps
;
1363 fr
.dt
= md
->delta_t
;
1364 fr
.nsum
= md
->ebin
->nsum
;
1365 fr
.nre
= (bEne
) ? md
->ebin
->nener
: 0;
1366 fr
.ener
= md
->ebin
->e
;
1367 ndisre
= bDR
? fcd
->disres
.npair
: 0;
1368 disre_rm3tav
= fcd
->disres
.rm3tav
;
1369 disre_rt
= fcd
->disres
.rt
;
1370 /* Optional additional old-style (real-only) blocks. */
1371 for (i
= 0; i
< enxNR
; i
++)
1375 if (fcd
->orires
.nr
> 0 && bOR
)
1377 diagonalize_orires_tensors(&(fcd
->orires
));
1378 nr
[enxOR
] = fcd
->orires
.nr
;
1379 block
[enxOR
] = fcd
->orires
.otav
;
1381 nr
[enxORI
] = (fcd
->orires
.oinsl
!= fcd
->orires
.otav
) ?
1383 block
[enxORI
] = fcd
->orires
.oinsl
;
1384 id
[enxORI
] = enxORI
;
1385 nr
[enxORT
] = fcd
->orires
.nex
*12;
1386 block
[enxORT
] = fcd
->orires
.eig
;
1387 id
[enxORT
] = enxORT
;
1390 /* whether we are going to wrte anything out: */
1391 if (fr
.nre
|| ndisre
|| nr
[enxOR
] || nr
[enxORI
])
1394 /* the old-style blocks go first */
1396 for (i
= 0; i
< enxNR
; i
++)
1403 add_blocks_enxframe(&fr
, fr
.nblock
);
1404 for (b
= 0; b
< fr
.nblock
; b
++)
1406 add_subblocks_enxblock(&(fr
.block
[b
]), 1);
1407 fr
.block
[b
].id
= id
[b
];
1408 fr
.block
[b
].sub
[0].nr
= nr
[b
];
1410 fr
.block
[b
].sub
[0].type
= xdr_datatype_float
;
1411 fr
.block
[b
].sub
[0].fval
= block
[b
];
1413 fr
.block
[b
].sub
[0].type
= xdr_datatype_double
;
1414 fr
.block
[b
].sub
[0].dval
= block
[b
];
1418 /* check for disre block & fill it. */
1423 add_blocks_enxframe(&fr
, fr
.nblock
);
1425 add_subblocks_enxblock(&(fr
.block
[db
]), 2);
1426 fr
.block
[db
].id
= enxDISRE
;
1427 fr
.block
[db
].sub
[0].nr
= ndisre
;
1428 fr
.block
[db
].sub
[1].nr
= ndisre
;
1430 fr
.block
[db
].sub
[0].type
= xdr_datatype_float
;
1431 fr
.block
[db
].sub
[1].type
= xdr_datatype_float
;
1432 fr
.block
[db
].sub
[0].fval
= disre_rt
;
1433 fr
.block
[db
].sub
[1].fval
= disre_rm3tav
;
1435 fr
.block
[db
].sub
[0].type
= xdr_datatype_double
;
1436 fr
.block
[db
].sub
[1].type
= xdr_datatype_double
;
1437 fr
.block
[db
].sub
[0].dval
= disre_rt
;
1438 fr
.block
[db
].sub
[1].dval
= disre_rm3tav
;
1441 /* here we can put new-style blocks */
1443 /* Free energy perturbation blocks */
1446 mde_delta_h_coll_handle_block(md
->dhc
, &fr
, fr
.nblock
);
1449 /* we can now free & reset the data in the blocks */
1452 mde_delta_h_coll_reset(md
->dhc
);
1455 /* do the actual I/O */
1456 do_enx(fp_ene
, &fr
);
1459 /* We have stored the sums, so reset the sum history */
1460 reset_ebin_sums(md
->ebin
);
1468 pprint(log
, "A V E R A G E S", md
);
1474 pprint(log
, "R M S - F L U C T U A T I O N S", md
);
1478 gmx_fatal(FARGS
, "Invalid print mode (%d)", mode
);
1483 for (i
= 0; i
< opts
->ngtc
; i
++)
1485 if (opts
->annealing
[i
] != eannNO
)
1487 fprintf(log
, "Current ref_t for group %s: %8.1f\n",
1488 *(groups
->grpname
[groups
->grps
[egcTC
].nm_ind
[i
]]),
1492 if (mode
== eprNORMAL
&& fcd
->orires
.nr
> 0)
1494 print_orires_log(log
, &(fcd
->orires
));
1496 fprintf(log
, " Energies (%s)\n", unit_energy
);
1497 pr_ebin(log
, md
->ebin
, md
->ie
, md
->f_nre
+md
->nCrmsd
, 5, mode
, TRUE
);
1504 pr_ebin(log
, md
->ebin
, md
->ib
, md
->bTricl
? NTRICLBOXS
: NBOXS
, 5,
1510 fprintf(log
, " Constraint Virial (%s)\n", unit_energy
);
1511 pr_ebin(log
, md
->ebin
, md
->isvir
, 9, 3, mode
, FALSE
);
1513 fprintf(log
, " Force Virial (%s)\n", unit_energy
);
1514 pr_ebin(log
, md
->ebin
, md
->ifvir
, 9, 3, mode
, FALSE
);
1519 fprintf(log
, " Total Virial (%s)\n", unit_energy
);
1520 pr_ebin(log
, md
->ebin
, md
->ivir
, 9, 3, mode
, FALSE
);
1525 fprintf(log
, " Pressure (%s)\n", unit_pres_bar
);
1526 pr_ebin(log
, md
->ebin
, md
->ipres
, 9, 3, mode
, FALSE
);
1531 fprintf(log
, " Total Dipole (%s)\n", unit_dipole_D
);
1532 pr_ebin(log
, md
->ebin
, md
->imu
, 3, 3, mode
, FALSE
);
1538 if (md
->print_grpnms
== NULL
)
1540 snew(md
->print_grpnms
, md
->nE
);
1542 for (i
= 0; (i
< md
->nEg
); i
++)
1544 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
1545 for (j
= i
; (j
< md
->nEg
); j
++)
1547 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
1548 sprintf(buf
, "%s-%s", *(groups
->grpname
[ni
]),
1549 *(groups
->grpname
[nj
]));
1550 md
->print_grpnms
[n
++] = gmx_strdup(buf
);
1554 sprintf(buf
, "Epot (%s)", unit_energy
);
1555 fprintf(log
, "%15s ", buf
);
1556 for (i
= 0; (i
< egNR
); i
++)
1560 fprintf(log
, "%12s ", egrp_nm
[i
]);
1564 for (i
= 0; (i
< md
->nE
); i
++)
1566 fprintf(log
, "%15s", md
->print_grpnms
[i
]);
1567 pr_ebin(log
, md
->ebin
, md
->igrp
[i
], md
->nEc
, md
->nEc
, mode
,
1574 pr_ebin(log
, md
->ebin
, md
->itemp
, md
->nTC
, 4, mode
, TRUE
);
1579 fprintf(log
, "%15s %12s %12s %12s\n",
1580 "Group", "Ux", "Uy", "Uz");
1581 for (i
= 0; (i
< md
->nU
); i
++)
1583 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
1584 fprintf(log
, "%15s", *groups
->grpname
[ni
]);
1585 pr_ebin(log
, md
->ebin
, md
->iu
+3*i
, 3, 3, mode
, FALSE
);
1594 void update_energyhistory(energyhistory_t
* enerhist
, t_mdebin
* mdebin
)
1598 enerhist
->nsteps
= mdebin
->ebin
->nsteps
;
1599 enerhist
->nsum
= mdebin
->ebin
->nsum
;
1600 enerhist
->nsteps_sim
= mdebin
->ebin
->nsteps_sim
;
1601 enerhist
->nsum_sim
= mdebin
->ebin
->nsum_sim
;
1602 enerhist
->nener
= mdebin
->ebin
->nener
;
1604 if (mdebin
->ebin
->nsum
> 0)
1606 /* Check if we need to allocate first */
1607 if (enerhist
->ener_ave
== NULL
)
1609 snew(enerhist
->ener_ave
, enerhist
->nener
);
1610 snew(enerhist
->ener_sum
, enerhist
->nener
);
1613 for (i
= 0; i
< enerhist
->nener
; i
++)
1615 enerhist
->ener_ave
[i
] = mdebin
->ebin
->e
[i
].eav
;
1616 enerhist
->ener_sum
[i
] = mdebin
->ebin
->e
[i
].esum
;
1620 if (mdebin
->ebin
->nsum_sim
> 0)
1622 /* Check if we need to allocate first */
1623 if (enerhist
->ener_sum_sim
== NULL
)
1625 snew(enerhist
->ener_sum_sim
, enerhist
->nener
);
1628 for (i
= 0; i
< enerhist
->nener
; i
++)
1630 enerhist
->ener_sum_sim
[i
] = mdebin
->ebin
->e_sim
[i
].esum
;
1635 mde_delta_h_coll_update_energyhistory(mdebin
->dhc
, enerhist
);
1639 void restore_energyhistory_from_state(t_mdebin
* mdebin
,
1640 energyhistory_t
* enerhist
)
1644 if ((enerhist
->nsum
> 0 || enerhist
->nsum_sim
> 0) &&
1645 mdebin
->ebin
->nener
!= enerhist
->nener
)
1647 gmx_fatal(FARGS
, "Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1648 mdebin
->ebin
->nener
, enerhist
->nener
);
1651 mdebin
->ebin
->nsteps
= enerhist
->nsteps
;
1652 mdebin
->ebin
->nsum
= enerhist
->nsum
;
1653 mdebin
->ebin
->nsteps_sim
= enerhist
->nsteps_sim
;
1654 mdebin
->ebin
->nsum_sim
= enerhist
->nsum_sim
;
1656 for (i
= 0; i
< mdebin
->ebin
->nener
; i
++)
1658 mdebin
->ebin
->e
[i
].eav
=
1659 (enerhist
->nsum
> 0 ? enerhist
->ener_ave
[i
] : 0);
1660 mdebin
->ebin
->e
[i
].esum
=
1661 (enerhist
->nsum
> 0 ? enerhist
->ener_sum
[i
] : 0);
1662 mdebin
->ebin
->e_sim
[i
].esum
=
1663 (enerhist
->nsum_sim
> 0 ? enerhist
->ener_sum_sim
[i
] : 0);
1667 mde_delta_h_coll_restore_energyhistory(mdebin
->dhc
, enerhist
);