File: | gromacs/mdlib/mdebin.c |
Location: | line 1237, column 17 |
Description: | Function call argument is an uninitialized value |
1 | /* | |||
2 | * This file is part of the GROMACS molecular simulation package. | |||
3 | * | |||
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. | |||
5 | * Copyright (c) 2001-2004, The GROMACS development team. | |||
6 | * Copyright (c) 2013,2014, 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. | |||
10 | * | |||
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. | |||
15 | * | |||
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. | |||
20 | * | |||
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. | |||
25 | * | |||
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. | |||
33 | * | |||
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. | |||
36 | */ | |||
37 | #ifdef HAVE_CONFIG_H1 | |||
38 | #include <config.h> | |||
39 | #endif | |||
40 | ||||
41 | #include <float.h> | |||
42 | #include <stdlib.h> | |||
43 | #include <string.h> | |||
44 | ||||
45 | #include "typedefs.h" | |||
46 | #include "mdebin.h" | |||
47 | #include "gromacs/utility/smalloc.h" | |||
48 | #include "physics.h" | |||
49 | #include "gromacs/fileio/enxio.h" | |||
50 | #include "gromacs/math/vec.h" | |||
51 | #include "disre.h" | |||
52 | #include "network.h" | |||
53 | #include "names.h" | |||
54 | #include "orires.h" | |||
55 | #include "constr.h" | |||
56 | #include "mtop_util.h" | |||
57 | #include "gromacs/fileio/xvgr.h" | |||
58 | #include "gromacs/fileio/gmxfio.h" | |||
59 | #include "macros.h" | |||
60 | #include "mdrun.h" | |||
61 | #include "mdebin_bar.h" | |||
62 | ||||
63 | ||||
64 | static const char *conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" }; | |||
65 | ||||
66 | static const char *boxs_nm[] = { "Box-X", "Box-Y", "Box-Z" }; | |||
67 | ||||
68 | static const char *tricl_boxs_nm[] = { | |||
69 | "Box-XX", "Box-YY", "Box-ZZ", | |||
70 | "Box-YX", "Box-ZX", "Box-ZY" | |||
71 | }; | |||
72 | ||||
73 | static const char *vol_nm[] = { "Volume" }; | |||
74 | ||||
75 | static const char *dens_nm[] = {"Density" }; | |||
76 | ||||
77 | static const char *pv_nm[] = {"pV" }; | |||
78 | ||||
79 | static const char *enthalpy_nm[] = {"Enthalpy" }; | |||
80 | ||||
81 | static const char *boxvel_nm[] = { | |||
82 | "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ", | |||
83 | "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY" | |||
84 | }; | |||
85 | ||||
86 | #define NBOXS((int)(sizeof(boxs_nm)/sizeof((boxs_nm)[0]))) asize(boxs_nm)((int)(sizeof(boxs_nm)/sizeof((boxs_nm)[0]))) | |||
87 | #define NTRICLBOXS((int)(sizeof(tricl_boxs_nm)/sizeof((tricl_boxs_nm)[0]))) asize(tricl_boxs_nm)((int)(sizeof(tricl_boxs_nm)/sizeof((tricl_boxs_nm)[0]))) | |||
88 | ||||
89 | t_mdebin *init_mdebin(ener_file_t fp_ene, | |||
90 | const gmx_mtop_t *mtop, | |||
91 | const t_inputrec *ir, | |||
92 | FILE *fp_dhdl) | |||
93 | { | |||
94 | const char *ener_nm[F_NRE]; | |||
95 | static const char *vir_nm[] = { | |||
96 | "Vir-XX", "Vir-XY", "Vir-XZ", | |||
97 | "Vir-YX", "Vir-YY", "Vir-YZ", | |||
98 | "Vir-ZX", "Vir-ZY", "Vir-ZZ" | |||
99 | }; | |||
100 | static const char *sv_nm[] = { | |||
101 | "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ", | |||
102 | "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ", | |||
103 | "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ" | |||
104 | }; | |||
105 | static const char *fv_nm[] = { | |||
106 | "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ", | |||
107 | "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ", | |||
108 | "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ" | |||
109 | }; | |||
110 | static const char *pres_nm[] = { | |||
111 | "Pres-XX", "Pres-XY", "Pres-XZ", | |||
112 | "Pres-YX", "Pres-YY", "Pres-YZ", | |||
113 | "Pres-ZX", "Pres-ZY", "Pres-ZZ" | |||
114 | }; | |||
115 | static const char *surft_nm[] = { | |||
116 | "#Surf*SurfTen" | |||
117 | }; | |||
118 | static const char *mu_nm[] = { | |||
119 | "Mu-X", "Mu-Y", "Mu-Z" | |||
120 | }; | |||
121 | static const char *vcos_nm[] = { | |||
122 | "2CosZ*Vel-X" | |||
123 | }; | |||
124 | static const char *visc_nm[] = { | |||
125 | "1/Viscosity" | |||
126 | }; | |||
127 | static const char *baro_nm[] = { | |||
128 | "Barostat" | |||
129 | }; | |||
130 | ||||
131 | char **grpnms; | |||
132 | const gmx_groups_t *groups; | |||
133 | char **gnm; | |||
134 | char buf[256]; | |||
135 | const char *bufi; | |||
136 | t_mdebin *md; | |||
137 | int i, j, ni, nj, n, nh, k, kk, ncon, nset; | |||
138 | gmx_bool bBHAM, bNoseHoover, b14; | |||
139 | ||||
140 | snew(md, 1)(md) = save_calloc("md", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 140, (1), sizeof(*(md))); | |||
141 | ||||
142 | md->bVir = TRUE1; | |||
143 | md->bPress = TRUE1; | |||
144 | md->bSurft = TRUE1; | |||
145 | md->bMu = TRUE1; | |||
146 | ||||
147 | if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || (ir->eI) == eiBD)) | |||
148 | { | |||
149 | md->delta_t = ir->delta_t; | |||
150 | } | |||
151 | else | |||
152 | { | |||
153 | md->delta_t = 0; | |||
154 | } | |||
155 | ||||
156 | groups = &mtop->groups; | |||
157 | ||||
158 | bBHAM = (mtop->ffparams.functype[0] == F_BHAM); | |||
159 | b14 = (gmx_mtop_ftype_count(mtop, F_LJ14) > 0 || | |||
160 | gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0); | |||
161 | ||||
162 | ncon = gmx_mtop_ftype_count(mtop, F_CONSTR); | |||
163 | nset = gmx_mtop_ftype_count(mtop, F_SETTLE); | |||
164 | md->bConstr = (ncon > 0 || nset > 0); | |||
165 | md->bConstrVir = FALSE0; | |||
166 | if (md->bConstr) | |||
167 | { | |||
168 | if (ncon > 0 && ir->eConstrAlg == econtLINCS) | |||
169 | { | |||
170 | if (ir->eI == eiSD2) | |||
171 | { | |||
172 | md->nCrmsd = 2; | |||
173 | } | |||
174 | else | |||
175 | { | |||
176 | md->nCrmsd = 1; | |||
177 | } | |||
178 | } | |||
179 | md->bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL((void*)0)); | |||
180 | } | |||
181 | else | |||
182 | { | |||
183 | md->nCrmsd = 0; | |||
184 | } | |||
185 | ||||
186 | /* Energy monitoring */ | |||
187 | for (i = 0; i < egNR; i++) | |||
188 | { | |||
189 | md->bEInd[i] = FALSE0; | |||
190 | } | |||
191 | ||||
192 | for (i = 0; i < F_NRE; i++) | |||
193 | { | |||
194 | md->bEner[i] = FALSE0; | |||
195 | if (i == F_LJ) | |||
196 | { | |||
197 | md->bEner[i] = !bBHAM; | |||
198 | } | |||
199 | else if (i == F_BHAM) | |||
200 | { | |||
201 | md->bEner[i] = bBHAM; | |||
202 | } | |||
203 | else if (i == F_EQM) | |||
204 | { | |||
205 | md->bEner[i] = ir->bQMMM; | |||
206 | } | |||
207 | else if (i == F_COUL_LR) | |||
208 | { | |||
209 | md->bEner[i] = (ir->rcoulomb > ir->rlist); | |||
210 | } | |||
211 | else if (i == F_LJ_LR) | |||
212 | { | |||
213 | md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist); | |||
214 | } | |||
215 | else if (i == F_BHAM_LR) | |||
216 | { | |||
217 | md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist); | |||
218 | } | |||
219 | else if (i == F_RF_EXCL) | |||
220 | { | |||
221 | md->bEner[i] = (EEL_RF(ir->coulombtype)((ir->coulombtype) == eelRF || (ir->coulombtype) == eelGRF || (ir->coulombtype) == eelRF_NEC || (ir->coulombtype) == eelRF_ZERO ) && ir->coulombtype != eelRF_NEC && ir->cutoff_scheme == ecutsGROUP); | |||
222 | } | |||
223 | else if (i == F_COUL_RECIP) | |||
224 | { | |||
225 | md->bEner[i] = EEL_FULL(ir->coulombtype)((((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH || (ir->coulombtype) == eelPMEUSER || (ir->coulombtype ) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD) || (ir->coulombtype) == eelEWALD) || (ir->coulombtype) == eelPOISSON); | |||
226 | } | |||
227 | else if (i == F_LJ_RECIP) | |||
228 | { | |||
229 | md->bEner[i] = EVDW_PME(ir->vdwtype)((ir->vdwtype) == evdwPME); | |||
230 | } | |||
231 | else if (i == F_LJ14) | |||
232 | { | |||
233 | md->bEner[i] = b14; | |||
234 | } | |||
235 | else if (i == F_COUL14) | |||
236 | { | |||
237 | md->bEner[i] = b14; | |||
238 | } | |||
239 | else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB) | |||
240 | { | |||
241 | md->bEner[i] = FALSE0; | |||
242 | } | |||
243 | else if ((i == F_DVDL_COUL && ir->fepvals->separate_dvdl[efptCOUL]) || | |||
244 | (i == F_DVDL_VDW && ir->fepvals->separate_dvdl[efptVDW]) || | |||
245 | (i == F_DVDL_BONDED && ir->fepvals->separate_dvdl[efptBONDED]) || | |||
246 | (i == F_DVDL_RESTRAINT && ir->fepvals->separate_dvdl[efptRESTRAINT]) || | |||
247 | (i == F_DKDL && ir->fepvals->separate_dvdl[efptMASS]) || | |||
248 | (i == F_DVDL && ir->fepvals->separate_dvdl[efptFEP])) | |||
249 | { | |||
250 | md->bEner[i] = (ir->efep != efepNO); | |||
251 | } | |||
252 | else if ((interaction_function[i].flags & IF_VSITE1<<1) || | |||
253 | (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE)) | |||
254 | { | |||
255 | md->bEner[i] = FALSE0; | |||
256 | } | |||
257 | else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES) || (i == F_EQM)) | |||
258 | { | |||
259 | md->bEner[i] = TRUE1; | |||
260 | } | |||
261 | else if ((i == F_GBPOL) && ir->implicit_solvent == eisGBSA) | |||
262 | { | |||
263 | md->bEner[i] = TRUE1; | |||
264 | } | |||
265 | else if ((i == F_NPSOLVATION) && ir->implicit_solvent == eisGBSA && (ir->sa_algorithm != esaNO)) | |||
266 | { | |||
267 | md->bEner[i] = TRUE1; | |||
268 | } | |||
269 | else if ((i == F_GB12) || (i == F_GB13) || (i == F_GB14)) | |||
270 | { | |||
271 | md->bEner[i] = FALSE0; | |||
272 | } | |||
273 | else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP)) | |||
274 | { | |||
275 | md->bEner[i] = EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || (ir->eI) == eiBD); | |||
276 | } | |||
277 | else if (i == F_DISPCORR || i == F_PDISPCORR) | |||
278 | { | |||
279 | md->bEner[i] = (ir->eDispCorr != edispcNO); | |||
280 | } | |||
281 | else if (i == F_DISRESVIOL) | |||
282 | { | |||
283 | md->bEner[i] = (gmx_mtop_ftype_count(mtop, F_DISRES) > 0); | |||
284 | } | |||
285 | else if (i == F_ORIRESDEV) | |||
286 | { | |||
287 | md->bEner[i] = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0); | |||
288 | } | |||
289 | else if (i == F_CONNBONDS) | |||
290 | { | |||
291 | md->bEner[i] = FALSE0; | |||
292 | } | |||
293 | else if (i == F_COM_PULL) | |||
294 | { | |||
295 | md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F || ir->bRot); | |||
296 | } | |||
297 | else if (i == F_ECONSERVED) | |||
298 | { | |||
299 | md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) && | |||
300 | (ir->epc == epcNO || ir->epc == epcMTTK)); | |||
301 | } | |||
302 | else | |||
303 | { | |||
304 | md->bEner[i] = (gmx_mtop_ftype_count(mtop, i) > 0); | |||
305 | } | |||
306 | } | |||
307 | ||||
308 | /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/ | |||
309 | if (ir->bAdress && !debug) | |||
310 | { | |||
311 | for (i = 0; i < F_NRE; i++) | |||
312 | { | |||
313 | md->bEner[i] = FALSE0; | |||
314 | if (i == F_EKIN) | |||
315 | { | |||
316 | md->bEner[i] = TRUE1; | |||
317 | } | |||
318 | if (i == F_TEMP) | |||
319 | { | |||
320 | md->bEner[i] = TRUE1; | |||
321 | } | |||
322 | } | |||
323 | md->bVir = FALSE0; | |||
324 | md->bPress = FALSE0; | |||
325 | md->bSurft = FALSE0; | |||
326 | md->bMu = FALSE0; | |||
327 | } | |||
328 | ||||
329 | md->f_nre = 0; | |||
330 | for (i = 0; i < F_NRE; i++) | |||
331 | { | |||
332 | if (md->bEner[i]) | |||
333 | { | |||
334 | ener_nm[md->f_nre] = interaction_function[i].longname; | |||
335 | md->f_nre++; | |||
336 | } | |||
337 | } | |||
338 | ||||
339 | md->epc = ir->epc; | |||
340 | md->bDiagPres = !TRICLINIC(ir->ref_p)(ir->ref_p[1][0] != 0 || ir->ref_p[2][0] != 0 || ir-> ref_p[2][1] != 0); | |||
341 | md->ref_p = (ir->ref_p[XX0][XX0]+ir->ref_p[YY1][YY1]+ir->ref_p[ZZ2][ZZ2])/DIM3; | |||
342 | md->bTricl = TRICLINIC(ir->compress)(ir->compress[1][0] != 0 || ir->compress[2][0] != 0 || ir ->compress[2][1] != 0) || TRICLINIC(ir->deform)(ir->deform[1][0] != 0 || ir->deform[2][0] != 0 || ir-> deform[2][1] != 0); | |||
343 | md->bDynBox = DYNAMIC_BOX(*ir)((*ir).epc != epcNO || (*ir).eI == eiTPI || ((*ir).deform[0][ 0] != 0 || (*ir).deform[1][1] != 0 || (*ir).deform[2][2] != 0 || (*ir).deform[1][0] != 0 || (*ir).deform[2][0] != 0 || (*ir ).deform[2][1] != 0)); | |||
344 | md->etc = ir->etc; | |||
345 | md->bNHC_trotter = IR_NVT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((!((ir)->epc == epcMTTK)) && ((ir)->etc == etcNOSEHOOVER ))); | |||
346 | md->bPrintNHChains = ir->bPrintNHChains; | |||
347 | md->bMTTK = (IR_NPT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && ((ir)->etc == etcNOSEHOOVER ))) || IR_NPH_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && (!(((ir)->etc == etcNOSEHOOVER )))))); | |||
348 | md->bMu = NEED_MUTOT(*ir)(((*ir).coulombtype == eelEWALD || (((*ir).coulombtype) == eelPME || ((*ir).coulombtype) == eelPMESWITCH || ((*ir).coulombtype ) == eelPMEUSER || ((*ir).coulombtype) == eelPMEUSERSWITCH || ((*ir).coulombtype) == eelP3M_AD)) && ((*ir).ewald_geometry == eewg3DC || (*ir).epsilon_surface != 0)); | |||
349 | ||||
350 | md->ebin = mk_ebin(); | |||
351 | /* Pass NULL for unit to let get_ebin_space determine the units | |||
352 | * for interaction_function[i].longname | |||
353 | */ | |||
354 | md->ie = get_ebin_space(md->ebin, md->f_nre, ener_nm, NULL((void*)0)); | |||
355 | if (md->nCrmsd) | |||
356 | { | |||
357 | /* This should be called directly after the call for md->ie, | |||
358 | * such that md->iconrmsd follows directly in the list. | |||
359 | */ | |||
360 | md->iconrmsd = get_ebin_space(md->ebin, md->nCrmsd, conrmsd_nm, ""); | |||
361 | } | |||
362 | if (md->bDynBox) | |||
363 | { | |||
364 | md->ib = get_ebin_space(md->ebin, | |||
365 | md->bTricl ? NTRICLBOXS((int)(sizeof(tricl_boxs_nm)/sizeof((tricl_boxs_nm)[0]))) : NBOXS((int)(sizeof(boxs_nm)/sizeof((boxs_nm)[0]))), | |||
366 | md->bTricl ? tricl_boxs_nm : boxs_nm, | |||
367 | unit_length"nm"); | |||
368 | md->ivol = get_ebin_space(md->ebin, 1, vol_nm, unit_volume"nm" "^3"); | |||
369 | md->idens = get_ebin_space(md->ebin, 1, dens_nm, unit_density_SI"kg" "/" "m" "^3"); | |||
370 | if (md->bDiagPres) | |||
371 | { | |||
372 | md->ipv = get_ebin_space(md->ebin, 1, pv_nm, unit_energy"kJ/mol"); | |||
373 | md->ienthalpy = get_ebin_space(md->ebin, 1, enthalpy_nm, unit_energy"kJ/mol"); | |||
374 | } | |||
375 | } | |||
376 | if (md->bConstrVir) | |||
377 | { | |||
378 | md->isvir = get_ebin_space(md->ebin, asize(sv_nm)((int)(sizeof(sv_nm)/sizeof((sv_nm)[0]))), sv_nm, unit_energy"kJ/mol"); | |||
379 | md->ifvir = get_ebin_space(md->ebin, asize(fv_nm)((int)(sizeof(fv_nm)/sizeof((fv_nm)[0]))), fv_nm, unit_energy"kJ/mol"); | |||
380 | } | |||
381 | if (md->bVir) | |||
382 | { | |||
383 | md->ivir = get_ebin_space(md->ebin, asize(vir_nm)((int)(sizeof(vir_nm)/sizeof((vir_nm)[0]))), vir_nm, unit_energy"kJ/mol"); | |||
384 | } | |||
385 | if (md->bPress) | |||
386 | { | |||
387 | md->ipres = get_ebin_space(md->ebin, asize(pres_nm)((int)(sizeof(pres_nm)/sizeof((pres_nm)[0]))), pres_nm, unit_pres_bar"bar"); | |||
388 | } | |||
389 | if (md->bSurft) | |||
390 | { | |||
391 | md->isurft = get_ebin_space(md->ebin, asize(surft_nm)((int)(sizeof(surft_nm)/sizeof((surft_nm)[0]))), surft_nm, | |||
392 | unit_surft_bar"bar" " " "nm"); | |||
393 | } | |||
394 | if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK) | |||
395 | { | |||
396 | md->ipc = get_ebin_space(md->ebin, md->bTricl ? 6 : 3, | |||
397 | boxvel_nm, unit_vel"nm" "/" "ps"); | |||
398 | } | |||
399 | if (md->bMu) | |||
400 | { | |||
401 | md->imu = get_ebin_space(md->ebin, asize(mu_nm)((int)(sizeof(mu_nm)/sizeof((mu_nm)[0]))), mu_nm, unit_dipole_D"D"); | |||
402 | } | |||
403 | if (ir->cos_accel != 0) | |||
404 | { | |||
405 | md->ivcos = get_ebin_space(md->ebin, asize(vcos_nm)((int)(sizeof(vcos_nm)/sizeof((vcos_nm)[0]))), vcos_nm, unit_vel"nm" "/" "ps"); | |||
406 | md->ivisc = get_ebin_space(md->ebin, asize(visc_nm)((int)(sizeof(visc_nm)/sizeof((visc_nm)[0]))), visc_nm, | |||
407 | unit_invvisc_SI"m" " " "s" "/" "kg"); | |||
408 | } | |||
409 | ||||
410 | /* Energy monitoring */ | |||
411 | for (i = 0; i < egNR; i++) | |||
412 | { | |||
413 | md->bEInd[i] = FALSE0; | |||
414 | } | |||
415 | md->bEInd[egCOULSR] = TRUE1; | |||
416 | md->bEInd[egLJSR ] = TRUE1; | |||
417 | ||||
418 | if (ir->rcoulomb > ir->rlist) | |||
419 | { | |||
420 | md->bEInd[egCOULLR] = TRUE1; | |||
421 | } | |||
422 | if (!bBHAM) | |||
423 | { | |||
424 | if (ir->rvdw > ir->rlist) | |||
425 | { | |||
426 | md->bEInd[egLJLR] = TRUE1; | |||
427 | } | |||
428 | } | |||
429 | else | |||
430 | { | |||
431 | md->bEInd[egLJSR] = FALSE0; | |||
432 | md->bEInd[egBHAMSR] = TRUE1; | |||
433 | if (ir->rvdw > ir->rlist) | |||
434 | { | |||
435 | md->bEInd[egBHAMLR] = TRUE1; | |||
436 | } | |||
437 | } | |||
438 | if (b14) | |||
439 | { | |||
440 | md->bEInd[egLJ14] = TRUE1; | |||
441 | md->bEInd[egCOUL14] = TRUE1; | |||
442 | } | |||
443 | md->nEc = 0; | |||
444 | for (i = 0; (i < egNR); i++) | |||
445 | { | |||
446 | if (md->bEInd[i]) | |||
447 | { | |||
448 | md->nEc++; | |||
449 | } | |||
450 | } | |||
451 | ||||
452 | n = groups->grps[egcENER].nr; | |||
453 | /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/ | |||
454 | if (!ir->bAdress) | |||
455 | { | |||
456 | /*standard simulation*/ | |||
457 | md->nEg = n; | |||
458 | md->nE = (n*(n+1))/2; | |||
459 | } | |||
460 | else if (!debug) | |||
461 | { | |||
462 | /*AdResS simulation*/ | |||
463 | md->nU = 0; | |||
464 | md->nEg = 0; | |||
465 | md->nE = 0; | |||
466 | md->nEc = 0; | |||
467 | md->isvir = FALSE0; | |||
468 | } | |||
469 | snew(md->igrp, md->nE)(md->igrp) = save_calloc("md->igrp", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 469, (md->nE), sizeof(*(md->igrp))); | |||
470 | if (md->nE > 1) | |||
471 | { | |||
472 | n = 0; | |||
473 | snew(gnm, md->nEc)(gnm) = save_calloc("gnm", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 473, (md->nEc), sizeof(*(gnm))); | |||
474 | for (k = 0; (k < md->nEc); k++) | |||
475 | { | |||
476 | snew(gnm[k], STRLEN)(gnm[k]) = save_calloc("gnm[k]", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 476, (4096), sizeof(*(gnm[k]))); | |||
477 | } | |||
478 | for (i = 0; (i < groups->grps[egcENER].nr); i++) | |||
479 | { | |||
480 | ni = groups->grps[egcENER].nm_ind[i]; | |||
481 | for (j = i; (j < groups->grps[egcENER].nr); j++) | |||
482 | { | |||
483 | nj = groups->grps[egcENER].nm_ind[j]; | |||
484 | for (k = kk = 0; (k < egNR); k++) | |||
485 | { | |||
486 | if (md->bEInd[k]) | |||
487 | { | |||
488 | sprintf(gnm[kk], "%s:%s-%s", egrp_nm[k], | |||
489 | *(groups->grpname[ni]), *(groups->grpname[nj])); | |||
490 | kk++; | |||
491 | } | |||
492 | } | |||
493 | md->igrp[n] = get_ebin_space(md->ebin, md->nEc, | |||
494 | (const char **)gnm, unit_energy"kJ/mol"); | |||
495 | n++; | |||
496 | } | |||
497 | } | |||
498 | for (k = 0; (k < md->nEc); k++) | |||
499 | { | |||
500 | sfree(gnm[k])save_free("gnm[k]", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 500, (gnm[k])); | |||
501 | } | |||
502 | sfree(gnm)save_free("gnm", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 502, (gnm)); | |||
503 | ||||
504 | if (n != md->nE) | |||
505 | { | |||
506 | gmx_incons("Number of energy terms wrong")_gmx_error("incons", "Number of energy terms wrong", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 506); | |||
507 | } | |||
508 | } | |||
509 | ||||
510 | md->nTC = groups->grps[egcTC].nr; | |||
511 | md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */ | |||
512 | if (md->bMTTK) | |||
513 | { | |||
514 | md->nTCP = 1; /* assume only one possible coupling system for barostat | |||
515 | for now */ | |||
516 | } | |||
517 | else | |||
518 | { | |||
519 | md->nTCP = 0; | |||
520 | } | |||
521 | if (md->etc == etcNOSEHOOVER) | |||
522 | { | |||
523 | if (md->bNHC_trotter) | |||
524 | { | |||
525 | md->mde_n = 2*md->nNHC*md->nTC; | |||
526 | } | |||
527 | else | |||
528 | { | |||
529 | md->mde_n = 2*md->nTC; | |||
530 | } | |||
531 | if (md->epc == epcMTTK) | |||
532 | { | |||
533 | md->mdeb_n = 2*md->nNHC*md->nTCP; | |||
534 | } | |||
535 | } | |||
536 | else | |||
537 | { | |||
538 | md->mde_n = md->nTC; | |||
539 | md->mdeb_n = 0; | |||
540 | } | |||
541 | ||||
542 | snew(md->tmp_r, md->mde_n)(md->tmp_r) = save_calloc("md->tmp_r", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 542, (md->mde_n), sizeof(*(md->tmp_r))); | |||
543 | snew(md->tmp_v, md->mde_n)(md->tmp_v) = save_calloc("md->tmp_v", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 543, (md->mde_n), sizeof(*(md->tmp_v))); | |||
544 | snew(md->grpnms, md->mde_n)(md->grpnms) = save_calloc("md->grpnms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 544, (md->mde_n), sizeof(*(md->grpnms))); | |||
545 | grpnms = md->grpnms; | |||
546 | ||||
547 | for (i = 0; (i < md->nTC); i++) | |||
548 | { | |||
549 | ni = groups->grps[egcTC].nm_ind[i]; | |||
550 | sprintf(buf, "T-%s", *(groups->grpname[ni])); | |||
551 | grpnms[i] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
552 | } | |||
553 | md->itemp = get_ebin_space(md->ebin, md->nTC, (const char **)grpnms, | |||
554 | unit_temp_K"K"); | |||
555 | ||||
556 | if (md->etc == etcNOSEHOOVER) | |||
557 | { | |||
558 | if (md->bPrintNHChains) | |||
559 | { | |||
560 | if (md->bNHC_trotter) | |||
561 | { | |||
562 | for (i = 0; (i < md->nTC); i++) | |||
563 | { | |||
564 | ni = groups->grps[egcTC].nm_ind[i]; | |||
565 | bufi = *(groups->grpname[ni]); | |||
566 | for (j = 0; (j < md->nNHC); j++) | |||
567 | { | |||
568 | sprintf(buf, "Xi-%d-%s", j, bufi); | |||
569 | grpnms[2*(i*md->nNHC+j)] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
570 | sprintf(buf, "vXi-%d-%s", j, bufi); | |||
571 | grpnms[2*(i*md->nNHC+j)+1] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
572 | } | |||
573 | } | |||
574 | md->itc = get_ebin_space(md->ebin, md->mde_n, | |||
575 | (const char **)grpnms, unit_invtime"1/" "ps"); | |||
576 | if (md->bMTTK) | |||
577 | { | |||
578 | for (i = 0; (i < md->nTCP); i++) | |||
579 | { | |||
580 | bufi = baro_nm[0]; /* All barostat DOF's together for now. */ | |||
581 | for (j = 0; (j < md->nNHC); j++) | |||
582 | { | |||
583 | sprintf(buf, "Xi-%d-%s", j, bufi); | |||
584 | grpnms[2*(i*md->nNHC+j)] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
585 | sprintf(buf, "vXi-%d-%s", j, bufi); | |||
586 | grpnms[2*(i*md->nNHC+j)+1] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
587 | } | |||
588 | } | |||
589 | md->itcb = get_ebin_space(md->ebin, md->mdeb_n, | |||
590 | (const char **)grpnms, unit_invtime"1/" "ps"); | |||
591 | } | |||
592 | } | |||
593 | else | |||
594 | { | |||
595 | for (i = 0; (i < md->nTC); i++) | |||
596 | { | |||
597 | ni = groups->grps[egcTC].nm_ind[i]; | |||
598 | bufi = *(groups->grpname[ni]); | |||
599 | sprintf(buf, "Xi-%s", bufi); | |||
600 | grpnms[2*i] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
601 | sprintf(buf, "vXi-%s", bufi); | |||
602 | grpnms[2*i+1] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
603 | } | |||
604 | md->itc = get_ebin_space(md->ebin, md->mde_n, | |||
605 | (const char **)grpnms, unit_invtime"1/" "ps"); | |||
606 | } | |||
607 | } | |||
608 | } | |||
609 | else if (md->etc == etcBERENDSEN || md->etc == etcYES || | |||
610 | md->etc == etcVRESCALE) | |||
611 | { | |||
612 | for (i = 0; (i < md->nTC); i++) | |||
613 | { | |||
614 | ni = groups->grps[egcTC].nm_ind[i]; | |||
615 | sprintf(buf, "Lamb-%s", *(groups->grpname[ni])); | |||
616 | grpnms[i] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
617 | } | |||
618 | md->itc = get_ebin_space(md->ebin, md->mde_n, (const char **)grpnms, ""); | |||
619 | } | |||
620 | ||||
621 | sfree(grpnms)save_free("grpnms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 621, (grpnms)); | |||
622 | ||||
623 | ||||
624 | md->nU = groups->grps[egcACC].nr; | |||
625 | if (md->nU > 1) | |||
626 | { | |||
627 | snew(grpnms, 3*md->nU)(grpnms) = save_calloc("grpnms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 627, (3*md->nU), sizeof(*(grpnms))); | |||
628 | for (i = 0; (i < md->nU); i++) | |||
629 | { | |||
630 | ni = groups->grps[egcACC].nm_ind[i]; | |||
631 | sprintf(buf, "Ux-%s", *(groups->grpname[ni])); | |||
632 | grpnms[3*i+XX0] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
633 | sprintf(buf, "Uy-%s", *(groups->grpname[ni])); | |||
634 | grpnms[3*i+YY1] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
635 | sprintf(buf, "Uz-%s", *(groups->grpname[ni])); | |||
636 | grpnms[3*i+ZZ2] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
637 | } | |||
638 | md->iu = get_ebin_space(md->ebin, 3*md->nU, (const char **)grpnms, unit_vel"nm" "/" "ps"); | |||
639 | sfree(grpnms)save_free("grpnms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 639, (grpnms)); | |||
640 | } | |||
641 | ||||
642 | if (fp_ene) | |||
643 | { | |||
644 | do_enxnms(fp_ene, &md->ebin->nener, &md->ebin->enm); | |||
645 | } | |||
646 | ||||
647 | md->print_grpnms = NULL((void*)0); | |||
648 | ||||
649 | /* check whether we're going to write dh histograms */ | |||
650 | md->dhc = NULL((void*)0); | |||
651 | if (ir->fepvals->separate_dhdl_file == esepdhdlfileNO) | |||
652 | { | |||
653 | /* Currently dh histograms are only written with dynamics */ | |||
654 | if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) || (ir->eI) == eiBD)) | |||
655 | { | |||
656 | snew(md->dhc, 1)(md->dhc) = save_calloc("md->dhc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 656, (1), sizeof(*(md->dhc))); | |||
657 | ||||
658 | mde_delta_h_coll_init(md->dhc, ir); | |||
659 | } | |||
660 | md->fp_dhdl = NULL((void*)0); | |||
661 | snew(md->dE, ir->fepvals->n_lambda)(md->dE) = save_calloc("md->dE", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 661, (ir->fepvals->n_lambda), sizeof(*(md->dE))); | |||
662 | } | |||
663 | else | |||
664 | { | |||
665 | md->fp_dhdl = fp_dhdl; | |||
666 | snew(md->dE, ir->fepvals->n_lambda)(md->dE) = save_calloc("md->dE", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 666, (ir->fepvals->n_lambda), sizeof(*(md->dE))); | |||
667 | } | |||
668 | if (ir->bSimTemp) | |||
669 | { | |||
670 | int i; | |||
671 | snew(md->temperatures, ir->fepvals->n_lambda)(md->temperatures) = save_calloc("md->temperatures", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 671, (ir->fepvals->n_lambda), sizeof(*(md->temperatures ))); | |||
672 | for (i = 0; i < ir->fepvals->n_lambda; i++) | |||
673 | { | |||
674 | md->temperatures[i] = ir->simtempvals->temperatures[i]; | |||
675 | } | |||
676 | } | |||
677 | return md; | |||
678 | } | |||
679 | ||||
680 | /* print a lambda vector to a string | |||
681 | fep = the inputrec's FEP input data | |||
682 | i = the index of the lambda vector | |||
683 | get_native_lambda = whether to print the native lambda | |||
684 | get_names = whether to print the names rather than the values | |||
685 | str = the pre-allocated string buffer to print to. */ | |||
686 | static void print_lambda_vector(t_lambda *fep, int i, | |||
687 | gmx_bool get_native_lambda, gmx_bool get_names, | |||
688 | char *str) | |||
689 | { | |||
690 | size_t nps = 0, np; | |||
691 | int j, k = 0; | |||
692 | int Nsep = 0; | |||
693 | ||||
694 | for (j = 0; j < efptNR; j++) | |||
695 | { | |||
696 | if (fep->separate_dvdl[j]) | |||
697 | { | |||
698 | Nsep++; | |||
699 | } | |||
700 | } | |||
701 | str[0] = 0; /* reset the string */ | |||
702 | if (Nsep > 1) | |||
703 | { | |||
704 | str += sprintf(str, "("); /* set the opening parenthesis*/ | |||
705 | } | |||
706 | for (j = 0; j < efptNR; j++) | |||
707 | { | |||
708 | if (fep->separate_dvdl[j]) | |||
709 | { | |||
710 | double lam; | |||
711 | if (!get_names) | |||
712 | { | |||
713 | if (get_native_lambda && fep->init_lambda >= 0) | |||
714 | { | |||
715 | str += sprintf(str, "%.4f", fep->init_lambda); | |||
716 | } | |||
717 | else | |||
718 | { | |||
719 | str += sprintf(str, "%.4f", fep->all_lambda[j][i]); | |||
720 | } | |||
721 | } | |||
722 | else | |||
723 | { | |||
724 | str += sprintf(str, "%s", efpt_singular_names[j]); | |||
725 | } | |||
726 | /* print comma for the next item */ | |||
727 | if (k < Nsep-1) | |||
728 | { | |||
729 | str += sprintf(str, ", "); | |||
730 | } | |||
731 | k++; | |||
732 | } | |||
733 | } | |||
734 | if (Nsep > 1) | |||
735 | { | |||
736 | /* and add the closing parenthesis */ | |||
737 | str += sprintf(str, ")"); | |||
738 | } | |||
739 | } | |||
740 | ||||
741 | ||||
742 | extern FILE *open_dhdl(const char *filename, const t_inputrec *ir, | |||
743 | const output_env_t oenv) | |||
744 | { | |||
745 | FILE *fp; | |||
746 | const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda", | |||
747 | *lambdastate = "\\lambda state", *remain = "remaining"; | |||
748 | char title[STRLEN4096], label_x[STRLEN4096], label_y[STRLEN4096]; | |||
749 | int i, np, 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; | |||
753 | char **setname; | |||
754 | char buf[STRLEN4096], lambda_vec_str[STRLEN4096], lambda_name_str[STRLEN4096]; | |||
755 | int bufplace = 0; | |||
756 | ||||
757 | int nsets_dhdl = 0; | |||
758 | int s = 0; | |||
759 | int nsetsextend; | |||
760 | gmx_bool write_pV = FALSE0; | |||
761 | ||||
762 | /* count the number of different lambda terms */ | |||
763 | for (i = 0; i < efptNR; i++) | |||
764 | { | |||
765 | if (fep->separate_dvdl[i]) | |||
766 | { | |||
767 | n_lambda_terms++; | |||
768 | } | |||
769 | } | |||
770 | ||||
771 | if (fep->n_lambda == 0) | |||
772 | { | |||
773 | sprintf(title, "%s", dhdl); | |||
774 | sprintf(label_x, "Time (ps)"); | |||
775 | sprintf(label_y, "%s (%s %s)", | |||
776 | dhdl, unit_energy"kJ/mol", "[\\lambda]\\S-1\\N"); | |||
777 | } | |||
778 | else | |||
779 | { | |||
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"kJ/mol", "[\\8l\\4]\\S-1\\N"); | |||
784 | } | |||
785 | fp = gmx_fio_fopen(filename, "w+"); | |||
786 | xvgr_header(fp, title, label_x, label_y, exvggtXNY, oenv); | |||
787 | ||||
788 | if (!(ir->bSimTemp)) | |||
789 | { | |||
790 | bufplace = sprintf(buf, "T = %g (K) ", | |||
791 | ir->opts.ref_t[0]); | |||
792 | } | |||
793 | if ((ir->efep != efepSLOWGROWTH) && (ir->efep != efepEXPANDED)) | |||
794 | { | |||
795 | if ( (fep->init_lambda >= 0) && (n_lambda_terms == 1 )) | |||
796 | { | |||
797 | /* compatibility output */ | |||
798 | sprintf(&(buf[bufplace]), "%s = %.4f", lambda, fep->init_lambda); | |||
799 | } | |||
800 | else | |||
801 | { | |||
802 | print_lambda_vector(fep, fep->init_fep_state, TRUE1, FALSE0, | |||
803 | lambda_vec_str); | |||
804 | print_lambda_vector(fep, fep->init_fep_state, TRUE1, TRUE1, | |||
805 | lambda_name_str); | |||
806 | sprintf(&(buf[bufplace]), "%s %d: %s = %s", | |||
807 | lambdastate, fep->init_fep_state, | |||
808 | lambda_name_str, lambda_vec_str); | |||
809 | } | |||
810 | } | |||
811 | xvgr_subtitle(fp, buf, oenv); | |||
812 | ||||
813 | ||||
814 | nsets_dhdl = 0; | |||
815 | if (fep->dhdl_derivatives == edhdlderivativesYES) | |||
816 | { | |||
817 | nsets_dhdl = n_lambda_terms; | |||
818 | } | |||
819 | /* count the number of delta_g states */ | |||
820 | nsets_de = fep->lambda_stop_n - fep->lambda_start_n; | |||
821 | ||||
822 | nsets = nsets_dhdl + nsets_de; /* dhdl + fep differences */ | |||
823 | ||||
824 | if (fep->n_lambda > 0 && (expand->elmcmove > elmcmoveNO)) | |||
825 | { | |||
826 | nsets += 1; /*add fep state for expanded ensemble */ | |||
827 | } | |||
828 | ||||
829 | if (fep->bPrintEnergy) | |||
830 | { | |||
831 | nsets += 1; /* add energy to the dhdl as well */ | |||
832 | } | |||
833 | ||||
834 | nsetsextend = nsets; | |||
835 | if ((ir->epc != epcNO) && (fep->n_lambda > 0) && (fep->init_lambda < 0)) | |||
836 | { | |||
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 | |||
841 | dhdl.xvg file) */ | |||
842 | write_pV = TRUE1; | |||
843 | } | |||
844 | snew(setname, nsetsextend)(setname) = save_calloc("setname", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 844, (nsetsextend), sizeof(*(setname))); | |||
845 | ||||
846 | if (expand->elmcmove > elmcmoveNO) | |||
847 | { | |||
848 | /* state for the fep_vals, if we have alchemical sampling */ | |||
849 | sprintf(buf, "%s", "Thermodynamic state"); | |||
850 | setname[s] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
851 | s += 1; | |||
852 | } | |||
853 | ||||
854 | if (fep->bPrintEnergy) | |||
855 | { | |||
856 | sprintf(buf, "%s (%s)", "Energy", unit_energy"kJ/mol"); | |||
857 | setname[s] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
858 | s += 1; | |||
859 | } | |||
860 | ||||
861 | if (fep->dhdl_derivatives == edhdlderivativesYES) | |||
862 | { | |||
863 | for (i = 0; i < efptNR; i++) | |||
864 | { | |||
865 | if (fep->separate_dvdl[i]) | |||
866 | { | |||
867 | ||||
868 | if ( (fep->init_lambda >= 0) && (n_lambda_terms == 1 )) | |||
869 | { | |||
870 | /* compatibility output */ | |||
871 | sprintf(buf, "%s %s %.4f", dhdl, lambda, fep->init_lambda); | |||
872 | } | |||
873 | else | |||
874 | { | |||
875 | double lam = fep->init_lambda; | |||
876 | if (fep->init_lambda < 0) | |||
877 | { | |||
878 | lam = fep->all_lambda[i][fep->init_fep_state]; | |||
879 | } | |||
880 | sprintf(buf, "%s %s = %.4f", dhdl, efpt_singular_names[i], | |||
881 | lam); | |||
882 | } | |||
883 | setname[s] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
884 | s += 1; | |||
885 | } | |||
886 | } | |||
887 | } | |||
888 | ||||
889 | if (fep->n_lambda > 0) | |||
890 | { | |||
891 | /* g_bar has to determine the lambda values used in this simulation | |||
892 | * from this xvg legend. | |||
893 | */ | |||
894 | ||||
895 | if (expand->elmcmove > elmcmoveNO) | |||
896 | { | |||
897 | nsetsbegin = 1; /* for including the expanded ensemble */ | |||
898 | } | |||
899 | else | |||
900 | { | |||
901 | nsetsbegin = 0; | |||
902 | } | |||
903 | ||||
904 | if (fep->bPrintEnergy) | |||
905 | { | |||
906 | nsetsbegin += 1; | |||
907 | } | |||
908 | nsetsbegin += nsets_dhdl; | |||
909 | ||||
910 | for (i = fep->lambda_start_n; i < fep->lambda_stop_n; i++) | |||
911 | { | |||
912 | print_lambda_vector(fep, i, FALSE0, FALSE0, lambda_vec_str); | |||
913 | if ( (fep->init_lambda >= 0) && (n_lambda_terms == 1 )) | |||
914 | { | |||
915 | /* for compatible dhdl.xvg files */ | |||
916 | nps = sprintf(buf, "%s %s %s", deltag, lambda, lambda_vec_str); | |||
917 | } | |||
918 | else | |||
919 | { | |||
920 | nps = sprintf(buf, "%s %s to %s", deltag, lambda, lambda_vec_str); | |||
921 | } | |||
922 | ||||
923 | if (ir->bSimTemp) | |||
924 | { | |||
925 | /* print the temperature for this state if doing simulated annealing */ | |||
926 | sprintf(&buf[nps], "T = %g (%s)", | |||
927 | ir->simtempvals->temperatures[s-(nsetsbegin)], | |||
928 | unit_temp_K"K"); | |||
929 | } | |||
930 | setname[s] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
931 | s++; | |||
932 | } | |||
933 | if (write_pV) | |||
934 | { | |||
935 | np = sprintf(buf, "pV (%s)", unit_energy"kJ/mol"); | |||
936 | setname[nsetsextend-1] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); /* the first entry after | |||
937 | nsets */ | |||
938 | } | |||
939 | ||||
940 | xvgr_legend(fp, nsetsextend, (const char **)setname, oenv); | |||
941 | ||||
942 | for (s = 0; s < nsetsextend; s++) | |||
943 | { | |||
944 | sfree(setname[s])save_free("setname[s]", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 944, (setname[s])); | |||
945 | } | |||
946 | sfree(setname)save_free("setname", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 946, (setname)); | |||
947 | } | |||
948 | ||||
949 | return fp; | |||
950 | } | |||
951 | ||||
952 | static void copy_energy(t_mdebin *md, real e[], real ecpy[]) | |||
953 | { | |||
954 | int i, j; | |||
955 | ||||
956 | for (i = j = 0; (i < F_NRE); i++) | |||
957 | { | |||
958 | if (md->bEner[i]) | |||
959 | { | |||
960 | ecpy[j++] = e[i]; | |||
961 | } | |||
962 | } | |||
963 | if (j != md->f_nre) | |||
964 | { | |||
965 | gmx_incons("Number of energy terms wrong")_gmx_error("incons", "Number of energy terms wrong", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 965); | |||
966 | } | |||
967 | } | |||
968 | ||||
969 | void upd_mdebin(t_mdebin *md, | |||
970 | gmx_bool bDoDHDL, | |||
971 | gmx_bool bSum, | |||
972 | double time, | |||
973 | real tmass, | |||
974 | gmx_enerdata_t *enerd, | |||
975 | t_state *state, | |||
976 | t_lambda *fep, | |||
977 | t_expanded *expand, | |||
978 | matrix box, | |||
979 | tensor svir, | |||
980 | tensor fvir, | |||
981 | tensor vir, | |||
982 | tensor pres, | |||
983 | gmx_ekindata_t *ekind, | |||
984 | rvec mu_tot, | |||
985 | gmx_constr_t constr) | |||
986 | { | |||
987 | int i, j, k, kk, m, n, gid; | |||
988 | real crmsd[2], tmp6[6]; | |||
989 | real bs[NTRICLBOXS((int)(sizeof(tricl_boxs_nm)/sizeof((tricl_boxs_nm)[0])))], vol, dens, pv, enthalpy; | |||
| ||||
990 | real eee[egNR]; | |||
991 | real ecopy[F_NRE]; | |||
992 | double store_dhdl[efptNR]; | |||
993 | real store_energy = 0; | |||
994 | real tmp; | |||
995 | ||||
996 | /* Do NOT use the box in the state variable, but the separate box provided | |||
997 | * as an argument. This is because we sometimes need to write the box from | |||
998 | * the last timestep to match the trajectory frames. | |||
999 | */ | |||
1000 | copy_energy(md, enerd->term, ecopy); | |||
1001 | add_ebin(md->ebin, md->ie, md->f_nre, ecopy, bSum); | |||
1002 | if (md->nCrmsd) | |||
1003 | { | |||
1004 | crmsd[0] = constr_rmsd(constr, FALSE0); | |||
1005 | if (md->nCrmsd > 1) | |||
1006 | { | |||
1007 | crmsd[1] = constr_rmsd(constr, TRUE1); | |||
1008 | } | |||
1009 | add_ebin(md->ebin, md->iconrmsd, md->nCrmsd, crmsd, FALSE0); | |||
1010 | } | |||
1011 | if (md->bDynBox) | |||
1012 | { | |||
1013 | int nboxs; | |||
1014 | if (md->bTricl) | |||
1015 | { | |||
1016 | bs[0] = box[XX0][XX0]; | |||
1017 | bs[1] = box[YY1][YY1]; | |||
1018 | bs[2] = box[ZZ2][ZZ2]; | |||
1019 | bs[3] = box[YY1][XX0]; | |||
1020 | bs[4] = box[ZZ2][XX0]; | |||
1021 | bs[5] = box[ZZ2][YY1]; | |||
1022 | nboxs = NTRICLBOXS((int)(sizeof(tricl_boxs_nm)/sizeof((tricl_boxs_nm)[0]))); | |||
1023 | } | |||
1024 | else | |||
1025 | { | |||
1026 | bs[0] = box[XX0][XX0]; | |||
1027 | bs[1] = box[YY1][YY1]; | |||
1028 | bs[2] = box[ZZ2][ZZ2]; | |||
1029 | nboxs = NBOXS((int)(sizeof(boxs_nm)/sizeof((boxs_nm)[0]))); | |||
1030 | } | |||
1031 | vol = box[XX0][XX0]*box[YY1][YY1]*box[ZZ2][ZZ2]; | |||
1032 | dens = (tmass*AMU(1.6605402e-27))/(vol*NANO(1e-9)*NANO(1e-9)*NANO(1e-9)); | |||
1033 | add_ebin(md->ebin, md->ib, nboxs, bs, bSum); | |||
1034 | add_ebin(md->ebin, md->ivol, 1, &vol, bSum); | |||
1035 | add_ebin(md->ebin, md->idens, 1, &dens, bSum); | |||
1036 | ||||
1037 | if (md->bDiagPres) | |||
1038 | { | |||
1039 | /* This is pV (in kJ/mol). The pressure is the reference pressure, | |||
1040 | not the instantaneous pressure */ | |||
1041 | pv = vol*md->ref_p/PRESFAC(16.6054); | |||
1042 | ||||
1043 | add_ebin(md->ebin, md->ipv, 1, &pv, bSum); | |||
1044 | enthalpy = pv + enerd->term[F_ETOT]; | |||
1045 | add_ebin(md->ebin, md->ienthalpy, 1, &enthalpy, bSum); | |||
1046 | } | |||
1047 | } | |||
1048 | if (md->bConstrVir) | |||
1049 | { | |||
1050 | add_ebin(md->ebin, md->isvir, 9, svir[0], bSum); | |||
1051 | add_ebin(md->ebin, md->ifvir, 9, fvir[0], bSum); | |||
1052 | } | |||
1053 | if (md->bVir) | |||
1054 | { | |||
1055 | add_ebin(md->ebin, md->ivir, 9, vir[0], bSum); | |||
1056 | } | |||
1057 | if (md->bPress) | |||
1058 | { | |||
1059 | add_ebin(md->ebin, md->ipres, 9, pres[0], bSum); | |||
1060 | } | |||
1061 | if (md->bSurft) | |||
1062 | { | |||
1063 | tmp = (pres[ZZ2][ZZ2]-(pres[XX0][XX0]+pres[YY1][YY1])*0.5)*box[ZZ2][ZZ2]; | |||
1064 | add_ebin(md->ebin, md->isurft, 1, &tmp, bSum); | |||
1065 | } | |||
1066 | if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK) | |||
1067 | { | |||
1068 | tmp6[0] = state->boxv[XX0][XX0]; | |||
1069 | tmp6[1] = state->boxv[YY1][YY1]; | |||
1070 | tmp6[2] = state->boxv[ZZ2][ZZ2]; | |||
1071 | tmp6[3] = state->boxv[YY1][XX0]; | |||
1072 | tmp6[4] = state->boxv[ZZ2][XX0]; | |||
1073 | tmp6[5] = state->boxv[ZZ2][YY1]; | |||
1074 | add_ebin(md->ebin, md->ipc, md->bTricl ? 6 : 3, tmp6, bSum); | |||
1075 | } | |||
1076 | if (md->bMu) | |||
1077 | { | |||
1078 | add_ebin(md->ebin, md->imu, 3, mu_tot, bSum); | |||
1079 | } | |||
1080 | if (ekind && ekind->cosacc.cos_accel != 0) | |||
1081 | { | |||
1082 | vol = box[XX0][XX0]*box[YY1][YY1]*box[ZZ2][ZZ2]; | |||
1083 | dens = (tmass*AMU(1.6605402e-27))/(vol*NANO(1e-9)*NANO(1e-9)*NANO(1e-9)); | |||
1084 | add_ebin(md->ebin, md->ivcos, 1, &(ekind->cosacc.vcos), bSum); | |||
1085 | /* 1/viscosity, unit 1/(kg m^-1 s^-1) */ | |||
1086 | tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO(1e-12)) | |||
1087 | *dens*sqr(box[ZZ2][ZZ2]*NANO(1e-9)/(2*M_PI3.14159265358979323846))); | |||
1088 | add_ebin(md->ebin, md->ivisc, 1, &tmp, bSum); | |||
1089 | } | |||
1090 | if (md->nE > 1) | |||
1091 | { | |||
1092 | n = 0; | |||
1093 | for (i = 0; (i < md->nEg); i++) | |||
1094 | { | |||
1095 | for (j = i; (j < md->nEg); j++) | |||
1096 | { | |||
1097 | gid = GID(i, j, md->nEg)((i < j) ? (i*md->nEg+j) : (j*md->nEg+i)); | |||
1098 | for (k = kk = 0; (k < egNR); k++) | |||
1099 | { | |||
1100 | if (md->bEInd[k]) | |||
1101 | { | |||
1102 | eee[kk++] = enerd->grpp.ener[k][gid]; | |||
1103 | } | |||
1104 | } | |||
1105 | add_ebin(md->ebin, md->igrp[n], md->nEc, eee, bSum); | |||
1106 | n++; | |||
1107 | } | |||
1108 | } | |||
1109 | } | |||
1110 | ||||
1111 | if (ekind) | |||
1112 | { | |||
1113 | for (i = 0; (i < md->nTC); i++) | |||
1114 | { | |||
1115 | md->tmp_r[i] = ekind->tcstat[i].T; | |||
1116 | } | |||
1117 | add_ebin(md->ebin, md->itemp, md->nTC, md->tmp_r, bSum); | |||
1118 | ||||
1119 | if (md->etc == etcNOSEHOOVER) | |||
1120 | { | |||
1121 | /* whether to print Nose-Hoover chains: */ | |||
1122 | if (md->bPrintNHChains) | |||
1123 | { | |||
1124 | if (md->bNHC_trotter) | |||
1125 | { | |||
1126 | for (i = 0; (i < md->nTC); i++) | |||
1127 | { | |||
1128 | for (j = 0; j < md->nNHC; j++) | |||
1129 | { | |||
1130 | k = i*md->nNHC+j; | |||
1131 | md->tmp_r[2*k] = state->nosehoover_xi[k]; | |||
1132 | md->tmp_r[2*k+1] = state->nosehoover_vxi[k]; | |||
1133 | } | |||
1134 | } | |||
1135 | add_ebin(md->ebin, md->itc, md->mde_n, md->tmp_r, bSum); | |||
1136 | ||||
1137 | if (md->bMTTK) | |||
1138 | { | |||
1139 | for (i = 0; (i < md->nTCP); i++) | |||
1140 | { | |||
1141 | for (j = 0; j < md->nNHC; j++) | |||
1142 | { | |||
1143 | k = i*md->nNHC+j; | |||
1144 | md->tmp_r[2*k] = state->nhpres_xi[k]; | |||
1145 | md->tmp_r[2*k+1] = state->nhpres_vxi[k]; | |||
1146 | } | |||
1147 | } | |||
1148 | add_ebin(md->ebin, md->itcb, md->mdeb_n, md->tmp_r, bSum); | |||
1149 | } | |||
1150 | } | |||
1151 | else | |||
1152 | { | |||
1153 | for (i = 0; (i < md->nTC); i++) | |||
1154 | { | |||
1155 | md->tmp_r[2*i] = state->nosehoover_xi[i]; | |||
1156 | md->tmp_r[2*i+1] = state->nosehoover_vxi[i]; | |||
1157 | } | |||
1158 | add_ebin(md->ebin, md->itc, md->mde_n, md->tmp_r, bSum); | |||
1159 | } | |||
1160 | } | |||
1161 | } | |||
1162 | else if (md->etc == etcBERENDSEN || md->etc == etcYES || | |||
1163 | md->etc == etcVRESCALE) | |||
1164 | { | |||
1165 | for (i = 0; (i < md->nTC); i++) | |||
1166 | { | |||
1167 | md->tmp_r[i] = ekind->tcstat[i].lambda; | |||
1168 | } | |||
1169 | add_ebin(md->ebin, md->itc, md->nTC, md->tmp_r, bSum); | |||
1170 | } | |||
1171 | } | |||
1172 | ||||
1173 | if (ekind && md->nU > 1) | |||
1174 | { | |||
1175 | for (i = 0; (i < md->nU); i++) | |||
1176 | { | |||
1177 | copy_rvec(ekind->grpstat[i].u, md->tmp_v[i]); | |||
1178 | } | |||
1179 | add_ebin(md->ebin, md->iu, 3*md->nU, md->tmp_v[0], bSum); | |||
1180 | } | |||
1181 | ||||
1182 | ebin_increase_count(md->ebin, bSum); | |||
1183 | ||||
1184 | /* BAR + thermodynamic integration values */ | |||
1185 | if ((md->fp_dhdl || md->dhc) && bDoDHDL) | |||
1186 | { | |||
1187 | for (i = 0; i < enerd->n_lambda-1; i++) | |||
1188 | { | |||
1189 | /* zero for simulated tempering */ | |||
1190 | md->dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0]; | |||
1191 | if (md->temperatures != NULL((void*)0)) | |||
1192 | { | |||
1193 | /* MRS: is this right, given the way we have defined the exchange probabilities? */ | |||
1194 | /* is this even useful to have at all? */ | |||
1195 | md->dE[i] += (md->temperatures[i]/ | |||
1196 | md->temperatures[state->fep_state]-1.0)* | |||
1197 | enerd->term[F_EKIN]; | |||
1198 | } | |||
1199 | } | |||
1200 | ||||
1201 | if (md->fp_dhdl) | |||
1202 | { | |||
1203 | fprintf(md->fp_dhdl, "%.4f", time); | |||
1204 | /* the current free energy state */ | |||
1205 | ||||
1206 | /* print the current state if we are doing expanded ensemble */ | |||
1207 | if (expand->elmcmove > elmcmoveNO) | |||
1208 | { | |||
1209 | fprintf(md->fp_dhdl, " %4d", state->fep_state); | |||
1210 | } | |||
1211 | /* total energy (for if the temperature changes */ | |||
1212 | if (fep->bPrintEnergy) | |||
1213 | { | |||
1214 | store_energy = enerd->term[F_ETOT]; | |||
1215 | fprintf(md->fp_dhdl, " %#.8g", store_energy); | |||
1216 | } | |||
1217 | ||||
1218 | if (fep->dhdl_derivatives == edhdlderivativesYES) | |||
1219 | { | |||
1220 | for (i = 0; i < efptNR; i++) | |||
1221 | { | |||
1222 | if (fep->separate_dvdl[i]) | |||
1223 | { | |||
1224 | /* assumes F_DVDL is first */ | |||
1225 | fprintf(md->fp_dhdl, " %#.8g", enerd->term[F_DVDL+i]); | |||
1226 | } | |||
1227 | } | |||
1228 | } | |||
1229 | for (i = fep->lambda_start_n; i < fep->lambda_stop_n; i++) | |||
1230 | { | |||
1231 | fprintf(md->fp_dhdl, " %#.8g", md->dE[i]); | |||
1232 | } | |||
1233 | if ((md->epc != epcNO) && | |||
1234 | (enerd->n_lambda > 0) && | |||
1235 | (fep->init_lambda < 0)) | |||
1236 | { | |||
1237 | fprintf(md->fp_dhdl, " %#.8g", pv); /* PV term only needed when | |||
| ||||
1238 | there are alternate state | |||
1239 | lambda and we're not in | |||
1240 | compatibility mode */ | |||
1241 | } | |||
1242 | fprintf(md->fp_dhdl, "\n"); | |||
1243 | /* and the binary free energy output */ | |||
1244 | } | |||
1245 | if (md->dhc && bDoDHDL) | |||
1246 | { | |||
1247 | int idhdl = 0; | |||
1248 | for (i = 0; i < efptNR; i++) | |||
1249 | { | |||
1250 | if (fep->separate_dvdl[i]) | |||
1251 | { | |||
1252 | /* assumes F_DVDL is first */ | |||
1253 | store_dhdl[idhdl] = enerd->term[F_DVDL+i]; | |||
1254 | idhdl += 1; | |||
1255 | } | |||
1256 | } | |||
1257 | store_energy = enerd->term[F_ETOT]; | |||
1258 | /* store_dh is dE */ | |||
1259 | mde_delta_h_coll_add_dh(md->dhc, | |||
1260 | (double)state->fep_state, | |||
1261 | store_energy, | |||
1262 | pv, | |||
1263 | store_dhdl, | |||
1264 | md->dE + fep->lambda_start_n, | |||
1265 | time); | |||
1266 | } | |||
1267 | } | |||
1268 | } | |||
1269 | ||||
1270 | ||||
1271 | void upd_mdebin_step(t_mdebin *md) | |||
1272 | { | |||
1273 | ebin_increase_count(md->ebin, FALSE0); | |||
1274 | } | |||
1275 | ||||
1276 | static void npr(FILE *log, int n, char c) | |||
1277 | { | |||
1278 | for (; (n > 0); n--) | |||
1279 | { | |||
1280 | fprintf(log, "%c", c); | |||
1281 | } | |||
1282 | } | |||
1283 | ||||
1284 | static void pprint(FILE *log, const char *s, t_mdebin *md) | |||
1285 | { | |||
1286 | char CHAR = '#'; | |||
1287 | int slen; | |||
1288 | char buf1[22], buf2[22]; | |||
1289 | ||||
1290 | slen = strlen(s); | |||
1291 | fprintf(log, "\t<====== "); | |||
1292 | npr(log, slen, CHAR); | |||
1293 | fprintf(log, " ==>\n"); | |||
1294 | fprintf(log, "\t<==== %s ====>\n", s); | |||
1295 | fprintf(log, "\t<== "); | |||
1296 | npr(log, slen, CHAR); | |||
1297 | fprintf(log, " ======>\n\n"); | |||
1298 | ||||
1299 | fprintf(log, "\tStatistics over %s steps using %s frames\n", | |||
1300 | gmx_step_str(md->ebin->nsteps_sim, buf1), | |||
1301 | gmx_step_str(md->ebin->nsum_sim, buf2)); | |||
1302 | fprintf(log, "\n"); | |||
1303 | } | |||
1304 | ||||
1305 | void print_ebin_header(FILE *log, gmx_int64_t steps, double time, real lambda) | |||
1306 | { | |||
1307 | char buf[22]; | |||
1308 | ||||
1309 | fprintf(log, " %12s %12s %12s\n" | |||
1310 | " %12s %12.5f %12.5f\n\n", | |||
1311 | "Step", "Time", "Lambda", gmx_step_str(steps, buf), time, lambda); | |||
1312 | } | |||
1313 | ||||
1314 | void print_ebin(ener_file_t fp_ene, gmx_bool bEne, gmx_bool bDR, gmx_bool bOR, | |||
1315 | FILE *log, | |||
1316 | gmx_int64_t step, double time, | |||
1317 | int mode, gmx_bool bCompact, | |||
1318 | t_mdebin *md, t_fcdata *fcd, | |||
1319 | gmx_groups_t *groups, t_grpopts *opts) | |||
1320 | { | |||
1321 | /*static char **grpnms=NULL;*/ | |||
1322 | char buf[246]; | |||
1323 | int i, j, n, ni, nj, ndr, nor, b; | |||
1324 | int ndisre = 0; | |||
1325 | real *disre_rm3tav, *disre_rt; | |||
1326 | ||||
1327 | /* these are for the old-style blocks (1 subblock, only reals), because | |||
1328 | there can be only one per ID for these */ | |||
1329 | int nr[enxNR]; | |||
1330 | int id[enxNR]; | |||
1331 | real *block[enxNR]; | |||
1332 | ||||
1333 | /* temporary arrays for the lambda values to write out */ | |||
1334 | double enxlambda_data[2]; | |||
1335 | ||||
1336 | t_enxframe fr; | |||
1337 | ||||
1338 | switch (mode) | |||
1339 | { | |||
1340 | case eprNORMAL: | |||
1341 | init_enxframe(&fr); | |||
1342 | fr.t = time; | |||
1343 | fr.step = step; | |||
1344 | fr.nsteps = md->ebin->nsteps; | |||
1345 | fr.dt = md->delta_t; | |||
1346 | fr.nsum = md->ebin->nsum; | |||
1347 | fr.nre = (bEne) ? md->ebin->nener : 0; | |||
1348 | fr.ener = md->ebin->e; | |||
1349 | ndisre = bDR ? fcd->disres.npair : 0; | |||
1350 | disre_rm3tav = fcd->disres.rm3tav; | |||
1351 | disre_rt = fcd->disres.rt; | |||
1352 | /* Optional additional old-style (real-only) blocks. */ | |||
1353 | for (i = 0; i < enxNR; i++) | |||
1354 | { | |||
1355 | nr[i] = 0; | |||
1356 | } | |||
1357 | if (fcd->orires.nr > 0 && bOR) | |||
1358 | { | |||
1359 | diagonalize_orires_tensors(&(fcd->orires)); | |||
1360 | nr[enxOR] = fcd->orires.nr; | |||
1361 | block[enxOR] = fcd->orires.otav; | |||
1362 | id[enxOR] = enxOR; | |||
1363 | nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ? | |||
1364 | fcd->orires.nr : 0; | |||
1365 | block[enxORI] = fcd->orires.oinsl; | |||
1366 | id[enxORI] = enxORI; | |||
1367 | nr[enxORT] = fcd->orires.nex*12; | |||
1368 | block[enxORT] = fcd->orires.eig; | |||
1369 | id[enxORT] = enxORT; | |||
1370 | } | |||
1371 | ||||
1372 | /* whether we are going to wrte anything out: */ | |||
1373 | if (fr.nre || ndisre || nr[enxOR] || nr[enxORI]) | |||
1374 | { | |||
1375 | ||||
1376 | /* the old-style blocks go first */ | |||
1377 | fr.nblock = 0; | |||
1378 | for (i = 0; i < enxNR; i++) | |||
1379 | { | |||
1380 | if (nr[i] > 0) | |||
1381 | { | |||
1382 | fr.nblock = i + 1; | |||
1383 | } | |||
1384 | } | |||
1385 | add_blocks_enxframe(&fr, fr.nblock); | |||
1386 | for (b = 0; b < fr.nblock; b++) | |||
1387 | { | |||
1388 | add_subblocks_enxblock(&(fr.block[b]), 1); | |||
1389 | fr.block[b].id = id[b]; | |||
1390 | fr.block[b].sub[0].nr = nr[b]; | |||
1391 | #ifndef GMX_DOUBLE | |||
1392 | fr.block[b].sub[0].type = xdr_datatype_float; | |||
1393 | fr.block[b].sub[0].fval = block[b]; | |||
1394 | #else | |||
1395 | fr.block[b].sub[0].type = xdr_datatype_double; | |||
1396 | fr.block[b].sub[0].dval = block[b]; | |||
1397 | #endif | |||
1398 | } | |||
1399 | ||||
1400 | /* check for disre block & fill it. */ | |||
1401 | if (ndisre > 0) | |||
1402 | { | |||
1403 | int db = fr.nblock; | |||
1404 | fr.nblock += 1; | |||
1405 | add_blocks_enxframe(&fr, fr.nblock); | |||
1406 | ||||
1407 | add_subblocks_enxblock(&(fr.block[db]), 2); | |||
1408 | fr.block[db].id = enxDISRE; | |||
1409 | fr.block[db].sub[0].nr = ndisre; | |||
1410 | fr.block[db].sub[1].nr = ndisre; | |||
1411 | #ifndef GMX_DOUBLE | |||
1412 | fr.block[db].sub[0].type = xdr_datatype_float; | |||
1413 | fr.block[db].sub[1].type = xdr_datatype_float; | |||
1414 | fr.block[db].sub[0].fval = disre_rt; | |||
1415 | fr.block[db].sub[1].fval = disre_rm3tav; | |||
1416 | #else | |||
1417 | fr.block[db].sub[0].type = xdr_datatype_double; | |||
1418 | fr.block[db].sub[1].type = xdr_datatype_double; | |||
1419 | fr.block[db].sub[0].dval = disre_rt; | |||
1420 | fr.block[db].sub[1].dval = disre_rm3tav; | |||
1421 | #endif | |||
1422 | } | |||
1423 | /* here we can put new-style blocks */ | |||
1424 | ||||
1425 | /* Free energy perturbation blocks */ | |||
1426 | if (md->dhc) | |||
1427 | { | |||
1428 | mde_delta_h_coll_handle_block(md->dhc, &fr, fr.nblock); | |||
1429 | } | |||
1430 | ||||
1431 | /* we can now free & reset the data in the blocks */ | |||
1432 | if (md->dhc) | |||
1433 | { | |||
1434 | mde_delta_h_coll_reset(md->dhc); | |||
1435 | } | |||
1436 | ||||
1437 | /* do the actual I/O */ | |||
1438 | do_enx(fp_ene, &fr); | |||
1439 | if (fr.nre) | |||
1440 | { | |||
1441 | /* We have stored the sums, so reset the sum history */ | |||
1442 | reset_ebin_sums(md->ebin); | |||
1443 | } | |||
1444 | } | |||
1445 | free_enxframe(&fr); | |||
1446 | break; | |||
1447 | case eprAVER: | |||
1448 | if (log) | |||
1449 | { | |||
1450 | pprint(log, "A V E R A G E S", md); | |||
1451 | } | |||
1452 | break; | |||
1453 | case eprRMS: | |||
1454 | if (log) | |||
1455 | { | |||
1456 | pprint(log, "R M S - F L U C T U A T I O N S", md); | |||
1457 | } | |||
1458 | break; | |||
1459 | default: | |||
1460 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c", 1460, "Invalid print mode (%d)", mode); | |||
1461 | } | |||
1462 | ||||
1463 | if (log) | |||
1464 | { | |||
1465 | for (i = 0; i < opts->ngtc; i++) | |||
1466 | { | |||
1467 | if (opts->annealing[i] != eannNO) | |||
1468 | { | |||
1469 | fprintf(log, "Current ref_t for group %s: %8.1f\n", | |||
1470 | *(groups->grpname[groups->grps[egcTC].nm_ind[i]]), | |||
1471 | opts->ref_t[i]); | |||
1472 | } | |||
1473 | } | |||
1474 | if (mode == eprNORMAL && fcd->orires.nr > 0) | |||
1475 | { | |||
1476 | print_orires_log(log, &(fcd->orires)); | |||
1477 | } | |||
1478 | fprintf(log, " Energies (%s)\n", unit_energy"kJ/mol"); | |||
1479 | pr_ebin(log, md->ebin, md->ie, md->f_nre+md->nCrmsd, 5, mode, TRUE1); | |||
1480 | fprintf(log, "\n"); | |||
1481 | ||||
1482 | if (!bCompact) | |||
1483 | { | |||
1484 | if (md->bDynBox) | |||
1485 | { | |||
1486 | pr_ebin(log, md->ebin, md->ib, md->bTricl ? NTRICLBOXS((int)(sizeof(tricl_boxs_nm)/sizeof((tricl_boxs_nm)[0]))) : NBOXS((int)(sizeof(boxs_nm)/sizeof((boxs_nm)[0]))), 5, | |||
1487 | mode, TRUE1); | |||
1488 | fprintf(log, "\n"); | |||
1489 | } | |||
1490 | if (md->bConstrVir) | |||
1491 | { | |||
1492 | fprintf(log, " Constraint Virial (%s)\n", unit_energy"kJ/mol"); | |||
1493 | pr_ebin(log, md->ebin, md->isvir, 9, 3, mode, FALSE0); | |||
1494 | fprintf(log, "\n"); | |||
1495 | fprintf(log, " Force Virial (%s)\n", unit_energy"kJ/mol"); | |||
1496 | pr_ebin(log, md->ebin, md->ifvir, 9, 3, mode, FALSE0); | |||
1497 | fprintf(log, "\n"); | |||
1498 | } | |||
1499 | if (md->bVir) | |||
1500 | { | |||
1501 | fprintf(log, " Total Virial (%s)\n", unit_energy"kJ/mol"); | |||
1502 | pr_ebin(log, md->ebin, md->ivir, 9, 3, mode, FALSE0); | |||
1503 | fprintf(log, "\n"); | |||
1504 | } | |||
1505 | if (md->bPress) | |||
1506 | { | |||
1507 | fprintf(log, " Pressure (%s)\n", unit_pres_bar"bar"); | |||
1508 | pr_ebin(log, md->ebin, md->ipres, 9, 3, mode, FALSE0); | |||
1509 | fprintf(log, "\n"); | |||
1510 | } | |||
1511 | if (md->bMu) | |||
1512 | { | |||
1513 | fprintf(log, " Total Dipole (%s)\n", unit_dipole_D"D"); | |||
1514 | pr_ebin(log, md->ebin, md->imu, 3, 3, mode, FALSE0); | |||
1515 | fprintf(log, "\n"); | |||
1516 | } | |||
1517 | ||||
1518 | if (md->nE > 1) | |||
1519 | { | |||
1520 | if (md->print_grpnms == NULL((void*)0)) | |||
1521 | { | |||
1522 | snew(md->print_grpnms, md->nE)(md->print_grpnms) = save_calloc("md->print_grpnms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c" , 1522, (md->nE), sizeof(*(md->print_grpnms))); | |||
1523 | n = 0; | |||
1524 | for (i = 0; (i < md->nEg); i++) | |||
1525 | { | |||
1526 | ni = groups->grps[egcENER].nm_ind[i]; | |||
1527 | for (j = i; (j < md->nEg); j++) | |||
1528 | { | |||
1529 | nj = groups->grps[egcENER].nm_ind[j]; | |||
1530 | sprintf(buf, "%s-%s", *(groups->grpname[ni]), | |||
1531 | *(groups->grpname[nj])); | |||
1532 | md->print_grpnms[n++] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); | |||
1533 | } | |||
1534 | } | |||
1535 | } | |||
1536 | sprintf(buf, "Epot (%s)", unit_energy"kJ/mol"); | |||
1537 | fprintf(log, "%15s ", buf); | |||
1538 | for (i = 0; (i < egNR); i++) | |||
1539 | { | |||
1540 | if (md->bEInd[i]) | |||
1541 | { | |||
1542 | fprintf(log, "%12s ", egrp_nm[i]); | |||
1543 | } | |||
1544 | } | |||
1545 | fprintf(log, "\n"); | |||
1546 | for (i = 0; (i < md->nE); i++) | |||
1547 | { | |||
1548 | fprintf(log, "%15s", md->print_grpnms[i]); | |||
1549 | pr_ebin(log, md->ebin, md->igrp[i], md->nEc, md->nEc, mode, | |||
1550 | FALSE0); | |||
1551 | } | |||
1552 | fprintf(log, "\n"); | |||
1553 | } | |||
1554 | if (md->nTC > 1) | |||
1555 | { | |||
1556 | pr_ebin(log, md->ebin, md->itemp, md->nTC, 4, mode, TRUE1); | |||
1557 | fprintf(log, "\n"); | |||
1558 | } | |||
1559 | if (md->nU > 1) | |||
1560 | { | |||
1561 | fprintf(log, "%15s %12s %12s %12s\n", | |||
1562 | "Group", "Ux", "Uy", "Uz"); | |||
1563 | for (i = 0; (i < md->nU); i++) | |||
1564 | { | |||
1565 | ni = groups->grps[egcACC].nm_ind[i]; | |||
1566 | fprintf(log, "%15s", *groups->grpname[ni]); | |||
1567 | pr_ebin(log, md->ebin, md->iu+3*i, 3, 3, mode, FALSE0); | |||
1568 | } | |||
1569 | fprintf(log, "\n"); | |||
1570 | } | |||
1571 | } | |||
1572 | } | |||
1573 | ||||
1574 | } | |||
1575 | ||||
1576 | void update_energyhistory(energyhistory_t * enerhist, t_mdebin * mdebin) | |||
1577 | { | |||
1578 | int i; | |||
1579 | ||||
1580 | enerhist->nsteps = mdebin->ebin->nsteps; | |||
1581 | enerhist->nsum = mdebin->ebin->nsum; | |||
1582 | enerhist->nsteps_sim = mdebin->ebin->nsteps_sim; | |||
1583 | enerhist->nsum_sim = mdebin->ebin->nsum_sim; | |||
1584 | enerhist->nener = mdebin->ebin->nener; | |||
1585 | ||||
1586 | if (mdebin->ebin->nsum > 0) | |||
1587 | { | |||
1588 | /* Check if we need to allocate first */ | |||
1589 | if (enerhist->ener_ave == NULL((void*)0)) | |||
1590 | { | |||
1591 | snew(enerhist->ener_ave, enerhist->nener)(enerhist->ener_ave) = save_calloc("enerhist->ener_ave" , "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c", 1591 , (enerhist->nener), sizeof(*(enerhist->ener_ave))); | |||
1592 | snew(enerhist->ener_sum, enerhist->nener)(enerhist->ener_sum) = save_calloc("enerhist->ener_sum" , "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c", 1592 , (enerhist->nener), sizeof(*(enerhist->ener_sum))); | |||
1593 | } | |||
1594 | ||||
1595 | for (i = 0; i < enerhist->nener; i++) | |||
1596 | { | |||
1597 | enerhist->ener_ave[i] = mdebin->ebin->e[i].eav; | |||
1598 | enerhist->ener_sum[i] = mdebin->ebin->e[i].esum; | |||
1599 | } | |||
1600 | } | |||
1601 | ||||
1602 | if (mdebin->ebin->nsum_sim > 0) | |||
1603 | { | |||
1604 | /* Check if we need to allocate first */ | |||
1605 | if (enerhist->ener_sum_sim == NULL((void*)0)) | |||
1606 | { | |||
1607 | snew(enerhist->ener_sum_sim, enerhist->nener)(enerhist->ener_sum_sim) = save_calloc("enerhist->ener_sum_sim" , "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c", 1607 , (enerhist->nener), sizeof(*(enerhist->ener_sum_sim))); | |||
1608 | } | |||
1609 | ||||
1610 | for (i = 0; i < enerhist->nener; i++) | |||
1611 | { | |||
1612 | enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum; | |||
1613 | } | |||
1614 | } | |||
1615 | if (mdebin->dhc) | |||
1616 | { | |||
1617 | mde_delta_h_coll_update_energyhistory(mdebin->dhc, enerhist); | |||
1618 | } | |||
1619 | } | |||
1620 | ||||
1621 | void restore_energyhistory_from_state(t_mdebin * mdebin, | |||
1622 | energyhistory_t * enerhist) | |||
1623 | { | |||
1624 | int i; | |||
1625 | ||||
1626 | if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) && | |||
1627 | mdebin->ebin->nener != enerhist->nener) | |||
1628 | { | |||
1629 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/mdebin.c", 1629, "Mismatch between number of energies in run input (%d) and checkpoint file (%d).", | |||
1630 | mdebin->ebin->nener, enerhist->nener); | |||
1631 | } | |||
1632 | ||||
1633 | mdebin->ebin->nsteps = enerhist->nsteps; | |||
1634 | mdebin->ebin->nsum = enerhist->nsum; | |||
1635 | mdebin->ebin->nsteps_sim = enerhist->nsteps_sim; | |||
1636 | mdebin->ebin->nsum_sim = enerhist->nsum_sim; | |||
1637 | ||||
1638 | for (i = 0; i < mdebin->ebin->nener; i++) | |||
1639 | { | |||
1640 | mdebin->ebin->e[i].eav = | |||
1641 | (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0); | |||
1642 | mdebin->ebin->e[i].esum = | |||
1643 | (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0); | |||
1644 | mdebin->ebin->e_sim[i].esum = | |||
1645 | (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0); | |||
1646 | } | |||
1647 | if (mdebin->dhc) | |||
1648 | { | |||
1649 | mde_delta_h_coll_restore_energyhistory(mdebin->dhc, enerhist); | |||
1650 | } | |||
1651 | } |