File: | gromacs/mdlib/mdebin.c |
Location: | line 737, column 9 |
Description: | Value stored to 'str' is never read |
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, ")"); |
Value stored to 'str' is never read | |
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 | } |