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/fcdata.h"
55 #include "gromacs/legacyheaders/types/group.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/constr.h"
59 #include "gromacs/mdlib/mdebin_bar.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/pulling/pull.h"
62 #include "gromacs/topology/mtop_util.h"
63 #include "gromacs/utility/arraysize.h"
64 #include "gromacs/utility/fatalerror.h"
65 #include "gromacs/utility/smalloc.h"
67 static const char *conrmsd_nm
[] = { "Constr. rmsd", "Constr.2 rmsd" };
69 static const char *boxs_nm
[] = { "Box-X", "Box-Y", "Box-Z" };
71 static const char *tricl_boxs_nm
[] = {
72 "Box-XX", "Box-YY", "Box-ZZ",
73 "Box-YX", "Box-ZX", "Box-ZY"
76 static const char *vol_nm
[] = { "Volume" };
78 static const char *dens_nm
[] = {"Density" };
80 static const char *pv_nm
[] = {"pV" };
82 static const char *enthalpy_nm
[] = {"Enthalpy" };
84 static const char *boxvel_nm
[] = {
85 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
86 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
89 #define NBOXS asize(boxs_nm)
90 #define NTRICLBOXS asize(tricl_boxs_nm)
92 t_mdebin
*init_mdebin(ener_file_t fp_ene
,
93 const gmx_mtop_t
*mtop
,
97 const char *ener_nm
[F_NRE
];
98 static const char *vir_nm
[] = {
99 "Vir-XX", "Vir-XY", "Vir-XZ",
100 "Vir-YX", "Vir-YY", "Vir-YZ",
101 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
103 static const char *sv_nm
[] = {
104 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
105 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
106 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
108 static const char *fv_nm
[] = {
109 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
110 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
111 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
113 static const char *pres_nm
[] = {
114 "Pres-XX", "Pres-XY", "Pres-XZ",
115 "Pres-YX", "Pres-YY", "Pres-YZ",
116 "Pres-ZX", "Pres-ZY", "Pres-ZZ"
118 static const char *surft_nm
[] = {
121 static const char *mu_nm
[] = {
122 "Mu-X", "Mu-Y", "Mu-Z"
124 static const char *vcos_nm
[] = {
127 static const char *visc_nm
[] = {
130 static const char *baro_nm
[] = {
135 const gmx_groups_t
*groups
;
140 int i
, j
, ni
, nj
, n
, k
, kk
, ncon
, nset
;
150 if (EI_DYNAMICS(ir
->eI
))
152 md
->delta_t
= ir
->delta_t
;
159 groups
= &mtop
->groups
;
161 bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
162 b14
= (gmx_mtop_ftype_count(mtop
, F_LJ14
) > 0 ||
163 gmx_mtop_ftype_count(mtop
, F_LJC14_Q
) > 0);
165 ncon
= gmx_mtop_ftype_count(mtop
, F_CONSTR
);
166 nset
= gmx_mtop_ftype_count(mtop
, F_SETTLE
);
167 md
->bConstr
= (ncon
> 0 || nset
> 0);
168 md
->bConstrVir
= FALSE
;
171 if (ncon
> 0 && ir
->eConstrAlg
== econtLINCS
)
182 md
->bConstrVir
= (getenv("GMX_CONSTRAINTVIR") != NULL
);
189 /* Energy monitoring */
190 for (i
= 0; i
< egNR
; i
++)
192 md
->bEInd
[i
] = FALSE
;
195 for (i
= 0; i
< F_NRE
; i
++)
197 md
->bEner
[i
] = FALSE
;
200 md
->bEner
[i
] = !bBHAM
;
202 else if (i
== F_BHAM
)
204 md
->bEner
[i
] = bBHAM
;
208 md
->bEner
[i
] = ir
->bQMMM
;
210 else if (i
== F_COUL_LR
)
212 md
->bEner
[i
] = (ir
->rcoulomb
> ir
->rlist
);
214 else if (i
== F_LJ_LR
)
216 md
->bEner
[i
] = (!bBHAM
&& ir
->rvdw
> ir
->rlist
);
218 else if (i
== F_BHAM_LR
)
220 md
->bEner
[i
] = (bBHAM
&& ir
->rvdw
> ir
->rlist
);
222 else if (i
== F_RF_EXCL
)
224 md
->bEner
[i
] = (EEL_RF(ir
->coulombtype
) && ir
->coulombtype
!= eelRF_NEC
&& ir
->cutoff_scheme
== ecutsGROUP
);
226 else if (i
== F_COUL_RECIP
)
228 md
->bEner
[i
] = EEL_FULL(ir
->coulombtype
);
230 else if (i
== F_LJ_RECIP
)
232 md
->bEner
[i
] = EVDW_PME(ir
->vdwtype
);
234 else if (i
== F_LJ14
)
238 else if (i
== F_COUL14
)
242 else if (i
== F_LJC14_Q
|| i
== F_LJC_PAIRS_NB
)
244 md
->bEner
[i
] = FALSE
;
246 else if ((i
== F_DVDL_COUL
&& ir
->fepvals
->separate_dvdl
[efptCOUL
]) ||
247 (i
== F_DVDL_VDW
&& ir
->fepvals
->separate_dvdl
[efptVDW
]) ||
248 (i
== F_DVDL_BONDED
&& ir
->fepvals
->separate_dvdl
[efptBONDED
]) ||
249 (i
== F_DVDL_RESTRAINT
&& ir
->fepvals
->separate_dvdl
[efptRESTRAINT
]) ||
250 (i
== F_DKDL
&& ir
->fepvals
->separate_dvdl
[efptMASS
]) ||
251 (i
== F_DVDL
&& ir
->fepvals
->separate_dvdl
[efptFEP
]))
253 md
->bEner
[i
] = (ir
->efep
!= efepNO
);
255 else if ((interaction_function
[i
].flags
& IF_VSITE
) ||
256 (i
== F_CONSTR
) || (i
== F_CONSTRNC
) || (i
== F_SETTLE
))
258 md
->bEner
[i
] = FALSE
;
260 else if ((i
== F_COUL_SR
) || (i
== F_EPOT
) || (i
== F_PRES
) || (i
== F_EQM
))
264 else if ((i
== F_GBPOL
) && ir
->implicit_solvent
== eisGBSA
)
268 else if ((i
== F_NPSOLVATION
) && ir
->implicit_solvent
== eisGBSA
&& (ir
->sa_algorithm
!= esaNO
))
272 else if ((i
== F_GB12
) || (i
== F_GB13
) || (i
== F_GB14
))
274 md
->bEner
[i
] = FALSE
;
276 else if ((i
== F_ETOT
) || (i
== F_EKIN
) || (i
== F_TEMP
))
278 md
->bEner
[i
] = EI_DYNAMICS(ir
->eI
);
280 else if (i
== F_DISPCORR
|| i
== F_PDISPCORR
)
282 md
->bEner
[i
] = (ir
->eDispCorr
!= edispcNO
);
284 else if (i
== F_DISRESVIOL
)
286 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_DISRES
) > 0);
288 else if (i
== F_ORIRESDEV
)
290 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, F_ORIRES
) > 0);
292 else if (i
== F_CONNBONDS
)
294 md
->bEner
[i
] = FALSE
;
296 else if (i
== F_COM_PULL
)
298 md
->bEner
[i
] = (ir
->bPull
&& pull_have_potential(ir
->pull_work
));
300 else if (i
== F_ECONSERVED
)
302 md
->bEner
[i
] = ((ir
->etc
== etcNOSEHOOVER
|| ir
->etc
== etcVRESCALE
) &&
303 (ir
->epc
== epcNO
|| ir
->epc
== epcMTTK
));
307 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
, i
) > 0);
311 /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
312 if (ir
->bAdress
&& !debug
)
314 for (i
= 0; i
< F_NRE
; i
++)
316 md
->bEner
[i
] = FALSE
;
333 for (i
= 0; i
< F_NRE
; i
++)
337 ener_nm
[md
->f_nre
] = interaction_function
[i
].longname
;
343 md
->bDiagPres
= !TRICLINIC(ir
->ref_p
);
344 md
->ref_p
= (ir
->ref_p
[XX
][XX
]+ir
->ref_p
[YY
][YY
]+ir
->ref_p
[ZZ
][ZZ
])/DIM
;
345 md
->bTricl
= TRICLINIC(ir
->compress
) || TRICLINIC(ir
->deform
);
346 md
->bDynBox
= DYNAMIC_BOX(*ir
);
348 md
->bNHC_trotter
= IR_NVT_TROTTER(ir
);
349 md
->bPrintNHChains
= ir
->bPrintNHChains
;
350 md
->bMTTK
= (IR_NPT_TROTTER(ir
) || IR_NPH_TROTTER(ir
));
351 md
->bMu
= NEED_MUTOT(*ir
);
353 md
->ebin
= mk_ebin();
354 /* Pass NULL for unit to let get_ebin_space determine the units
355 * for interaction_function[i].longname
357 md
->ie
= get_ebin_space(md
->ebin
, md
->f_nre
, ener_nm
, NULL
);
360 /* This should be called directly after the call for md->ie,
361 * such that md->iconrmsd follows directly in the list.
363 md
->iconrmsd
= get_ebin_space(md
->ebin
, md
->nCrmsd
, conrmsd_nm
, "");
367 md
->ib
= get_ebin_space(md
->ebin
,
368 md
->bTricl
? NTRICLBOXS
: NBOXS
,
369 md
->bTricl
? tricl_boxs_nm
: boxs_nm
,
371 md
->ivol
= get_ebin_space(md
->ebin
, 1, vol_nm
, unit_volume
);
372 md
->idens
= get_ebin_space(md
->ebin
, 1, dens_nm
, unit_density_SI
);
375 md
->ipv
= get_ebin_space(md
->ebin
, 1, pv_nm
, unit_energy
);
376 md
->ienthalpy
= get_ebin_space(md
->ebin
, 1, enthalpy_nm
, unit_energy
);
381 md
->isvir
= get_ebin_space(md
->ebin
, asize(sv_nm
), sv_nm
, unit_energy
);
382 md
->ifvir
= get_ebin_space(md
->ebin
, asize(fv_nm
), fv_nm
, unit_energy
);
386 md
->ivir
= get_ebin_space(md
->ebin
, asize(vir_nm
), vir_nm
, unit_energy
);
390 md
->ipres
= get_ebin_space(md
->ebin
, asize(pres_nm
), pres_nm
, unit_pres_bar
);
394 md
->isurft
= get_ebin_space(md
->ebin
, asize(surft_nm
), surft_nm
,
397 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
399 md
->ipc
= get_ebin_space(md
->ebin
, md
->bTricl
? 6 : 3,
400 boxvel_nm
, unit_vel
);
404 md
->imu
= get_ebin_space(md
->ebin
, asize(mu_nm
), mu_nm
, unit_dipole_D
);
406 if (ir
->cos_accel
!= 0)
408 md
->ivcos
= get_ebin_space(md
->ebin
, asize(vcos_nm
), vcos_nm
, unit_vel
);
409 md
->ivisc
= get_ebin_space(md
->ebin
, asize(visc_nm
), visc_nm
,
413 /* Energy monitoring */
414 for (i
= 0; i
< egNR
; i
++)
416 md
->bEInd
[i
] = FALSE
;
418 md
->bEInd
[egCOULSR
] = TRUE
;
419 md
->bEInd
[egLJSR
] = TRUE
;
421 if (ir
->rcoulomb
> ir
->rlist
)
423 md
->bEInd
[egCOULLR
] = TRUE
;
427 if (ir
->rvdw
> ir
->rlist
)
429 md
->bEInd
[egLJLR
] = TRUE
;
434 md
->bEInd
[egLJSR
] = FALSE
;
435 md
->bEInd
[egBHAMSR
] = TRUE
;
436 if (ir
->rvdw
> ir
->rlist
)
438 md
->bEInd
[egBHAMLR
] = TRUE
;
443 md
->bEInd
[egLJ14
] = TRUE
;
444 md
->bEInd
[egCOUL14
] = TRUE
;
447 for (i
= 0; (i
< egNR
); i
++)
455 n
= groups
->grps
[egcENER
].nr
;
456 /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
459 /*standard simulation*/
461 md
->nE
= (n
*(n
+1))/2;
465 /*AdResS simulation*/
472 snew(md
->igrp
, md
->nE
);
477 for (k
= 0; (k
< md
->nEc
); k
++)
479 snew(gnm
[k
], STRLEN
);
481 for (i
= 0; (i
< groups
->grps
[egcENER
].nr
); i
++)
483 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
484 for (j
= i
; (j
< groups
->grps
[egcENER
].nr
); j
++)
486 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
487 for (k
= kk
= 0; (k
< egNR
); k
++)
491 sprintf(gnm
[kk
], "%s:%s-%s", egrp_nm
[k
],
492 *(groups
->grpname
[ni
]), *(groups
->grpname
[nj
]));
496 md
->igrp
[n
] = get_ebin_space(md
->ebin
, md
->nEc
,
497 (const char **)gnm
, unit_energy
);
501 for (k
= 0; (k
< md
->nEc
); k
++)
509 gmx_incons("Number of energy terms wrong");
513 md
->nTC
= groups
->grps
[egcTC
].nr
;
514 md
->nNHC
= ir
->opts
.nhchainlength
; /* shorthand for number of NH chains */
517 md
->nTCP
= 1; /* assume only one possible coupling system for barostat
524 if (md
->etc
== etcNOSEHOOVER
)
526 if (md
->bNHC_trotter
)
528 md
->mde_n
= 2*md
->nNHC
*md
->nTC
;
532 md
->mde_n
= 2*md
->nTC
;
534 if (md
->epc
== epcMTTK
)
536 md
->mdeb_n
= 2*md
->nNHC
*md
->nTCP
;
545 snew(md
->tmp_r
, md
->mde_n
);
546 snew(md
->tmp_v
, md
->mde_n
);
547 snew(md
->grpnms
, md
->mde_n
);
550 for (i
= 0; (i
< md
->nTC
); i
++)
552 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
553 sprintf(buf
, "T-%s", *(groups
->grpname
[ni
]));
554 grpnms
[i
] = gmx_strdup(buf
);
556 md
->itemp
= get_ebin_space(md
->ebin
, md
->nTC
, (const char **)grpnms
,
559 if (md
->etc
== etcNOSEHOOVER
)
561 if (md
->bPrintNHChains
)
563 if (md
->bNHC_trotter
)
565 for (i
= 0; (i
< md
->nTC
); i
++)
567 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
568 bufi
= *(groups
->grpname
[ni
]);
569 for (j
= 0; (j
< md
->nNHC
); j
++)
571 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
572 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
573 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
574 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
577 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
578 (const char **)grpnms
, unit_invtime
);
581 for (i
= 0; (i
< md
->nTCP
); i
++)
583 bufi
= baro_nm
[0]; /* All barostat DOF's together for now. */
584 for (j
= 0; (j
< md
->nNHC
); j
++)
586 sprintf(buf
, "Xi-%d-%s", j
, bufi
);
587 grpnms
[2*(i
*md
->nNHC
+j
)] = gmx_strdup(buf
);
588 sprintf(buf
, "vXi-%d-%s", j
, bufi
);
589 grpnms
[2*(i
*md
->nNHC
+j
)+1] = gmx_strdup(buf
);
592 md
->itcb
= get_ebin_space(md
->ebin
, md
->mdeb_n
,
593 (const char **)grpnms
, unit_invtime
);
598 for (i
= 0; (i
< md
->nTC
); i
++)
600 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
601 bufi
= *(groups
->grpname
[ni
]);
602 sprintf(buf
, "Xi-%s", bufi
);
603 grpnms
[2*i
] = gmx_strdup(buf
);
604 sprintf(buf
, "vXi-%s", bufi
);
605 grpnms
[2*i
+1] = gmx_strdup(buf
);
607 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
,
608 (const char **)grpnms
, unit_invtime
);
612 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
613 md
->etc
== etcVRESCALE
)
615 for (i
= 0; (i
< md
->nTC
); i
++)
617 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
618 sprintf(buf
, "Lamb-%s", *(groups
->grpname
[ni
]));
619 grpnms
[i
] = gmx_strdup(buf
);
621 md
->itc
= get_ebin_space(md
->ebin
, md
->mde_n
, (const char **)grpnms
, "");
627 md
->nU
= groups
->grps
[egcACC
].nr
;
630 snew(grpnms
, 3*md
->nU
);
631 for (i
= 0; (i
< md
->nU
); i
++)
633 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
634 sprintf(buf
, "Ux-%s", *(groups
->grpname
[ni
]));
635 grpnms
[3*i
+XX
] = gmx_strdup(buf
);
636 sprintf(buf
, "Uy-%s", *(groups
->grpname
[ni
]));
637 grpnms
[3*i
+YY
] = gmx_strdup(buf
);
638 sprintf(buf
, "Uz-%s", *(groups
->grpname
[ni
]));
639 grpnms
[3*i
+ZZ
] = gmx_strdup(buf
);
641 md
->iu
= get_ebin_space(md
->ebin
, 3*md
->nU
, (const char **)grpnms
, unit_vel
);
647 do_enxnms(fp_ene
, &md
->ebin
->nener
, &md
->ebin
->enm
);
650 md
->print_grpnms
= NULL
;
652 /* check whether we're going to write dh histograms */
654 if (ir
->fepvals
->separate_dhdl_file
== esepdhdlfileNO
)
656 /* Currently dh histograms are only written with dynamics */
657 if (EI_DYNAMICS(ir
->eI
))
661 mde_delta_h_coll_init(md
->dhc
, ir
);
664 snew(md
->dE
, ir
->fepvals
->n_lambda
);
668 md
->fp_dhdl
= fp_dhdl
;
669 snew(md
->dE
, ir
->fepvals
->n_lambda
);
674 snew(md
->temperatures
, ir
->fepvals
->n_lambda
);
675 for (i
= 0; i
< ir
->fepvals
->n_lambda
; i
++)
677 md
->temperatures
[i
] = ir
->simtempvals
->temperatures
[i
];
683 /* print a lambda vector to a string
684 fep = the inputrec's FEP input data
685 i = the index of the lambda vector
686 get_native_lambda = whether to print the native lambda
687 get_names = whether to print the names rather than the values
688 str = the pre-allocated string buffer to print to. */
689 static void print_lambda_vector(t_lambda
*fep
, int i
,
690 gmx_bool get_native_lambda
, gmx_bool get_names
,
696 for (j
= 0; j
< efptNR
; j
++)
698 if (fep
->separate_dvdl
[j
])
703 str
[0] = 0; /* reset the string */
706 str
+= sprintf(str
, "("); /* set the opening parenthesis*/
708 for (j
= 0; j
< efptNR
; j
++)
710 if (fep
->separate_dvdl
[j
])
714 if (get_native_lambda
&& fep
->init_lambda
>= 0)
716 str
+= sprintf(str
, "%.4f", fep
->init_lambda
);
720 str
+= sprintf(str
, "%.4f", fep
->all_lambda
[j
][i
]);
725 str
+= sprintf(str
, "%s", efpt_singular_names
[j
]);
727 /* print comma for the next item */
730 str
+= sprintf(str
, ", ");
737 /* and add the closing parenthesis */
743 extern FILE *open_dhdl(const char *filename
, const t_inputrec
*ir
,
744 const gmx_output_env_t
*oenv
)
747 const char *dhdl
= "dH/d\\lambda", *deltag
= "\\DeltaH", *lambda
= "\\lambda",
748 *lambdastate
= "\\lambda state";
749 char title
[STRLEN
], label_x
[STRLEN
], label_y
[STRLEN
];
750 int i
, nps
, nsets
, nsets_de
, nsetsbegin
;
751 int n_lambda_terms
= 0;
752 t_lambda
*fep
= ir
->fepvals
; /* for simplicity */
753 t_expanded
*expand
= ir
->expandedvals
;
755 char buf
[STRLEN
], lambda_vec_str
[STRLEN
], lambda_name_str
[STRLEN
];
761 gmx_bool write_pV
= FALSE
;
763 /* count the number of different lambda terms */
764 for (i
= 0; i
< efptNR
; i
++)
766 if (fep
->separate_dvdl
[i
])
772 if (fep
->n_lambda
== 0)
774 sprintf(title
, "%s", dhdl
);
775 sprintf(label_x
, "Time (ps)");
776 sprintf(label_y
, "%s (%s %s)",
777 dhdl
, unit_energy
, "[\\lambda]\\S-1\\N");
781 sprintf(title
, "%s and %s", dhdl
, deltag
);
782 sprintf(label_x
, "Time (ps)");
783 sprintf(label_y
, "%s and %s (%s %s)",
784 dhdl
, deltag
, unit_energy
, "[\\8l\\4]\\S-1\\N");
786 fp
= gmx_fio_fopen(filename
, "w+");
787 xvgr_header(fp
, title
, label_x
, label_y
, exvggtXNY
, oenv
);
791 bufplace
= sprintf(buf
, "T = %g (K) ",
794 if ((ir
->efep
!= efepSLOWGROWTH
) && (ir
->efep
!= efepEXPANDED
))
796 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
798 /* compatibility output */
799 sprintf(&(buf
[bufplace
]), "%s = %.4f", lambda
, fep
->init_lambda
);
803 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, FALSE
,
805 print_lambda_vector(fep
, fep
->init_fep_state
, TRUE
, TRUE
,
807 sprintf(&(buf
[bufplace
]), "%s %d: %s = %s",
808 lambdastate
, fep
->init_fep_state
,
809 lambda_name_str
, lambda_vec_str
);
812 xvgr_subtitle(fp
, buf
, oenv
);
816 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
818 nsets_dhdl
= n_lambda_terms
;
820 /* count the number of delta_g states */
821 nsets_de
= fep
->lambda_stop_n
- fep
->lambda_start_n
;
823 nsets
= nsets_dhdl
+ nsets_de
; /* dhdl + fep differences */
825 if (fep
->n_lambda
> 0 && (expand
->elmcmove
> elmcmoveNO
))
827 nsets
+= 1; /*add fep state for expanded ensemble */
830 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
832 nsets
+= 1; /* add energy to the dhdl as well */
836 if ((ir
->epc
!= epcNO
) && (fep
->n_lambda
> 0) && (fep
->init_lambda
< 0))
838 nsetsextend
+= 1; /* for PV term, other terms possible if required for
839 the reduced potential (only needed with foreign
840 lambda, and only output when init_lambda is not
841 set in order to maintain compatibility of the
845 snew(setname
, nsetsextend
);
847 if (expand
->elmcmove
> elmcmoveNO
)
849 /* state for the fep_vals, if we have alchemical sampling */
850 sprintf(buf
, "%s", "Thermodynamic state");
851 setname
[s
] = gmx_strdup(buf
);
855 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
857 switch (fep
->edHdLPrintEnergy
)
859 case edHdLPrintEnergyPOTENTIAL
:
860 sprintf(buf
, "%s (%s)", "Potential Energy", unit_energy
);
862 case edHdLPrintEnergyTOTAL
:
863 case edHdLPrintEnergyYES
:
865 sprintf(buf
, "%s (%s)", "Total Energy", unit_energy
);
867 setname
[s
] = gmx_strdup(buf
);
871 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
873 for (i
= 0; i
< efptNR
; i
++)
875 if (fep
->separate_dvdl
[i
])
878 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
880 /* compatibility output */
881 sprintf(buf
, "%s %s %.4f", dhdl
, lambda
, fep
->init_lambda
);
885 double lam
= fep
->init_lambda
;
886 if (fep
->init_lambda
< 0)
888 lam
= fep
->all_lambda
[i
][fep
->init_fep_state
];
890 sprintf(buf
, "%s %s = %.4f", dhdl
, efpt_singular_names
[i
],
893 setname
[s
] = gmx_strdup(buf
);
899 if (fep
->n_lambda
> 0)
901 /* g_bar has to determine the lambda values used in this simulation
902 * from this xvg legend.
905 if (expand
->elmcmove
> elmcmoveNO
)
907 nsetsbegin
= 1; /* for including the expanded ensemble */
914 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
918 nsetsbegin
+= nsets_dhdl
;
920 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
922 print_lambda_vector(fep
, i
, FALSE
, FALSE
, lambda_vec_str
);
923 if ( (fep
->init_lambda
>= 0) && (n_lambda_terms
== 1 ))
925 /* for compatible dhdl.xvg files */
926 nps
= sprintf(buf
, "%s %s %s", deltag
, lambda
, lambda_vec_str
);
930 nps
= sprintf(buf
, "%s %s to %s", deltag
, lambda
, lambda_vec_str
);
935 /* print the temperature for this state if doing simulated annealing */
936 sprintf(&buf
[nps
], "T = %g (%s)",
937 ir
->simtempvals
->temperatures
[s
-(nsetsbegin
)],
940 setname
[s
] = gmx_strdup(buf
);
945 sprintf(buf
, "pV (%s)", unit_energy
);
946 setname
[nsetsextend
-1] = gmx_strdup(buf
); /* the first entry after
950 xvgr_legend(fp
, nsetsextend
, (const char **)setname
, oenv
);
952 for (s
= 0; s
< nsetsextend
; s
++)
962 static void copy_energy(t_mdebin
*md
, real e
[], real ecpy
[])
966 for (i
= j
= 0; (i
< F_NRE
); i
++)
975 gmx_incons("Number of energy terms wrong");
979 void upd_mdebin(t_mdebin
*md
,
984 gmx_enerdata_t
*enerd
,
993 gmx_ekindata_t
*ekind
,
997 int i
, j
, k
, kk
, n
, gid
;
998 real crmsd
[2], tmp6
[6];
999 real bs
[NTRICLBOXS
], vol
, dens
, pv
, enthalpy
;
1002 double store_dhdl
[efptNR
];
1003 real store_energy
= 0;
1006 /* Do NOT use the box in the state variable, but the separate box provided
1007 * as an argument. This is because we sometimes need to write the box from
1008 * the last timestep to match the trajectory frames.
1010 copy_energy(md
, enerd
->term
, ecopy
);
1011 add_ebin(md
->ebin
, md
->ie
, md
->f_nre
, ecopy
, bSum
);
1014 crmsd
[0] = constr_rmsd(constr
, FALSE
);
1017 crmsd
[1] = constr_rmsd(constr
, TRUE
);
1019 add_ebin(md
->ebin
, md
->iconrmsd
, md
->nCrmsd
, crmsd
, FALSE
);
1026 bs
[0] = box
[XX
][XX
];
1027 bs
[1] = box
[YY
][YY
];
1028 bs
[2] = box
[ZZ
][ZZ
];
1029 bs
[3] = box
[YY
][XX
];
1030 bs
[4] = box
[ZZ
][XX
];
1031 bs
[5] = box
[ZZ
][YY
];
1036 bs
[0] = box
[XX
][XX
];
1037 bs
[1] = box
[YY
][YY
];
1038 bs
[2] = box
[ZZ
][ZZ
];
1041 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
1042 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
1043 add_ebin(md
->ebin
, md
->ib
, nboxs
, bs
, bSum
);
1044 add_ebin(md
->ebin
, md
->ivol
, 1, &vol
, bSum
);
1045 add_ebin(md
->ebin
, md
->idens
, 1, &dens
, bSum
);
1049 /* This is pV (in kJ/mol). The pressure is the reference pressure,
1050 not the instantaneous pressure */
1051 pv
= vol
*md
->ref_p
/PRESFAC
;
1053 add_ebin(md
->ebin
, md
->ipv
, 1, &pv
, bSum
);
1054 enthalpy
= pv
+ enerd
->term
[F_ETOT
];
1055 add_ebin(md
->ebin
, md
->ienthalpy
, 1, &enthalpy
, bSum
);
1060 add_ebin(md
->ebin
, md
->isvir
, 9, svir
[0], bSum
);
1061 add_ebin(md
->ebin
, md
->ifvir
, 9, fvir
[0], bSum
);
1065 add_ebin(md
->ebin
, md
->ivir
, 9, vir
[0], bSum
);
1069 add_ebin(md
->ebin
, md
->ipres
, 9, pres
[0], bSum
);
1073 tmp
= (pres
[ZZ
][ZZ
]-(pres
[XX
][XX
]+pres
[YY
][YY
])*0.5)*box
[ZZ
][ZZ
];
1074 add_ebin(md
->ebin
, md
->isurft
, 1, &tmp
, bSum
);
1076 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
1078 tmp6
[0] = state
->boxv
[XX
][XX
];
1079 tmp6
[1] = state
->boxv
[YY
][YY
];
1080 tmp6
[2] = state
->boxv
[ZZ
][ZZ
];
1081 tmp6
[3] = state
->boxv
[YY
][XX
];
1082 tmp6
[4] = state
->boxv
[ZZ
][XX
];
1083 tmp6
[5] = state
->boxv
[ZZ
][YY
];
1084 add_ebin(md
->ebin
, md
->ipc
, md
->bTricl
? 6 : 3, tmp6
, bSum
);
1088 add_ebin(md
->ebin
, md
->imu
, 3, mu_tot
, bSum
);
1090 if (ekind
&& ekind
->cosacc
.cos_accel
!= 0)
1092 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
1093 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
1094 add_ebin(md
->ebin
, md
->ivcos
, 1, &(ekind
->cosacc
.vcos
), bSum
);
1095 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
1096 tmp
= 1/(ekind
->cosacc
.cos_accel
/(ekind
->cosacc
.vcos
*PICO
)
1097 *dens
*sqr(box
[ZZ
][ZZ
]*NANO
/(2*M_PI
)));
1098 add_ebin(md
->ebin
, md
->ivisc
, 1, &tmp
, bSum
);
1103 for (i
= 0; (i
< md
->nEg
); i
++)
1105 for (j
= i
; (j
< md
->nEg
); j
++)
1107 gid
= GID(i
, j
, md
->nEg
);
1108 for (k
= kk
= 0; (k
< egNR
); k
++)
1112 eee
[kk
++] = enerd
->grpp
.ener
[k
][gid
];
1115 add_ebin(md
->ebin
, md
->igrp
[n
], md
->nEc
, eee
, bSum
);
1123 for (i
= 0; (i
< md
->nTC
); i
++)
1125 md
->tmp_r
[i
] = ekind
->tcstat
[i
].T
;
1127 add_ebin(md
->ebin
, md
->itemp
, md
->nTC
, md
->tmp_r
, bSum
);
1129 if (md
->etc
== etcNOSEHOOVER
)
1131 /* whether to print Nose-Hoover chains: */
1132 if (md
->bPrintNHChains
)
1134 if (md
->bNHC_trotter
)
1136 for (i
= 0; (i
< md
->nTC
); i
++)
1138 for (j
= 0; j
< md
->nNHC
; j
++)
1141 md
->tmp_r
[2*k
] = state
->nosehoover_xi
[k
];
1142 md
->tmp_r
[2*k
+1] = state
->nosehoover_vxi
[k
];
1145 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1149 for (i
= 0; (i
< md
->nTCP
); i
++)
1151 for (j
= 0; j
< md
->nNHC
; j
++)
1154 md
->tmp_r
[2*k
] = state
->nhpres_xi
[k
];
1155 md
->tmp_r
[2*k
+1] = state
->nhpres_vxi
[k
];
1158 add_ebin(md
->ebin
, md
->itcb
, md
->mdeb_n
, md
->tmp_r
, bSum
);
1163 for (i
= 0; (i
< md
->nTC
); i
++)
1165 md
->tmp_r
[2*i
] = state
->nosehoover_xi
[i
];
1166 md
->tmp_r
[2*i
+1] = state
->nosehoover_vxi
[i
];
1168 add_ebin(md
->ebin
, md
->itc
, md
->mde_n
, md
->tmp_r
, bSum
);
1172 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
1173 md
->etc
== etcVRESCALE
)
1175 for (i
= 0; (i
< md
->nTC
); i
++)
1177 md
->tmp_r
[i
] = ekind
->tcstat
[i
].lambda
;
1179 add_ebin(md
->ebin
, md
->itc
, md
->nTC
, md
->tmp_r
, bSum
);
1183 if (ekind
&& md
->nU
> 1)
1185 for (i
= 0; (i
< md
->nU
); i
++)
1187 copy_rvec(ekind
->grpstat
[i
].u
, md
->tmp_v
[i
]);
1189 add_ebin(md
->ebin
, md
->iu
, 3*md
->nU
, md
->tmp_v
[0], bSum
);
1192 ebin_increase_count(md
->ebin
, bSum
);
1194 /* BAR + thermodynamic integration values */
1195 if ((md
->fp_dhdl
|| md
->dhc
) && bDoDHDL
)
1197 for (i
= 0; i
< enerd
->n_lambda
-1; i
++)
1199 /* zero for simulated tempering */
1200 md
->dE
[i
] = enerd
->enerpart_lambda
[i
+1]-enerd
->enerpart_lambda
[0];
1201 if (md
->temperatures
!= NULL
)
1203 /* MRS: is this right, given the way we have defined the exchange probabilities? */
1204 /* is this even useful to have at all? */
1205 md
->dE
[i
] += (md
->temperatures
[i
]/
1206 md
->temperatures
[state
->fep_state
]-1.0)*
1207 enerd
->term
[F_EKIN
];
1213 fprintf(md
->fp_dhdl
, "%.4f", time
);
1214 /* the current free energy state */
1216 /* print the current state if we are doing expanded ensemble */
1217 if (expand
->elmcmove
> elmcmoveNO
)
1219 fprintf(md
->fp_dhdl
, " %4d", state
->fep_state
);
1221 /* total energy (for if the temperature changes */
1223 if (fep
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
1225 switch (fep
->edHdLPrintEnergy
)
1227 case edHdLPrintEnergyPOTENTIAL
:
1228 store_energy
= enerd
->term
[F_EPOT
];
1230 case edHdLPrintEnergyTOTAL
:
1231 case edHdLPrintEnergyYES
:
1233 store_energy
= enerd
->term
[F_ETOT
];
1235 fprintf(md
->fp_dhdl
, " %#.8g", store_energy
);
1238 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
1240 for (i
= 0; i
< efptNR
; i
++)
1242 if (fep
->separate_dvdl
[i
])
1244 /* assumes F_DVDL is first */
1245 fprintf(md
->fp_dhdl
, " %#.8g", enerd
->term
[F_DVDL
+i
]);
1249 for (i
= fep
->lambda_start_n
; i
< fep
->lambda_stop_n
; i
++)
1251 fprintf(md
->fp_dhdl
, " %#.8g", md
->dE
[i
]);
1255 (md
->epc
!= epcNO
) &&
1256 (enerd
->n_lambda
> 0) &&
1257 (fep
->init_lambda
< 0))
1259 fprintf(md
->fp_dhdl
, " %#.8g", pv
); /* PV term only needed when
1260 there are alternate state
1261 lambda and we're not in
1262 compatibility mode */
1264 fprintf(md
->fp_dhdl
, "\n");
1265 /* and the binary free energy output */
1267 if (md
->dhc
&& bDoDHDL
)
1270 for (i
= 0; i
< efptNR
; i
++)
1272 if (fep
->separate_dvdl
[i
])
1274 /* assumes F_DVDL is first */
1275 store_dhdl
[idhdl
] = enerd
->term
[F_DVDL
+i
];
1279 store_energy
= enerd
->term
[F_ETOT
];
1280 /* store_dh is dE */
1281 mde_delta_h_coll_add_dh(md
->dhc
,
1282 (double)state
->fep_state
,
1286 md
->dE
+ fep
->lambda_start_n
,
1293 void upd_mdebin_step(t_mdebin
*md
)
1295 ebin_increase_count(md
->ebin
, FALSE
);
1298 static void npr(FILE *log
, int n
, char c
)
1300 for (; (n
> 0); n
--)
1302 fprintf(log
, "%c", c
);
1306 static void pprint(FILE *log
, const char *s
, t_mdebin
*md
)
1310 char buf1
[22], buf2
[22];
1313 fprintf(log
, "\t<====== ");
1314 npr(log
, slen
, CHAR
);
1315 fprintf(log
, " ==>\n");
1316 fprintf(log
, "\t<==== %s ====>\n", s
);
1317 fprintf(log
, "\t<== ");
1318 npr(log
, slen
, CHAR
);
1319 fprintf(log
, " ======>\n\n");
1321 fprintf(log
, "\tStatistics over %s steps using %s frames\n",
1322 gmx_step_str(md
->ebin
->nsteps_sim
, buf1
),
1323 gmx_step_str(md
->ebin
->nsum_sim
, buf2
));
1327 void print_ebin_header(FILE *log
, gmx_int64_t steps
, double time
, real lambda
)
1331 fprintf(log
, " %12s %12s %12s\n"
1332 " %12s %12.5f %12.5f\n\n",
1333 "Step", "Time", "Lambda", gmx_step_str(steps
, buf
), time
, lambda
);
1336 void print_ebin(ener_file_t fp_ene
, gmx_bool bEne
, gmx_bool bDR
, gmx_bool bOR
,
1338 gmx_int64_t step
, double time
,
1339 int mode
, gmx_bool bCompact
,
1340 t_mdebin
*md
, t_fcdata
*fcd
,
1341 gmx_groups_t
*groups
, t_grpopts
*opts
)
1343 /*static char **grpnms=NULL;*/
1345 int i
, j
, n
, ni
, nj
, b
;
1347 real
*disre_rm3tav
, *disre_rt
;
1349 /* these are for the old-style blocks (1 subblock, only reals), because
1350 there can be only one per ID for these */
1363 fr
.nsteps
= md
->ebin
->nsteps
;
1364 fr
.dt
= md
->delta_t
;
1365 fr
.nsum
= md
->ebin
->nsum
;
1366 fr
.nre
= (bEne
) ? md
->ebin
->nener
: 0;
1367 fr
.ener
= md
->ebin
->e
;
1368 ndisre
= bDR
? fcd
->disres
.npair
: 0;
1369 disre_rm3tav
= fcd
->disres
.rm3tav
;
1370 disre_rt
= fcd
->disres
.rt
;
1371 /* Optional additional old-style (real-only) blocks. */
1372 for (i
= 0; i
< enxNR
; i
++)
1376 if (fcd
->orires
.nr
> 0 && bOR
)
1378 diagonalize_orires_tensors(&(fcd
->orires
));
1379 nr
[enxOR
] = fcd
->orires
.nr
;
1380 block
[enxOR
] = fcd
->orires
.otav
;
1382 nr
[enxORI
] = (fcd
->orires
.oinsl
!= fcd
->orires
.otav
) ?
1384 block
[enxORI
] = fcd
->orires
.oinsl
;
1385 id
[enxORI
] = enxORI
;
1386 nr
[enxORT
] = fcd
->orires
.nex
*12;
1387 block
[enxORT
] = fcd
->orires
.eig
;
1388 id
[enxORT
] = enxORT
;
1391 /* whether we are going to wrte anything out: */
1392 if (fr
.nre
|| ndisre
|| nr
[enxOR
] || nr
[enxORI
])
1395 /* the old-style blocks go first */
1397 for (i
= 0; i
< enxNR
; i
++)
1404 add_blocks_enxframe(&fr
, fr
.nblock
);
1405 for (b
= 0; b
< fr
.nblock
; b
++)
1407 add_subblocks_enxblock(&(fr
.block
[b
]), 1);
1408 fr
.block
[b
].id
= id
[b
];
1409 fr
.block
[b
].sub
[0].nr
= nr
[b
];
1411 fr
.block
[b
].sub
[0].type
= xdr_datatype_float
;
1412 fr
.block
[b
].sub
[0].fval
= block
[b
];
1414 fr
.block
[b
].sub
[0].type
= xdr_datatype_double
;
1415 fr
.block
[b
].sub
[0].dval
= block
[b
];
1419 /* check for disre block & fill it. */
1424 add_blocks_enxframe(&fr
, fr
.nblock
);
1426 add_subblocks_enxblock(&(fr
.block
[db
]), 2);
1427 fr
.block
[db
].id
= enxDISRE
;
1428 fr
.block
[db
].sub
[0].nr
= ndisre
;
1429 fr
.block
[db
].sub
[1].nr
= ndisre
;
1431 fr
.block
[db
].sub
[0].type
= xdr_datatype_float
;
1432 fr
.block
[db
].sub
[1].type
= xdr_datatype_float
;
1433 fr
.block
[db
].sub
[0].fval
= disre_rt
;
1434 fr
.block
[db
].sub
[1].fval
= disre_rm3tav
;
1436 fr
.block
[db
].sub
[0].type
= xdr_datatype_double
;
1437 fr
.block
[db
].sub
[1].type
= xdr_datatype_double
;
1438 fr
.block
[db
].sub
[0].dval
= disre_rt
;
1439 fr
.block
[db
].sub
[1].dval
= disre_rm3tav
;
1442 /* here we can put new-style blocks */
1444 /* Free energy perturbation blocks */
1447 mde_delta_h_coll_handle_block(md
->dhc
, &fr
, fr
.nblock
);
1450 /* we can now free & reset the data in the blocks */
1453 mde_delta_h_coll_reset(md
->dhc
);
1456 /* do the actual I/O */
1457 do_enx(fp_ene
, &fr
);
1460 /* We have stored the sums, so reset the sum history */
1461 reset_ebin_sums(md
->ebin
);
1469 pprint(log
, "A V E R A G E S", md
);
1475 pprint(log
, "R M S - F L U C T U A T I O N S", md
);
1479 gmx_fatal(FARGS
, "Invalid print mode (%d)", mode
);
1484 for (i
= 0; i
< opts
->ngtc
; i
++)
1486 if (opts
->annealing
[i
] != eannNO
)
1488 fprintf(log
, "Current ref_t for group %s: %8.1f\n",
1489 *(groups
->grpname
[groups
->grps
[egcTC
].nm_ind
[i
]]),
1493 if (mode
== eprNORMAL
&& fcd
->orires
.nr
> 0)
1495 print_orires_log(log
, &(fcd
->orires
));
1497 fprintf(log
, " Energies (%s)\n", unit_energy
);
1498 pr_ebin(log
, md
->ebin
, md
->ie
, md
->f_nre
+md
->nCrmsd
, 5, mode
, TRUE
);
1505 pr_ebin(log
, md
->ebin
, md
->ib
, md
->bTricl
? NTRICLBOXS
: NBOXS
, 5,
1511 fprintf(log
, " Constraint Virial (%s)\n", unit_energy
);
1512 pr_ebin(log
, md
->ebin
, md
->isvir
, 9, 3, mode
, FALSE
);
1514 fprintf(log
, " Force Virial (%s)\n", unit_energy
);
1515 pr_ebin(log
, md
->ebin
, md
->ifvir
, 9, 3, mode
, FALSE
);
1520 fprintf(log
, " Total Virial (%s)\n", unit_energy
);
1521 pr_ebin(log
, md
->ebin
, md
->ivir
, 9, 3, mode
, FALSE
);
1526 fprintf(log
, " Pressure (%s)\n", unit_pres_bar
);
1527 pr_ebin(log
, md
->ebin
, md
->ipres
, 9, 3, mode
, FALSE
);
1532 fprintf(log
, " Total Dipole (%s)\n", unit_dipole_D
);
1533 pr_ebin(log
, md
->ebin
, md
->imu
, 3, 3, mode
, FALSE
);
1539 if (md
->print_grpnms
== NULL
)
1541 snew(md
->print_grpnms
, md
->nE
);
1543 for (i
= 0; (i
< md
->nEg
); i
++)
1545 ni
= groups
->grps
[egcENER
].nm_ind
[i
];
1546 for (j
= i
; (j
< md
->nEg
); j
++)
1548 nj
= groups
->grps
[egcENER
].nm_ind
[j
];
1549 sprintf(buf
, "%s-%s", *(groups
->grpname
[ni
]),
1550 *(groups
->grpname
[nj
]));
1551 md
->print_grpnms
[n
++] = gmx_strdup(buf
);
1555 sprintf(buf
, "Epot (%s)", unit_energy
);
1556 fprintf(log
, "%15s ", buf
);
1557 for (i
= 0; (i
< egNR
); i
++)
1561 fprintf(log
, "%12s ", egrp_nm
[i
]);
1565 for (i
= 0; (i
< md
->nE
); i
++)
1567 fprintf(log
, "%15s", md
->print_grpnms
[i
]);
1568 pr_ebin(log
, md
->ebin
, md
->igrp
[i
], md
->nEc
, md
->nEc
, mode
,
1575 pr_ebin(log
, md
->ebin
, md
->itemp
, md
->nTC
, 4, mode
, TRUE
);
1580 fprintf(log
, "%15s %12s %12s %12s\n",
1581 "Group", "Ux", "Uy", "Uz");
1582 for (i
= 0; (i
< md
->nU
); i
++)
1584 ni
= groups
->grps
[egcACC
].nm_ind
[i
];
1585 fprintf(log
, "%15s", *groups
->grpname
[ni
]);
1586 pr_ebin(log
, md
->ebin
, md
->iu
+3*i
, 3, 3, mode
, FALSE
);
1595 void update_energyhistory(energyhistory_t
* enerhist
, t_mdebin
* mdebin
)
1599 enerhist
->nsteps
= mdebin
->ebin
->nsteps
;
1600 enerhist
->nsum
= mdebin
->ebin
->nsum
;
1601 enerhist
->nsteps_sim
= mdebin
->ebin
->nsteps_sim
;
1602 enerhist
->nsum_sim
= mdebin
->ebin
->nsum_sim
;
1603 enerhist
->nener
= mdebin
->ebin
->nener
;
1605 if (mdebin
->ebin
->nsum
> 0)
1607 /* Check if we need to allocate first */
1608 if (enerhist
->ener_ave
== NULL
)
1610 snew(enerhist
->ener_ave
, enerhist
->nener
);
1611 snew(enerhist
->ener_sum
, enerhist
->nener
);
1614 for (i
= 0; i
< enerhist
->nener
; i
++)
1616 enerhist
->ener_ave
[i
] = mdebin
->ebin
->e
[i
].eav
;
1617 enerhist
->ener_sum
[i
] = mdebin
->ebin
->e
[i
].esum
;
1621 if (mdebin
->ebin
->nsum_sim
> 0)
1623 /* Check if we need to allocate first */
1624 if (enerhist
->ener_sum_sim
== NULL
)
1626 snew(enerhist
->ener_sum_sim
, enerhist
->nener
);
1629 for (i
= 0; i
< enerhist
->nener
; i
++)
1631 enerhist
->ener_sum_sim
[i
] = mdebin
->ebin
->e_sim
[i
].esum
;
1636 mde_delta_h_coll_update_energyhistory(mdebin
->dhc
, enerhist
);
1640 void restore_energyhistory_from_state(t_mdebin
* mdebin
,
1641 energyhistory_t
* enerhist
)
1645 if ((enerhist
->nsum
> 0 || enerhist
->nsum_sim
> 0) &&
1646 mdebin
->ebin
->nener
!= enerhist
->nener
)
1648 gmx_fatal(FARGS
, "Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1649 mdebin
->ebin
->nener
, enerhist
->nener
);
1652 mdebin
->ebin
->nsteps
= enerhist
->nsteps
;
1653 mdebin
->ebin
->nsum
= enerhist
->nsum
;
1654 mdebin
->ebin
->nsteps_sim
= enerhist
->nsteps_sim
;
1655 mdebin
->ebin
->nsum_sim
= enerhist
->nsum_sim
;
1657 for (i
= 0; i
< mdebin
->ebin
->nener
; i
++)
1659 mdebin
->ebin
->e
[i
].eav
=
1660 (enerhist
->nsum
> 0 ? enerhist
->ener_ave
[i
] : 0);
1661 mdebin
->ebin
->e
[i
].esum
=
1662 (enerhist
->nsum
> 0 ? enerhist
->ener_sum
[i
] : 0);
1663 mdebin
->ebin
->e_sim
[i
].esum
=
1664 (enerhist
->nsum_sim
> 0 ? enerhist
->ener_sum_sim
[i
] : 0);
1668 mde_delta_h_coll_restore_energyhistory(mdebin
->dhc
, enerhist
);