1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
55 #include "mtop_util.h"
59 #include "mdebin_bar.h"
62 static const char *conrmsd_nm
[] = { "Constr. rmsd", "Constr.2 rmsd" };
64 static const char *boxs_nm
[] = { "Box-X", "Box-Y", "Box-Z" };
66 static const char *tricl_boxs_nm
[] = {
67 "Box-XX", "Box-YY", "Box-ZZ",
68 "Box-YX", "Box-ZX", "Box-ZY"
71 static const char *vol_nm
[] = { "Volume" };
73 static const char *dens_nm
[] = {"Density" };
75 static const char *pv_nm
[] = {"pV" };
77 static const char *enthalpy_nm
[] = {"Enthalpy" };
79 static const char *boxvel_nm
[] = {
80 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
81 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
84 #define NBOXS asize(boxs_nm)
85 #define NTRICLBOXS asize(tricl_boxs_nm)
87 t_mdebin
*init_mdebin(ener_file_t fp_ene
,
88 const gmx_mtop_t
*mtop
,
92 const char *ener_nm
[F_NRE
];
93 static const char *vir_nm
[] = {
94 "Vir-XX", "Vir-XY", "Vir-XZ",
95 "Vir-YX", "Vir-YY", "Vir-YZ",
96 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
98 static const char *sv_nm
[] = {
99 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
100 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
101 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
103 static const char *fv_nm
[] = {
104 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
105 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
106 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
108 static const char *pres_nm
[] = {
109 "Pres-XX","Pres-XY","Pres-XZ",
110 "Pres-YX","Pres-YY","Pres-YZ",
111 "Pres-ZX","Pres-ZY","Pres-ZZ"
113 static const char *surft_nm
[] = {
116 static const char *mu_nm
[] = {
117 "Mu-X", "Mu-Y", "Mu-Z"
119 static const char *vcos_nm
[] = {
122 static const char *visc_nm
[] = {
125 static const char *baro_nm
[] = {
130 const gmx_groups_t
*groups
;
135 int i
,j
,ni
,nj
,n
,nh
,k
,kk
,ncon
,nset
;
136 gmx_bool bBHAM
,bNoseHoover
,b14
;
145 if (EI_DYNAMICS(ir
->eI
))
147 md
->delta_t
= ir
->delta_t
;
154 groups
= &mtop
->groups
;
156 bBHAM
= (mtop
->ffparams
.functype
[0] == F_BHAM
);
157 b14
= (gmx_mtop_ftype_count(mtop
,F_LJ14
) > 0 ||
158 gmx_mtop_ftype_count(mtop
,F_LJC14_Q
) > 0);
160 ncon
= gmx_mtop_ftype_count(mtop
,F_CONSTR
);
161 nset
= gmx_mtop_ftype_count(mtop
,F_SETTLE
);
162 md
->bConstr
= (ncon
> 0 || nset
> 0);
163 md
->bConstrVir
= FALSE
;
165 if (ncon
> 0 && ir
->eConstrAlg
== econtLINCS
) {
171 md
->bConstrVir
= (getenv("GMX_CONSTRAINTVIR") != NULL
);
176 /* Energy monitoring */
183 for(i
=0; i
<F_NRE
; i
++)
185 md
->bEner
[i
] = FALSE
;
187 md
->bEner
[i
] = !bBHAM
;
188 else if (i
== F_BHAM
)
189 md
->bEner
[i
] = bBHAM
;
191 md
->bEner
[i
] = ir
->bQMMM
;
192 else if (i
== F_COUL_LR
)
193 md
->bEner
[i
] = (ir
->rcoulomb
> ir
->rlist
);
194 else if (i
== F_LJ_LR
)
195 md
->bEner
[i
] = (!bBHAM
&& ir
->rvdw
> ir
->rlist
);
196 else if (i
== F_BHAM_LR
)
197 md
->bEner
[i
] = (bBHAM
&& ir
->rvdw
> ir
->rlist
);
198 else if (i
== F_RF_EXCL
)
199 md
->bEner
[i
] = (EEL_RF(ir
->coulombtype
) && ir
->coulombtype
!= eelRF_NEC
&& ir
->cutoff_scheme
== ecutsGROUP
);
200 else if (i
== F_COUL_RECIP
)
201 md
->bEner
[i
] = EEL_FULL(ir
->coulombtype
);
202 else if (i
== F_LJ14
)
204 else if (i
== F_COUL14
)
206 else if (i
== F_LJC14_Q
|| i
== F_LJC_PAIRS_NB
)
207 md
->bEner
[i
] = FALSE
;
208 else if ((i
== F_DVDL_COUL
&& ir
->fepvals
->separate_dvdl
[efptCOUL
]) ||
209 (i
== F_DVDL_VDW
&& ir
->fepvals
->separate_dvdl
[efptVDW
]) ||
210 (i
== F_DVDL_BONDED
&& ir
->fepvals
->separate_dvdl
[efptBONDED
]) ||
211 (i
== F_DVDL_RESTRAINT
&& ir
->fepvals
->separate_dvdl
[efptRESTRAINT
]) ||
212 (i
== F_DKDL
&& ir
->fepvals
->separate_dvdl
[efptMASS
]) ||
213 (i
== F_DVDL
&& ir
->fepvals
->separate_dvdl
[efptFEP
]))
214 md
->bEner
[i
] = (ir
->efep
!= efepNO
);
215 else if ((interaction_function
[i
].flags
& IF_VSITE
) ||
216 (i
== F_CONSTR
) || (i
== F_CONSTRNC
) || (i
== F_SETTLE
))
217 md
->bEner
[i
] = FALSE
;
218 else if ((i
== F_COUL_SR
) || (i
== F_EPOT
) || (i
== F_PRES
) || (i
==F_EQM
))
220 else if ((i
== F_GBPOL
) && ir
->implicit_solvent
==eisGBSA
)
222 else if ((i
== F_NPSOLVATION
) && ir
->implicit_solvent
==eisGBSA
&& (ir
->sa_algorithm
!= esaNO
))
224 else if ((i
== F_GB12
) || (i
== F_GB13
) || (i
== F_GB14
))
225 md
->bEner
[i
] = FALSE
;
226 else if ((i
== F_ETOT
) || (i
== F_EKIN
) || (i
== F_TEMP
))
227 md
->bEner
[i
] = EI_DYNAMICS(ir
->eI
);
229 md
->bEner
[i
] = (EI_DYNAMICS(ir
->eI
) && getenv("GMX_VIRIAL_TEMPERATURE"));
230 else if (i
== F_DISPCORR
|| i
== F_PDISPCORR
)
231 md
->bEner
[i
] = (ir
->eDispCorr
!= edispcNO
);
232 else if (i
== F_DISRESVIOL
)
233 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
,F_DISRES
) > 0);
234 else if (i
== F_ORIRESDEV
)
235 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
,F_ORIRES
) > 0);
236 else if (i
== F_CONNBONDS
)
237 md
->bEner
[i
] = FALSE
;
238 else if (i
== F_COM_PULL
)
239 md
->bEner
[i
] = (ir
->ePull
== epullUMBRELLA
|| ir
->ePull
== epullCONST_F
|| ir
->bRot
);
240 else if (i
== F_ECONSERVED
)
241 md
->bEner
[i
] = ((ir
->etc
== etcNOSEHOOVER
|| ir
->etc
== etcVRESCALE
) &&
242 (ir
->epc
== epcNO
|| ir
->epc
==epcMTTK
));
244 md
->bEner
[i
] = (gmx_mtop_ftype_count(mtop
,i
) > 0);
247 /* OpenMM always produces only the following 4 energy terms */
248 md
->bEner
[F_EPOT
] = TRUE
;
249 md
->bEner
[F_EKIN
] = TRUE
;
250 md
->bEner
[F_ETOT
] = TRUE
;
251 md
->bEner
[F_TEMP
] = TRUE
;
254 /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
255 if (ir
->bAdress
&& !debug
) {
256 for (i
= 0; i
< F_NRE
; i
++) {
257 md
->bEner
[i
] = FALSE
;
258 if(i
== F_EKIN
){ md
->bEner
[i
] = TRUE
;}
259 if(i
== F_TEMP
){ md
->bEner
[i
] = TRUE
;}
268 for(i
=0; i
<F_NRE
; i
++)
272 ener_nm
[md
->f_nre
]=interaction_function
[i
].longname
;
278 md
->bDiagPres
= !TRICLINIC(ir
->ref_p
);
279 md
->ref_p
= (ir
->ref_p
[XX
][XX
]+ir
->ref_p
[YY
][YY
]+ir
->ref_p
[ZZ
][ZZ
])/DIM
;
280 md
->bTricl
= TRICLINIC(ir
->compress
) || TRICLINIC(ir
->deform
);
281 md
->bDynBox
= DYNAMIC_BOX(*ir
);
283 md
->bNHC_trotter
= IR_NVT_TROTTER(ir
);
284 md
->bPrintNHChains
= ir
-> bPrintNHChains
;
285 md
->bMTTK
= (IR_NPT_TROTTER(ir
) || IR_NPH_TROTTER(ir
));
286 md
->bMu
= NEED_MUTOT(*ir
);
288 md
->ebin
= mk_ebin();
289 /* Pass NULL for unit to let get_ebin_space determine the units
290 * for interaction_function[i].longname
292 md
->ie
= get_ebin_space(md
->ebin
,md
->f_nre
,ener_nm
,NULL
);
295 /* This should be called directly after the call for md->ie,
296 * such that md->iconrmsd follows directly in the list.
298 md
->iconrmsd
= get_ebin_space(md
->ebin
,md
->nCrmsd
,conrmsd_nm
,"");
302 md
->ib
= get_ebin_space(md
->ebin
,
303 md
->bTricl
? NTRICLBOXS
: NBOXS
,
304 md
->bTricl
? tricl_boxs_nm
: boxs_nm
,
306 md
->ivol
= get_ebin_space(md
->ebin
, 1, vol_nm
, unit_volume
);
307 md
->idens
= get_ebin_space(md
->ebin
, 1, dens_nm
, unit_density_SI
);
310 md
->ipv
= get_ebin_space(md
->ebin
, 1, pv_nm
, unit_energy
);
311 md
->ienthalpy
= get_ebin_space(md
->ebin
, 1, enthalpy_nm
, unit_energy
);
316 md
->isvir
= get_ebin_space(md
->ebin
,asize(sv_nm
),sv_nm
,unit_energy
);
317 md
->ifvir
= get_ebin_space(md
->ebin
,asize(fv_nm
),fv_nm
,unit_energy
);
320 md
->ivir
= get_ebin_space(md
->ebin
,asize(vir_nm
),vir_nm
,unit_energy
);
322 md
->ipres
= get_ebin_space(md
->ebin
,asize(pres_nm
),pres_nm
,unit_pres_bar
);
324 md
->isurft
= get_ebin_space(md
->ebin
,asize(surft_nm
),surft_nm
,
326 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
328 md
->ipc
= get_ebin_space(md
->ebin
,md
->bTricl
? 6 : 3,
333 md
->imu
= get_ebin_space(md
->ebin
,asize(mu_nm
),mu_nm
,unit_dipole_D
);
335 if (ir
->cos_accel
!= 0)
337 md
->ivcos
= get_ebin_space(md
->ebin
,asize(vcos_nm
),vcos_nm
,unit_vel
);
338 md
->ivisc
= get_ebin_space(md
->ebin
,asize(visc_nm
),visc_nm
,
342 /* Energy monitoring */
345 md
->bEInd
[i
] = FALSE
;
347 md
->bEInd
[egCOULSR
] = TRUE
;
348 md
->bEInd
[egLJSR
] = TRUE
;
350 if (ir
->rcoulomb
> ir
->rlist
)
352 md
->bEInd
[egCOULLR
] = TRUE
;
356 if (ir
->rvdw
> ir
->rlist
)
358 md
->bEInd
[egLJLR
] = TRUE
;
363 md
->bEInd
[egLJSR
] = FALSE
;
364 md
->bEInd
[egBHAMSR
] = TRUE
;
365 if (ir
->rvdw
> ir
->rlist
)
367 md
->bEInd
[egBHAMLR
] = TRUE
;
372 md
->bEInd
[egLJ14
] = TRUE
;
373 md
->bEInd
[egCOUL14
] = TRUE
;
376 for(i
=0; (i
<egNR
); i
++)
384 n
=groups
->grps
[egcENER
].nr
;
385 /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
387 /*standard simulation*/
392 /*AdResS simulation*/
399 snew(md
->igrp
,md
->nE
);
404 for(k
=0; (k
<md
->nEc
); k
++)
408 for(i
=0; (i
<groups
->grps
[egcENER
].nr
); i
++)
410 ni
=groups
->grps
[egcENER
].nm_ind
[i
];
411 for(j
=i
; (j
<groups
->grps
[egcENER
].nr
); j
++)
413 nj
=groups
->grps
[egcENER
].nm_ind
[j
];
414 for(k
=kk
=0; (k
<egNR
); k
++)
418 sprintf(gnm
[kk
],"%s:%s-%s",egrp_nm
[k
],
419 *(groups
->grpname
[ni
]),*(groups
->grpname
[nj
]));
423 md
->igrp
[n
]=get_ebin_space(md
->ebin
,md
->nEc
,
424 (const char **)gnm
,unit_energy
);
428 for(k
=0; (k
<md
->nEc
); k
++)
436 gmx_incons("Number of energy terms wrong");
440 md
->nTC
=groups
->grps
[egcTC
].nr
;
441 md
->nNHC
= ir
->opts
.nhchainlength
; /* shorthand for number of NH chains */
444 md
->nTCP
= 1; /* assume only one possible coupling system for barostat
451 if (md
->etc
== etcNOSEHOOVER
)
453 if (md
->bNHC_trotter
)
455 md
->mde_n
= 2*md
->nNHC
*md
->nTC
;
459 md
->mde_n
= 2*md
->nTC
;
461 if (md
->epc
== epcMTTK
)
463 md
->mdeb_n
= 2*md
->nNHC
*md
->nTCP
;
470 snew(md
->tmp_r
,md
->mde_n
);
471 snew(md
->tmp_v
,md
->mde_n
);
472 snew(md
->grpnms
,md
->mde_n
);
475 for(i
=0; (i
<md
->nTC
); i
++)
477 ni
=groups
->grps
[egcTC
].nm_ind
[i
];
478 sprintf(buf
,"T-%s",*(groups
->grpname
[ni
]));
479 grpnms
[i
]=strdup(buf
);
481 md
->itemp
=get_ebin_space(md
->ebin
,md
->nTC
,(const char **)grpnms
,
484 if (md
->etc
== etcNOSEHOOVER
)
486 if (md
->bPrintNHChains
)
488 if (md
->bNHC_trotter
)
490 for(i
=0; (i
<md
->nTC
); i
++)
492 ni
=groups
->grps
[egcTC
].nm_ind
[i
];
493 bufi
= *(groups
->grpname
[ni
]);
494 for(j
=0; (j
<md
->nNHC
); j
++)
496 sprintf(buf
,"Xi-%d-%s",j
,bufi
);
497 grpnms
[2*(i
*md
->nNHC
+j
)]=strdup(buf
);
498 sprintf(buf
,"vXi-%d-%s",j
,bufi
);
499 grpnms
[2*(i
*md
->nNHC
+j
)+1]=strdup(buf
);
502 md
->itc
=get_ebin_space(md
->ebin
,md
->mde_n
,
503 (const char **)grpnms
,unit_invtime
);
506 for(i
=0; (i
<md
->nTCP
); i
++)
508 bufi
= baro_nm
[0]; /* All barostat DOF's together for now. */
509 for(j
=0; (j
<md
->nNHC
); j
++)
511 sprintf(buf
,"Xi-%d-%s",j
,bufi
);
512 grpnms
[2*(i
*md
->nNHC
+j
)]=strdup(buf
);
513 sprintf(buf
,"vXi-%d-%s",j
,bufi
);
514 grpnms
[2*(i
*md
->nNHC
+j
)+1]=strdup(buf
);
517 md
->itcb
=get_ebin_space(md
->ebin
,md
->mdeb_n
,
518 (const char **)grpnms
,unit_invtime
);
523 for(i
=0; (i
<md
->nTC
); i
++)
525 ni
=groups
->grps
[egcTC
].nm_ind
[i
];
526 bufi
= *(groups
->grpname
[ni
]);
527 sprintf(buf
,"Xi-%s",bufi
);
528 grpnms
[2*i
]=strdup(buf
);
529 sprintf(buf
,"vXi-%s",bufi
);
530 grpnms
[2*i
+1]=strdup(buf
);
532 md
->itc
=get_ebin_space(md
->ebin
,md
->mde_n
,
533 (const char **)grpnms
,unit_invtime
);
537 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
538 md
->etc
== etcVRESCALE
)
540 for(i
=0; (i
<md
->nTC
); i
++)
542 ni
=groups
->grps
[egcTC
].nm_ind
[i
];
543 sprintf(buf
,"Lamb-%s",*(groups
->grpname
[ni
]));
544 grpnms
[i
]=strdup(buf
);
546 md
->itc
=get_ebin_space(md
->ebin
,md
->mde_n
,(const char **)grpnms
,"");
552 md
->nU
=groups
->grps
[egcACC
].nr
;
555 snew(grpnms
,3*md
->nU
);
556 for(i
=0; (i
<md
->nU
); i
++)
558 ni
=groups
->grps
[egcACC
].nm_ind
[i
];
559 sprintf(buf
,"Ux-%s",*(groups
->grpname
[ni
]));
560 grpnms
[3*i
+XX
]=strdup(buf
);
561 sprintf(buf
,"Uy-%s",*(groups
->grpname
[ni
]));
562 grpnms
[3*i
+YY
]=strdup(buf
);
563 sprintf(buf
,"Uz-%s",*(groups
->grpname
[ni
]));
564 grpnms
[3*i
+ZZ
]=strdup(buf
);
566 md
->iu
=get_ebin_space(md
->ebin
,3*md
->nU
,(const char **)grpnms
,unit_vel
);
572 do_enxnms(fp_ene
,&md
->ebin
->nener
,&md
->ebin
->enm
);
575 md
->print_grpnms
=NULL
;
577 /* check whether we're going to write dh histograms */
579 if (ir
->fepvals
->separate_dhdl_file
== esepdhdlfileNO
)
581 /* Currently dh histograms are only written with dynamics */
582 if (EI_DYNAMICS(ir
->eI
))
586 mde_delta_h_coll_init(md
->dhc
, ir
);
592 md
->fp_dhdl
= fp_dhdl
;
596 snew(md
->temperatures
,ir
->fepvals
->n_lambda
);
597 for (i
=0;i
<ir
->fepvals
->n_lambda
;i
++)
599 md
->temperatures
[i
] = ir
->simtempvals
->temperatures
[i
];
605 extern FILE *open_dhdl(const char *filename
,const t_inputrec
*ir
,
606 const output_env_t oenv
)
609 const char *dhdl
="dH/d\\lambda",*deltag
="\\DeltaH",*lambda
="\\lambda",
610 *lambdastate
="\\lambda state",*remain
="remaining";
611 char title
[STRLEN
],label_x
[STRLEN
],label_y
[STRLEN
];
612 int i
,np
,nps
,nsets
,nsets_de
,nsetsbegin
;
625 if (fep
->n_lambda
== 0)
627 sprintf(title
,"%s",dhdl
);
628 sprintf(label_x
,"Time (ps)");
629 sprintf(label_y
,"%s (%s %s)",
630 dhdl
,unit_energy
,"[\\lambda]\\S-1\\N");
634 sprintf(title
,"%s and %s",dhdl
,deltag
);
635 sprintf(label_x
,"Time (ps)");
636 sprintf(label_y
,"%s and %s (%s %s)",
637 dhdl
,deltag
,unit_energy
,"[\\8l\\4]\\S-1\\N");
639 fp
= gmx_fio_fopen(filename
,"w+");
640 xvgr_header(fp
,title
,label_x
,label_y
,exvggtXNY
,oenv
);
644 bufplace
= sprintf(buf
,"T = %g (K) ",
647 if (ir
->efep
!= efepSLOWGROWTH
)
649 if (fep
->n_lambda
== 0)
651 sprintf(&(buf
[bufplace
]),"%s = %g",
652 lambda
,fep
->init_lambda
);
656 sprintf(&(buf
[bufplace
]),"%s = %d",
657 lambdastate
,fep
->init_fep_state
);
660 xvgr_subtitle(fp
,buf
,oenv
);
662 for (i
=0;i
<efptNR
;i
++)
664 if (fep
->separate_dvdl
[i
]) {nsets_dhdl
++;}
667 /* count the number of delta_g states */
668 nsets_de
= fep
->n_lambda
;
670 nsets
= nsets_dhdl
+ nsets_de
; /* dhdl + fep differences */
672 if (fep
->n_lambda
>0 && ir
->bExpanded
)
674 nsets
+= 1; /*add fep state for expanded ensemble */
677 if (fep
->bPrintEnergy
)
679 nsets
+= 1; /* add energy to the dhdl as well */
683 if ((ir
->epc
!=epcNO
) && (fep
->n_lambda
>0))
685 nsetsextend
+= 1; /* for PV term, other terms possible if required for the reduced potential (only needed with foreign lambda) */
687 snew(setname
,nsetsextend
);
691 /* state for the fep_vals, if we have alchemical sampling */
692 sprintf(buf
,"%s","Thermodynamic state");
693 setname
[s
] = strdup(buf
);
697 if (fep
->bPrintEnergy
)
699 sprintf(buf
,"%s (%s)","Energy",unit_energy
);
700 setname
[s
] = strdup(buf
);
704 for (i
=0;i
<efptNR
;i
++)
706 if (fep
->separate_dvdl
[i
]) {
707 sprintf(buf
,"%s (%s)",dhdl
,efpt_names
[i
]);
708 setname
[s
] = strdup(buf
);
713 if (fep
->n_lambda
> 0)
715 /* g_bar has to determine the lambda values used in this simulation
716 * from this xvg legend.
720 nsetsbegin
= 1; /* for including the expanded ensemble */
725 if (fep
->bPrintEnergy
)
729 nsetsbegin
+= nsets_dhdl
;
731 for(s
=nsetsbegin
; s
<nsets
; s
++)
733 nps
= sprintf(buf
,"%s %s (",deltag
,lambda
);
734 for (i
=0;i
<efptNR
;i
++)
736 if (fep
->separate_dvdl
[i
])
738 np
= sprintf(&buf
[nps
],"%g,",fep
->all_lambda
[i
][s
-(nsetsbegin
)]);
744 /* print the temperature for this state if doing simulated annealing */
745 sprintf(&buf
[nps
],"T = %g (%s))",ir
->simtempvals
->temperatures
[s
-(nsetsbegin
)],unit_temp_K
);
749 sprintf(&buf
[nps
-1],")"); /* -1 to overwrite the last comma */
751 setname
[s
] = strdup(buf
);
753 if (ir
->epc
!=epcNO
) {
754 np
= sprintf(buf
,"pV (%s)",unit_energy
);
755 setname
[nsetsextend
-1] = strdup(buf
); /* the first entry after nsets */
758 xvgr_legend(fp
,nsetsextend
,(const char **)setname
,oenv
);
760 for(s
=0; s
<nsetsextend
; s
++)
770 static void copy_energy(t_mdebin
*md
, real e
[],real ecpy
[])
774 for(i
=j
=0; (i
<F_NRE
); i
++)
778 gmx_incons("Number of energy terms wrong");
781 void upd_mdebin(t_mdebin
*md
,
786 gmx_enerdata_t
*enerd
,
795 gmx_ekindata_t
*ekind
,
799 int i
,j
,k
,kk
,m
,n
,gid
;
800 real crmsd
[2],tmp6
[6];
801 real bs
[NTRICLBOXS
],vol
,dens
,pv
,enthalpy
;
804 double store_dhdl
[efptNR
];
809 /* Do NOT use the box in the state variable, but the separate box provided
810 * as an argument. This is because we sometimes need to write the box from
811 * the last timestep to match the trajectory frames.
813 copy_energy(md
, enerd
->term
,ecopy
);
814 add_ebin(md
->ebin
,md
->ie
,md
->f_nre
,ecopy
,bSum
);
817 crmsd
[0] = constr_rmsd(constr
,FALSE
);
820 crmsd
[1] = constr_rmsd(constr
,TRUE
);
822 add_ebin(md
->ebin
,md
->iconrmsd
,md
->nCrmsd
,crmsd
,FALSE
);
844 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
845 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
846 add_ebin(md
->ebin
,md
->ib
,nboxs
,bs
,bSum
);
847 add_ebin(md
->ebin
,md
->ivol
,1 ,&vol
,bSum
);
848 add_ebin(md
->ebin
,md
->idens
,1 ,&dens
,bSum
);
852 /* This is pV (in kJ/mol). The pressure is the reference pressure,
853 not the instantaneous pressure */
854 pv
= vol
*md
->ref_p
/PRESFAC
;
856 add_ebin(md
->ebin
,md
->ipv
,1 ,&pv
,bSum
);
857 enthalpy
= pv
+ enerd
->term
[F_ETOT
];
858 add_ebin(md
->ebin
,md
->ienthalpy
,1 ,&enthalpy
,bSum
);
863 add_ebin(md
->ebin
,md
->isvir
,9,svir
[0],bSum
);
864 add_ebin(md
->ebin
,md
->ifvir
,9,fvir
[0],bSum
);
867 add_ebin(md
->ebin
,md
->ivir
,9,vir
[0],bSum
);
869 add_ebin(md
->ebin
,md
->ipres
,9,pres
[0],bSum
);
871 tmp
= (pres
[ZZ
][ZZ
]-(pres
[XX
][XX
]+pres
[YY
][YY
])*0.5)*box
[ZZ
][ZZ
];
872 add_ebin(md
->ebin
,md
->isurft
,1,&tmp
,bSum
);
874 if (md
->epc
== epcPARRINELLORAHMAN
|| md
->epc
== epcMTTK
)
876 tmp6
[0] = state
->boxv
[XX
][XX
];
877 tmp6
[1] = state
->boxv
[YY
][YY
];
878 tmp6
[2] = state
->boxv
[ZZ
][ZZ
];
879 tmp6
[3] = state
->boxv
[YY
][XX
];
880 tmp6
[4] = state
->boxv
[ZZ
][XX
];
881 tmp6
[5] = state
->boxv
[ZZ
][YY
];
882 add_ebin(md
->ebin
,md
->ipc
,md
->bTricl
? 6 : 3,tmp6
,bSum
);
886 add_ebin(md
->ebin
,md
->imu
,3,mu_tot
,bSum
);
888 if (ekind
&& ekind
->cosacc
.cos_accel
!= 0)
890 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
891 dens
= (tmass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
892 add_ebin(md
->ebin
,md
->ivcos
,1,&(ekind
->cosacc
.vcos
),bSum
);
893 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
894 tmp
= 1/(ekind
->cosacc
.cos_accel
/(ekind
->cosacc
.vcos
*PICO
)
895 *dens
*vol
*sqr(box
[ZZ
][ZZ
]*NANO
/(2*M_PI
)));
896 add_ebin(md
->ebin
,md
->ivisc
,1,&tmp
,bSum
);
901 for(i
=0; (i
<md
->nEg
); i
++)
903 for(j
=i
; (j
<md
->nEg
); j
++)
905 gid
=GID(i
,j
,md
->nEg
);
906 for(k
=kk
=0; (k
<egNR
); k
++)
910 eee
[kk
++] = enerd
->grpp
.ener
[k
][gid
];
913 add_ebin(md
->ebin
,md
->igrp
[n
],md
->nEc
,eee
,bSum
);
921 for(i
=0; (i
<md
->nTC
); i
++)
923 md
->tmp_r
[i
] = ekind
->tcstat
[i
].T
;
925 add_ebin(md
->ebin
,md
->itemp
,md
->nTC
,md
->tmp_r
,bSum
);
927 if (md
->etc
== etcNOSEHOOVER
)
929 /* whether to print Nose-Hoover chains: */
930 if (md
->bPrintNHChains
)
932 if (md
->bNHC_trotter
)
934 for(i
=0; (i
<md
->nTC
); i
++)
936 for (j
=0;j
<md
->nNHC
;j
++)
939 md
->tmp_r
[2*k
] = state
->nosehoover_xi
[k
];
940 md
->tmp_r
[2*k
+1] = state
->nosehoover_vxi
[k
];
943 add_ebin(md
->ebin
,md
->itc
,md
->mde_n
,md
->tmp_r
,bSum
);
946 for(i
=0; (i
<md
->nTCP
); i
++)
948 for (j
=0;j
<md
->nNHC
;j
++)
951 md
->tmp_r
[2*k
] = state
->nhpres_xi
[k
];
952 md
->tmp_r
[2*k
+1] = state
->nhpres_vxi
[k
];
955 add_ebin(md
->ebin
,md
->itcb
,md
->mdeb_n
,md
->tmp_r
,bSum
);
960 for(i
=0; (i
<md
->nTC
); i
++)
962 md
->tmp_r
[2*i
] = state
->nosehoover_xi
[i
];
963 md
->tmp_r
[2*i
+1] = state
->nosehoover_vxi
[i
];
965 add_ebin(md
->ebin
,md
->itc
,md
->mde_n
,md
->tmp_r
,bSum
);
969 else if (md
->etc
== etcBERENDSEN
|| md
->etc
== etcYES
||
970 md
->etc
== etcVRESCALE
)
972 for(i
=0; (i
<md
->nTC
); i
++)
974 md
->tmp_r
[i
] = ekind
->tcstat
[i
].lambda
;
976 add_ebin(md
->ebin
,md
->itc
,md
->nTC
,md
->tmp_r
,bSum
);
980 if (ekind
&& md
->nU
> 1)
982 for(i
=0; (i
<md
->nU
); i
++)
984 copy_rvec(ekind
->grpstat
[i
].u
,md
->tmp_v
[i
]);
986 add_ebin(md
->ebin
,md
->iu
,3*md
->nU
,md
->tmp_v
[0],bSum
);
989 ebin_increase_count(md
->ebin
,bSum
);
991 /* BAR + thermodynamic integration values */
992 if ((md
->fp_dhdl
|| md
->dhc
) && bDoDHDL
&& (enerd
->n_lambda
> 0))
994 snew(dE
,enerd
->n_lambda
-1);
995 for(i
=0; i
<enerd
->n_lambda
-1; i
++) {
996 dE
[i
] = enerd
->enerpart_lambda
[i
+1]-enerd
->enerpart_lambda
[0]; /* zero for simulated tempering */
997 if (md
->temperatures
!=NULL
)
999 /* MRS: is this right, given the way we have defined the exchange probabilities? */
1000 /* is this even useful to have at all? */
1001 dE
[i
] += (md
->temperatures
[i
]/md
->temperatures
[state
->fep_state
]-1.0)*enerd
->term
[F_EKIN
];
1006 if (md
->fp_dhdl
&& bDoDHDL
)
1008 fprintf(md
->fp_dhdl
,"%.4f",time
);
1009 /* the current free energy state */
1011 /* print the current state if we are doing expanded ensemble */
1012 if (expand
->elmcmove
> elmcmoveNO
) {
1013 fprintf(md
->fp_dhdl
," %4d",state
->fep_state
);
1015 /* total energy (for if the temperature changes */
1016 if (fep
->bPrintEnergy
)
1018 store_energy
= enerd
->term
[F_ETOT
];
1019 fprintf(md
->fp_dhdl
," %#.8g",store_energy
);
1022 for (i
=0;i
<efptNR
;i
++)
1024 if (fep
->separate_dvdl
[i
])
1026 fprintf(md
->fp_dhdl
," %#.8g",enerd
->term
[F_DVDL
+i
]); /* assumes F_DVDL is first */
1029 for(i
=1; i
<enerd
->n_lambda
; i
++)
1031 fprintf(md
->fp_dhdl
," %#.8g",dE
[i
-1]);
1034 if ((md
->epc
!=epcNO
) && (enerd
->n_lambda
> 0))
1036 fprintf(md
->fp_dhdl
," %#.8g",pv
); /* PV term only needed when there are alternate state lambda */
1038 fprintf(md
->fp_dhdl
,"\n");
1039 /* and the binary free energy output */
1041 if (md
->dhc
&& bDoDHDL
)
1044 for (i
=0;i
<efptNR
;i
++)
1046 if (fep
->separate_dvdl
[i
])
1048 store_dhdl
[idhdl
] = enerd
->term
[F_DVDL
+i
]; /* assumes F_DVDL is first */
1052 /* store_dh is dE */
1053 mde_delta_h_coll_add_dh(md
->dhc
,
1054 (double)state
->fep_state
,
1057 (expand
->elamstats
>elamstatsNO
),
1058 (fep
->bPrintEnergy
),
1066 if ((md
->fp_dhdl
|| md
->dhc
) && bDoDHDL
&& (enerd
->n_lambda
>0))
1073 void upd_mdebin_step(t_mdebin
*md
)
1075 ebin_increase_count(md
->ebin
,FALSE
);
1078 static void npr(FILE *log
,int n
,char c
)
1080 for(; (n
>0); n
--) fprintf(log
,"%c",c
);
1083 static void pprint(FILE *log
,const char *s
,t_mdebin
*md
)
1087 char buf1
[22],buf2
[22];
1090 fprintf(log
,"\t<====== ");
1092 fprintf(log
," ==>\n");
1093 fprintf(log
,"\t<==== %s ====>\n",s
);
1094 fprintf(log
,"\t<== ");
1096 fprintf(log
," ======>\n\n");
1098 fprintf(log
,"\tStatistics over %s steps using %s frames\n",
1099 gmx_step_str(md
->ebin
->nsteps_sim
,buf1
),
1100 gmx_step_str(md
->ebin
->nsum_sim
,buf2
));
1104 void print_ebin_header(FILE *log
,gmx_large_int_t steps
,double time
,real lambda
)
1108 fprintf(log
," %12s %12s %12s\n"
1109 " %12s %12.5f %12.5f\n\n",
1110 "Step","Time","Lambda",gmx_step_str(steps
,buf
),time
,lambda
);
1113 void print_ebin(ener_file_t fp_ene
,gmx_bool bEne
,gmx_bool bDR
,gmx_bool bOR
,
1115 gmx_large_int_t step
,double time
,
1116 int mode
,gmx_bool bCompact
,
1117 t_mdebin
*md
,t_fcdata
*fcd
,
1118 gmx_groups_t
*groups
,t_grpopts
*opts
)
1120 /*static char **grpnms=NULL;*/
1122 int i
,j
,n
,ni
,nj
,ndr
,nor
,b
;
1124 real
*disre_rm3tav
, *disre_rt
;
1126 /* these are for the old-style blocks (1 subblock, only reals), because
1127 there can be only one per ID for these */
1132 /* temporary arrays for the lambda values to write out */
1133 double enxlambda_data
[2];
1143 fr
.nsteps
= md
->ebin
->nsteps
;
1144 fr
.dt
= md
->delta_t
;
1145 fr
.nsum
= md
->ebin
->nsum
;
1146 fr
.nre
= (bEne
) ? md
->ebin
->nener
: 0;
1147 fr
.ener
= md
->ebin
->e
;
1148 ndisre
= bDR
? fcd
->disres
.npair
: 0;
1149 disre_rm3tav
= fcd
->disres
.rm3tav
;
1150 disre_rt
= fcd
->disres
.rt
;
1151 /* Optional additional old-style (real-only) blocks. */
1152 for(i
=0; i
<enxNR
; i
++)
1156 if (fcd
->orires
.nr
> 0 && bOR
)
1158 diagonalize_orires_tensors(&(fcd
->orires
));
1159 nr
[enxOR
] = fcd
->orires
.nr
;
1160 block
[enxOR
] = fcd
->orires
.otav
;
1162 nr
[enxORI
] = (fcd
->orires
.oinsl
!= fcd
->orires
.otav
) ?
1164 block
[enxORI
] = fcd
->orires
.oinsl
;
1165 id
[enxORI
] = enxORI
;
1166 nr
[enxORT
] = fcd
->orires
.nex
*12;
1167 block
[enxORT
] = fcd
->orires
.eig
;
1168 id
[enxORT
] = enxORT
;
1171 /* whether we are going to wrte anything out: */
1172 if (fr
.nre
|| ndisre
|| nr
[enxOR
] || nr
[enxORI
])
1175 /* the old-style blocks go first */
1177 for(i
=0; i
<enxNR
; i
++)
1184 add_blocks_enxframe(&fr
, fr
.nblock
);
1185 for(b
=0;b
<fr
.nblock
;b
++)
1187 add_subblocks_enxblock(&(fr
.block
[b
]), 1);
1188 fr
.block
[b
].id
=id
[b
];
1189 fr
.block
[b
].sub
[0].nr
= nr
[b
];
1191 fr
.block
[b
].sub
[0].type
= xdr_datatype_float
;
1192 fr
.block
[b
].sub
[0].fval
= block
[b
];
1194 fr
.block
[b
].sub
[0].type
= xdr_datatype_double
;
1195 fr
.block
[b
].sub
[0].dval
= block
[b
];
1199 /* check for disre block & fill it. */
1204 add_blocks_enxframe(&fr
, fr
.nblock
);
1206 add_subblocks_enxblock(&(fr
.block
[db
]), 2);
1207 fr
.block
[db
].id
=enxDISRE
;
1208 fr
.block
[db
].sub
[0].nr
=ndisre
;
1209 fr
.block
[db
].sub
[1].nr
=ndisre
;
1211 fr
.block
[db
].sub
[0].type
=xdr_datatype_float
;
1212 fr
.block
[db
].sub
[1].type
=xdr_datatype_float
;
1213 fr
.block
[db
].sub
[0].fval
=disre_rt
;
1214 fr
.block
[db
].sub
[1].fval
=disre_rm3tav
;
1216 fr
.block
[db
].sub
[0].type
=xdr_datatype_double
;
1217 fr
.block
[db
].sub
[1].type
=xdr_datatype_double
;
1218 fr
.block
[db
].sub
[0].dval
=disre_rt
;
1219 fr
.block
[db
].sub
[1].dval
=disre_rm3tav
;
1222 /* here we can put new-style blocks */
1224 /* Free energy perturbation blocks */
1227 mde_delta_h_coll_handle_block(md
->dhc
, &fr
, fr
.nblock
);
1230 /* we can now free & reset the data in the blocks */
1233 mde_delta_h_coll_reset(md
->dhc
);
1236 /* do the actual I/O */
1238 gmx_fio_check_file_position(enx_file_pointer(fp_ene
));
1241 /* We have stored the sums, so reset the sum history */
1242 reset_ebin_sums(md
->ebin
);
1250 pprint(log
,"A V E R A G E S",md
);
1256 pprint(log
,"R M S - F L U C T U A T I O N S",md
);
1260 gmx_fatal(FARGS
,"Invalid print mode (%d)",mode
);
1265 for(i
=0;i
<opts
->ngtc
;i
++)
1267 if(opts
->annealing
[i
]!=eannNO
)
1269 fprintf(log
,"Current ref_t for group %s: %8.1f\n",
1270 *(groups
->grpname
[groups
->grps
[egcTC
].nm_ind
[i
]]),
1274 if (mode
==eprNORMAL
&& fcd
->orires
.nr
>0)
1276 print_orires_log(log
,&(fcd
->orires
));
1278 fprintf(log
," Energies (%s)\n",unit_energy
);
1279 pr_ebin(log
,md
->ebin
,md
->ie
,md
->f_nre
+md
->nCrmsd
,5,mode
,TRUE
);
1286 pr_ebin(log
,md
->ebin
,md
->ib
, md
->bTricl
? NTRICLBOXS
: NBOXS
,5,
1292 fprintf(log
," Constraint Virial (%s)\n",unit_energy
);
1293 pr_ebin(log
,md
->ebin
,md
->isvir
,9,3,mode
,FALSE
);
1295 fprintf(log
," Force Virial (%s)\n",unit_energy
);
1296 pr_ebin(log
,md
->ebin
,md
->ifvir
,9,3,mode
,FALSE
);
1301 fprintf(log
," Total Virial (%s)\n",unit_energy
);
1302 pr_ebin(log
,md
->ebin
,md
->ivir
,9,3,mode
,FALSE
);
1307 fprintf(log
," Pressure (%s)\n",unit_pres_bar
);
1308 pr_ebin(log
,md
->ebin
,md
->ipres
,9,3,mode
,FALSE
);
1313 fprintf(log
," Total Dipole (%s)\n",unit_dipole_D
);
1314 pr_ebin(log
,md
->ebin
,md
->imu
,3,3,mode
,FALSE
);
1320 if (md
->print_grpnms
==NULL
)
1322 snew(md
->print_grpnms
,md
->nE
);
1324 for(i
=0; (i
<md
->nEg
); i
++)
1326 ni
=groups
->grps
[egcENER
].nm_ind
[i
];
1327 for(j
=i
; (j
<md
->nEg
); j
++)
1329 nj
=groups
->grps
[egcENER
].nm_ind
[j
];
1330 sprintf(buf
,"%s-%s",*(groups
->grpname
[ni
]),
1331 *(groups
->grpname
[nj
]));
1332 md
->print_grpnms
[n
++]=strdup(buf
);
1336 sprintf(buf
,"Epot (%s)",unit_energy
);
1337 fprintf(log
,"%15s ",buf
);
1338 for(i
=0; (i
<egNR
); i
++)
1342 fprintf(log
,"%12s ",egrp_nm
[i
]);
1346 for(i
=0; (i
<md
->nE
); i
++)
1348 fprintf(log
,"%15s",md
->print_grpnms
[i
]);
1349 pr_ebin(log
,md
->ebin
,md
->igrp
[i
],md
->nEc
,md
->nEc
,mode
,
1356 pr_ebin(log
,md
->ebin
,md
->itemp
,md
->nTC
,4,mode
,TRUE
);
1361 fprintf(log
,"%15s %12s %12s %12s\n",
1362 "Group","Ux","Uy","Uz");
1363 for(i
=0; (i
<md
->nU
); i
++)
1365 ni
=groups
->grps
[egcACC
].nm_ind
[i
];
1366 fprintf(log
,"%15s",*groups
->grpname
[ni
]);
1367 pr_ebin(log
,md
->ebin
,md
->iu
+3*i
,3,3,mode
,FALSE
);
1376 void update_energyhistory(energyhistory_t
* enerhist
,t_mdebin
* mdebin
)
1380 enerhist
->nsteps
= mdebin
->ebin
->nsteps
;
1381 enerhist
->nsum
= mdebin
->ebin
->nsum
;
1382 enerhist
->nsteps_sim
= mdebin
->ebin
->nsteps_sim
;
1383 enerhist
->nsum_sim
= mdebin
->ebin
->nsum_sim
;
1384 enerhist
->nener
= mdebin
->ebin
->nener
;
1386 if (mdebin
->ebin
->nsum
> 0)
1388 /* Check if we need to allocate first */
1389 if(enerhist
->ener_ave
== NULL
)
1391 snew(enerhist
->ener_ave
,enerhist
->nener
);
1392 snew(enerhist
->ener_sum
,enerhist
->nener
);
1395 for(i
=0;i
<enerhist
->nener
;i
++)
1397 enerhist
->ener_ave
[i
] = mdebin
->ebin
->e
[i
].eav
;
1398 enerhist
->ener_sum
[i
] = mdebin
->ebin
->e
[i
].esum
;
1402 if (mdebin
->ebin
->nsum_sim
> 0)
1404 /* Check if we need to allocate first */
1405 if(enerhist
->ener_sum_sim
== NULL
)
1407 snew(enerhist
->ener_sum_sim
,enerhist
->nener
);
1410 for(i
=0;i
<enerhist
->nener
;i
++)
1412 enerhist
->ener_sum_sim
[i
] = mdebin
->ebin
->e_sim
[i
].esum
;
1417 mde_delta_h_coll_update_energyhistory(mdebin
->dhc
, enerhist
);
1421 void restore_energyhistory_from_state(t_mdebin
* mdebin
,
1422 energyhistory_t
* enerhist
)
1426 if ((enerhist
->nsum
> 0 || enerhist
->nsum_sim
> 0) &&
1427 mdebin
->ebin
->nener
!= enerhist
->nener
)
1429 gmx_fatal(FARGS
,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1430 mdebin
->ebin
->nener
,enerhist
->nener
);
1433 mdebin
->ebin
->nsteps
= enerhist
->nsteps
;
1434 mdebin
->ebin
->nsum
= enerhist
->nsum
;
1435 mdebin
->ebin
->nsteps_sim
= enerhist
->nsteps_sim
;
1436 mdebin
->ebin
->nsum_sim
= enerhist
->nsum_sim
;
1438 for(i
=0; i
<mdebin
->ebin
->nener
; i
++)
1440 mdebin
->ebin
->e
[i
].eav
=
1441 (enerhist
->nsum
> 0 ? enerhist
->ener_ave
[i
] : 0);
1442 mdebin
->ebin
->e
[i
].esum
=
1443 (enerhist
->nsum
> 0 ? enerhist
->ener_sum
[i
] : 0);
1444 mdebin
->ebin
->e_sim
[i
].esum
=
1445 (enerhist
->nsum_sim
> 0 ? enerhist
->ener_sum_sim
[i
] : 0);
1449 mde_delta_h_coll_restore_energyhistory(mdebin
->dhc
, enerhist
);