?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Disable instruction fusion on Power8
[gromacs.git] / src / gromacs / fileio / tpxio.cpp
blobb4526ddb818dff0137941c8354cf90ea7a865306
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,2015,2016,2017,2018, 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 #include "gmxpre.h"
38
39 /* This file is completely threadsafe - keep it that way! */
40
41 #include "tpxio.h"
42
43 #include <cstdio>
44 #include <cstdlib>
45 #include <cstring>
46
47 #include <algorithm>
48 #include <vector>
49
50 #include "gromacs/fileio/filetypes.h"
51 #include "gromacs/fileio/gmxfio.h"
52 #include "gromacs/fileio/gmxfio-xdr.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/mdtypes/awh-history.h"
56 #include "gromacs/mdtypes/awh-params.h"
57 #include "gromacs/mdtypes/inputrec.h"
58 #include "gromacs/mdtypes/md_enums.h"
59 #include "gromacs/mdtypes/pull-params.h"
60 #include "gromacs/mdtypes/state.h"
61 #include "gromacs/pbcutil/boxutilities.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/topology/block.h"
64 #include "gromacs/topology/ifunc.h"
65 #include "gromacs/topology/mtop_util.h"
66 #include "gromacs/topology/symtab.h"
67 #include "gromacs/topology/topology.h"
68 #include "gromacs/utility/arraysize.h"
69 #include "gromacs/utility/baseversion.h"
70 #include "gromacs/utility/cstringutil.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/futil.h"
73 #include "gromacs/utility/gmxassert.h"
74 #include "gromacs/utility/keyvaluetreebuilder.h"
75 #include "gromacs/utility/keyvaluetreeserializer.h"
76 #include "gromacs/utility/smalloc.h"
77 #include "gromacs/utility/snprintf.h"
78 #include "gromacs/utility/txtdump.h"
79
80 #define TPX_TAG_RELEASE "release"
81
82 /*! \brief Tag string for the file format written to run input files
83 * written by this version of the code.
84 *
85 * Change this if you want to change the run input format in a feature
86 * branch. This ensures that there will not be different run input
87 * formats around which cannot be distinguished, while not causing
88 * problems rebasing the feature branch onto upstream changes. When
89 * merging with mainstream GROMACS, set this tag string back to
90 * TPX_TAG_RELEASE, and instead add an element to tpxv.
91 */
92 static const char *tpx_tag = TPX_TAG_RELEASE;
93
94 /*! \brief Enum of values that describe the contents of a tpr file
95 * whose format matches a version number
96 *
97 * The enum helps the code be more self-documenting and ensure merges
98 * do not silently resolve when two patches make the same bump. When
99 * adding new functionality, add a new element just above tpxv_Count
100 * in this enumeration, and write code below that does the right thing
101 * according to the value of file_version.
102 */
103 enum tpxv {
104 tpxv_ComputationalElectrophysiology = 96, /**< support for ion/water position swaps (computational electrophysiology) */
105 tpxv_Use64BitRandomSeed, /**< change ld_seed from int to gmx_int64_t */
106 tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, /**< potentials for supporting coarse-grained force fields */
107 tpxv_InteractiveMolecularDynamics, /**< interactive molecular dynamics (IMD) */
108 tpxv_RemoveObsoleteParameters1, /**< remove optimize_fft, dihre_fc, nstcheckpoint */
109 tpxv_PullCoordTypeGeom, /**< add pull type and geometry per group and flat-bottom */
110 tpxv_PullGeomDirRel, /**< add pull geometry direction-relative */
111 tpxv_IntermolecularBondeds, /**< permit inter-molecular bonded interactions in the topology */
112 tpxv_CompElWithSwapLayerOffset, /**< added parameters for improved CompEl setups */
113 tpxv_CompElPolyatomicIonsAndMultipleIonTypes, /**< CompEl now can handle polyatomic ions and more than two types of ions */
114 tpxv_RemoveAdress, /**< removed support for AdResS */
115 tpxv_PullCoordNGroup, /**< add ngroup to pull coord */
116 tpxv_RemoveTwinRange, /**< removed support for twin-range interactions */
117 tpxv_ReplacePullPrintCOM12, /**< Replaced print-com-1, 2 with pull-print-com */
118 tpxv_PullExternalPotential, /**< Added pull type external potential */
119 tpxv_GenericParamsForElectricField, /**< Introduced KeyValueTree and moved electric field parameters */
120 tpxv_AcceleratedWeightHistogram, /**< sampling with accelerated weight histogram method (AWH) */
121 tpxv_Count /**< the total number of tpxv versions */
122 };
123
124 /*! \brief Version number of the file format written to run input
125 * files by this version of the code.
126 *
127 * The tpx_version increases whenever the file format in the main
128 * development branch changes, due to an extension of the tpxv enum above.
129 * Backward compatibility for reading old run input files is maintained
130 * by checking this version number against that of the file and then using
131 * the correct code path.
132 *
133 * When developing a feature branch that needs to change the run input
134 * file format, change tpx_tag instead. */
135 static const int tpx_version = tpxv_Count - 1;
136
137
138 /* This number should only be increased when you edit the TOPOLOGY section
139 * or the HEADER of the tpx format.
140 * This way we can maintain forward compatibility too for all analysis tools
141 * and/or external programs that only need to know the atom/residue names,
142 * charges, and bond connectivity.
143 *
144 * It first appeared in tpx version 26, when I also moved the inputrecord
145 * to the end of the tpx file, so we can just skip it if we only
146 * want the topology.
147 *
148 * In particular, it must be increased when adding new elements to
149 * ftupd, so that old code can read new .tpr files.
150 */
151 static const int tpx_generation = 26;
152
153 /* This number should be the most recent backwards incompatible version
154 * I.e., if this number is 9, we cannot read tpx version 9 with this code.
155 */
156 static const int tpx_incompatible_version = 30; // GMX3.2 has version 31
157
158
159
160 /* Struct used to maintain tpx compatibility when function types are added */
161 typedef struct {
162 int fvnr; /* file version number in which the function type first appeared */
163 int ftype; /* function type */
164 } t_ftupd;
165
166 /*
167 * The entries should be ordered in:
168 * 1. ascending function type number
169 * 2. ascending file version number
170 */
171 static const t_ftupd ftupd[] = {
172 { 34, F_FENEBONDS },
173 { 43, F_TABBONDS },
174 { 43, F_TABBONDSNC },
175 { 70, F_RESTRBONDS },
176 { tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_RESTRANGLES },
177 { 76, F_LINEAR_ANGLES },
178 { 34, F_QUARTIC_ANGLES },
179 { 43, F_TABANGLES },
180 { tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_RESTRDIHS },
181 { tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_CBTDIHS },
182 { 43, F_TABDIHS },
183 { 65, F_CMAP },
184 { 60, F_GB12 },
185 { 61, F_GB13 },
186 { 61, F_GB14 },
187 { 72, F_GBPOL },
188 { 72, F_NPSOLVATION },
189 { 41, F_LJC14_Q },
190 { 41, F_LJC_PAIRS_NB },
191 { 32, F_BHAM_LR_NOLONGERUSED },
192 { 32, F_RF_EXCL },
193 { 32, F_COUL_RECIP },
194 { 93, F_LJ_RECIP },
195 { 46, F_DPD },
196 { 36, F_THOLE_POL },
197 { 90, F_FBPOSRES },
198 { 49, F_VSITE4FDN },
199 { 50, F_VSITEN },
200 { 46, F_COM_PULL },
201 { 46, F_ECONSERVED },
202 { 69, F_VTEMP_NOLONGERUSED},
203 { 66, F_PDISPCORR },
204 { 54, F_DVDL_CONSTR },
205 { 76, F_ANHARM_POL },
206 { 79, F_DVDL_COUL },
207 { 79, F_DVDL_VDW, },
208 { 79, F_DVDL_BONDED, },
209 { 79, F_DVDL_RESTRAINT },
210 { 79, F_DVDL_TEMPERATURE },
211 };
212 #define NFTUPD asize(ftupd)
213
214 /* Needed for backward compatibility */
215 #define MAXNODES 256
216
217 /**************************************************************
218 *
219 * Now the higer level routines that do io of the structures and arrays
220 *
221 **************************************************************/
222 static void do_pullgrp_tpx_pre95(t_fileio *fio,
223 t_pull_group *pgrp,
224 t_pull_coord *pcrd,
225 gmx_bool bRead,
226 int file_version)
227 {
228 rvec tmp;
229
230 gmx_fio_do_int(fio, pgrp->nat);
231 if (bRead)
232 {
233 snew(pgrp->ind, pgrp->nat);
234 }
235 gmx_fio_ndo_int(fio, pgrp->ind, pgrp->nat);
236 gmx_fio_do_int(fio, pgrp->nweight);
237 if (bRead)
238 {
239 snew(pgrp->weight, pgrp->nweight);
240 }
241 gmx_fio_ndo_real(fio, pgrp->weight, pgrp->nweight);
242 gmx_fio_do_int(fio, pgrp->pbcatom);
243 gmx_fio_do_rvec(fio, pcrd->vec);
244 clear_rvec(pcrd->origin);
245 gmx_fio_do_rvec(fio, tmp);
246 pcrd->init = tmp[0];
247 gmx_fio_do_real(fio, pcrd->rate);
248 gmx_fio_do_real(fio, pcrd->k);
249 if (file_version >= 56)
250 {
251 gmx_fio_do_real(fio, pcrd->kB);
252 }
253 else
254 {
255 pcrd->kB = pcrd->k;
256 }
257 }
258
259 static void do_pull_group(t_fileio *fio, t_pull_group *pgrp, gmx_bool bRead)
260 {
261 gmx_fio_do_int(fio, pgrp->nat);
262 if (bRead)
263 {
264 snew(pgrp->ind, pgrp->nat);
265 }
266 gmx_fio_ndo_int(fio, pgrp->ind, pgrp->nat);
267 gmx_fio_do_int(fio, pgrp->nweight);
268 if (bRead)
269 {
270 snew(pgrp->weight, pgrp->nweight);
271 }
272 gmx_fio_ndo_real(fio, pgrp->weight, pgrp->nweight);
273 gmx_fio_do_int(fio, pgrp->pbcatom);
274 }
275
276 static void do_pull_coord(t_fileio *fio, t_pull_coord *pcrd,
277 gmx_bool bRead, int file_version,
278 int ePullOld, int eGeomOld, ivec dimOld)
279 {
280 if (file_version >= tpxv_PullCoordNGroup)
281 {
282 gmx_fio_do_int(fio, pcrd->eType);
283 if (file_version >= tpxv_PullExternalPotential)
284 {
285 if (pcrd->eType == epullEXTERNAL)
286 {
287 if (bRead)
288 {
289 char buf[STRLEN];
290
291 gmx_fio_do_string(fio, buf);
292 pcrd->externalPotentialProvider = gmx_strdup(buf);
293 }
294 else
295 {
296 gmx_fio_do_string(fio, pcrd->externalPotentialProvider);
297 }
298 }
299 else
300 {
301 pcrd->externalPotentialProvider = nullptr;
302 }
303 }
304 else
305 {
306 if (bRead)
307 {
308 pcrd->externalPotentialProvider = nullptr;
309 }
310 }
311 /* Note that we try to support adding new geometries without
312 * changing the tpx version. This requires checks when printing the
313 * geometry string and a check and fatal_error in init_pull.
314 */
315 gmx_fio_do_int(fio, pcrd->eGeom);
316 gmx_fio_do_int(fio, pcrd->ngroup);
317 if (pcrd->ngroup <= c_pullCoordNgroupMax)
318 {
319 gmx_fio_ndo_int(fio, pcrd->group, pcrd->ngroup);
320 }
321 else
322 {
323 /* More groups in file than supported, this must be a new geometry
324 * that is not supported by our current code. Since we will not
325 * use the groups for this coord (checks in the pull and WHAM code
326 * ensure this), we can ignore the groups and set ngroup=0.
327 */
328 int *dum;
329 snew(dum, pcrd->ngroup);
330 gmx_fio_ndo_int(fio, dum, pcrd->ngroup);
331 sfree(dum);
332
333 pcrd->ngroup = 0;
334 }
335 gmx_fio_do_ivec(fio, pcrd->dim);
336 }
337 else
338 {
339 pcrd->ngroup = 2;
340 gmx_fio_do_int(fio, pcrd->group[0]);
341 gmx_fio_do_int(fio, pcrd->group[1]);
342 if (file_version >= tpxv_PullCoordTypeGeom)
343 {
344 pcrd->ngroup = (pcrd->eGeom == epullgDIRRELATIVE ? 4 : 2);
345 gmx_fio_do_int(fio, pcrd->eType);
346 gmx_fio_do_int(fio, pcrd->eGeom);
347 if (pcrd->ngroup == 4)
348 {
349 gmx_fio_do_int(fio, pcrd->group[2]);
350 gmx_fio_do_int(fio, pcrd->group[3]);
351 }
352 gmx_fio_do_ivec(fio, pcrd->dim);
353 }
354 else
355 {
356 pcrd->eType = ePullOld;
357 pcrd->eGeom = eGeomOld;
358 copy_ivec(dimOld, pcrd->dim);
359 }
360 }
361 gmx_fio_do_rvec(fio, pcrd->origin);
362 gmx_fio_do_rvec(fio, pcrd->vec);
363 if (file_version >= tpxv_PullCoordTypeGeom)
364 {
365 gmx_fio_do_gmx_bool(fio, pcrd->bStart);
366 }
367 else
368 {
369 /* This parameter is only printed, but not actually used by mdrun */
370 pcrd->bStart = FALSE;
371 }
372 gmx_fio_do_real(fio, pcrd->init);
373 gmx_fio_do_real(fio, pcrd->rate);
374 gmx_fio_do_real(fio, pcrd->k);
375 gmx_fio_do_real(fio, pcrd->kB);
376 }
377
378 static void do_expandedvals(t_fileio *fio, t_expanded *expand, t_lambda *fepvals, gmx_bool bRead, int file_version)
379 {
380 int n_lambda = fepvals->n_lambda;
381
382 /* reset the lambda calculation window */
383 fepvals->lambda_start_n = 0;
384 fepvals->lambda_stop_n = n_lambda;
385 if (file_version >= 79)
386 {
387 if (n_lambda > 0)
388 {
389 if (bRead)
390 {
391 snew(expand->init_lambda_weights, n_lambda);
392 }
393 gmx_fio_ndo_real(fio, expand->init_lambda_weights, n_lambda);
394 gmx_fio_do_gmx_bool(fio, expand->bInit_weights);
395 }
396
397 gmx_fio_do_int(fio, expand->nstexpanded);
398 gmx_fio_do_int(fio, expand->elmcmove);
399 gmx_fio_do_int(fio, expand->elamstats);
400 gmx_fio_do_int(fio, expand->lmc_repeats);
401 gmx_fio_do_int(fio, expand->gibbsdeltalam);
402 gmx_fio_do_int(fio, expand->lmc_forced_nstart);
403 gmx_fio_do_int(fio, expand->lmc_seed);
404 gmx_fio_do_real(fio, expand->mc_temp);
405 gmx_fio_do_int(fio, expand->bSymmetrizedTMatrix);
406 gmx_fio_do_int(fio, expand->nstTij);
407 gmx_fio_do_int(fio, expand->minvarmin);
408 gmx_fio_do_int(fio, expand->c_range);
409 gmx_fio_do_real(fio, expand->wl_scale);
410 gmx_fio_do_real(fio, expand->wl_ratio);
411 gmx_fio_do_real(fio, expand->init_wl_delta);
412 gmx_fio_do_gmx_bool(fio, expand->bWLoneovert);
413 gmx_fio_do_int(fio, expand->elmceq);
414 gmx_fio_do_int(fio, expand->equil_steps);
415 gmx_fio_do_int(fio, expand->equil_samples);
416 gmx_fio_do_int(fio, expand->equil_n_at_lam);
417 gmx_fio_do_real(fio, expand->equil_wl_delta);
418 gmx_fio_do_real(fio, expand->equil_ratio);
419 }
420 }
421
422 static void do_simtempvals(t_fileio *fio, t_simtemp *simtemp, int n_lambda, gmx_bool bRead,
423 int file_version)
424 {
425 if (file_version >= 79)
426 {
427 gmx_fio_do_int(fio, simtemp->eSimTempScale);
428 gmx_fio_do_real(fio, simtemp->simtemp_high);
429 gmx_fio_do_real(fio, simtemp->simtemp_low);
430 if (n_lambda > 0)
431 {
432 if (bRead)
433 {
434 snew(simtemp->temperatures, n_lambda);
435 }
436 gmx_fio_ndo_real(fio, simtemp->temperatures, n_lambda);
437 }
438 }
439 }
440
441 static void do_imd(t_fileio *fio, t_IMD *imd, gmx_bool bRead)
442 {
443 gmx_fio_do_int(fio, imd->nat);
444 if (bRead)
445 {
446 snew(imd->ind, imd->nat);
447 }
448 gmx_fio_ndo_int(fio, imd->ind, imd->nat);
449 }
450
451 static void do_fepvals(t_fileio *fio, t_lambda *fepvals, gmx_bool bRead, int file_version)
452 {
453 /* i is defined in the ndo_double macro; use g to iterate. */
454 int g;
455 real rdum;
456
457 /* free energy values */
458
459 if (file_version >= 79)
460 {
461 gmx_fio_do_int(fio, fepvals->init_fep_state);
462 gmx_fio_do_double(fio, fepvals->init_lambda);
463 gmx_fio_do_double(fio, fepvals->delta_lambda);
464 }
465 else if (file_version >= 59)
466 {
467 gmx_fio_do_double(fio, fepvals->init_lambda);
468 gmx_fio_do_double(fio, fepvals->delta_lambda);
469 }
470 else
471 {
472 gmx_fio_do_real(fio, rdum);
473 fepvals->init_lambda = rdum;
474 gmx_fio_do_real(fio, rdum);
475 fepvals->delta_lambda = rdum;
476 }
477 if (file_version >= 79)
478 {
479 gmx_fio_do_int(fio, fepvals->n_lambda);
480 if (bRead)
481 {
482 snew(fepvals->all_lambda, efptNR);
483 }
484 for (g = 0; g < efptNR; g++)
485 {
486 if (fepvals->n_lambda > 0)
487 {
488 if (bRead)
489 {
490 snew(fepvals->all_lambda[g], fepvals->n_lambda);
491 }
492 gmx_fio_ndo_double(fio, fepvals->all_lambda[g], fepvals->n_lambda);
493 gmx_fio_ndo_int(fio, fepvals->separate_dvdl, efptNR);
494 }
495 else if (fepvals->init_lambda >= 0)
496 {
497 fepvals->separate_dvdl[efptFEP] = TRUE;
498 }
499 }
500 }
501 else if (file_version >= 64)
502 {
503 gmx_fio_do_int(fio, fepvals->n_lambda);
504 if (bRead)
505 {
506 int g;
507
508 snew(fepvals->all_lambda, efptNR);
509 /* still allocate the all_lambda array's contents. */
510 for (g = 0; g < efptNR; g++)
511 {
512 if (fepvals->n_lambda > 0)
513 {
514 snew(fepvals->all_lambda[g], fepvals->n_lambda);
515 }
516 }
517 }
518 gmx_fio_ndo_double(fio, fepvals->all_lambda[efptFEP],
519 fepvals->n_lambda);
520 if (fepvals->init_lambda >= 0)
521 {
522 int g, h;
523
524 fepvals->separate_dvdl[efptFEP] = TRUE;
525
526 if (bRead)
527 {
528 /* copy the contents of the efptFEP lambda component to all
529 the other components */
530 for (g = 0; g < efptNR; g++)
531 {
532 for (h = 0; h < fepvals->n_lambda; h++)
533 {
534 if (g != efptFEP)
535 {
536 fepvals->all_lambda[g][h] =
537 fepvals->all_lambda[efptFEP][h];
538 }
539 }
540 }
541 }
542 }
543 }
544 else
545 {
546 fepvals->n_lambda = 0;
547 fepvals->all_lambda = nullptr;
548 if (fepvals->init_lambda >= 0)
549 {
550 fepvals->separate_dvdl[efptFEP] = TRUE;
551 }
552 }
553 gmx_fio_do_real(fio, fepvals->sc_alpha);
554 if (file_version >= 38)
555 {
556 gmx_fio_do_int(fio, fepvals->sc_power);
557 }
558 else
559 {
560 fepvals->sc_power = 2;
561 }
562 if (file_version >= 79)
563 {
564 gmx_fio_do_real(fio, fepvals->sc_r_power);
565 }
566 else
567 {
568 fepvals->sc_r_power = 6.0;
569 }
570 gmx_fio_do_real(fio, fepvals->sc_sigma);
571 if (bRead)
572 {
573 if (file_version >= 71)
574 {
575 fepvals->sc_sigma_min = fepvals->sc_sigma;
576 }
577 else
578 {
579 fepvals->sc_sigma_min = 0;
580 }
581 }
582 if (file_version >= 79)
583 {
584 gmx_fio_do_int(fio, fepvals->bScCoul);
585 }
586 else
587 {
588 fepvals->bScCoul = TRUE;
589 }
590 if (file_version >= 64)
591 {
592 gmx_fio_do_int(fio, fepvals->nstdhdl);
593 }
594 else
595 {
596 fepvals->nstdhdl = 1;
597 }
598
599 if (file_version >= 73)
600 {
601 gmx_fio_do_int(fio, fepvals->separate_dhdl_file);
602 gmx_fio_do_int(fio, fepvals->dhdl_derivatives);
603 }
604 else
605 {
606 fepvals->separate_dhdl_file = esepdhdlfileYES;
607 fepvals->dhdl_derivatives = edhdlderivativesYES;
608 }
609 if (file_version >= 71)
610 {
611 gmx_fio_do_int(fio, fepvals->dh_hist_size);
612 gmx_fio_do_double(fio, fepvals->dh_hist_spacing);
613 }
614 else
615 {
616 fepvals->dh_hist_size = 0;
617 fepvals->dh_hist_spacing = 0.1;
618 }
619 if (file_version >= 79)
620 {
621 gmx_fio_do_int(fio, fepvals->edHdLPrintEnergy);
622 }
623 else
624 {
625 fepvals->edHdLPrintEnergy = edHdLPrintEnergyNO;
626 }
627
628 /* handle lambda_neighbors */
629 if ((file_version >= 83 && file_version < 90) || file_version >= 92)
630 {
631 gmx_fio_do_int(fio, fepvals->lambda_neighbors);
632 if ( (fepvals->lambda_neighbors >= 0) && (fepvals->init_fep_state >= 0) &&
633 (fepvals->init_lambda < 0) )
634 {
635 fepvals->lambda_start_n = (fepvals->init_fep_state -
636 fepvals->lambda_neighbors);
637 fepvals->lambda_stop_n = (fepvals->init_fep_state +
638 fepvals->lambda_neighbors + 1);
639 if (fepvals->lambda_start_n < 0)
640 {
641 fepvals->lambda_start_n = 0;;
642 }
643 if (fepvals->lambda_stop_n >= fepvals->n_lambda)
644 {
645 fepvals->lambda_stop_n = fepvals->n_lambda;
646 }
647 }
648 else
649 {
650 fepvals->lambda_start_n = 0;
651 fepvals->lambda_stop_n = fepvals->n_lambda;
652 }
653 }
654 else
655 {
656 fepvals->lambda_start_n = 0;
657 fepvals->lambda_stop_n = fepvals->n_lambda;
658 }
659 }
660
661 static void do_awhBias(t_fileio *fio, gmx::AwhBiasParams *awhBiasParams, gmx_bool bRead,
662 int gmx_unused file_version)
663 {
664 gmx_fio_do_int(fio, awhBiasParams->eTarget);
665 gmx_fio_do_double(fio, awhBiasParams->targetBetaScaling);
666 gmx_fio_do_double(fio, awhBiasParams->targetCutoff);
667 gmx_fio_do_int(fio, awhBiasParams->eGrowth);
668 gmx_fio_do_int(fio, awhBiasParams->bUserData);
669 gmx_fio_do_double(fio, awhBiasParams->errorInitial);
670 gmx_fio_do_int(fio, awhBiasParams->ndim);
671 gmx_fio_do_int(fio, awhBiasParams->shareGroup);
672 gmx_fio_do_gmx_bool(fio, awhBiasParams->equilibrateHistogram);
673
674 if (bRead)
675 {
676 snew(awhBiasParams->dimParams, awhBiasParams->ndim);
677 }
678
679 for (int d = 0; d < awhBiasParams->ndim; d++)
680 {
681 gmx::AwhDimParams *dimParams = &awhBiasParams->dimParams[d];
682
683 gmx_fio_do_int(fio, dimParams->eCoordProvider);
684 gmx_fio_do_int(fio, dimParams->coordIndex);
685 gmx_fio_do_double(fio, dimParams->origin);
686 gmx_fio_do_double(fio, dimParams->end);
687 gmx_fio_do_double(fio, dimParams->period);
688 gmx_fio_do_double(fio, dimParams->forceConstant);
689 gmx_fio_do_double(fio, dimParams->diffusion);
690 gmx_fio_do_double(fio, dimParams->coordValueInit);
691 gmx_fio_do_double(fio, dimParams->coverDiameter);
692 }
693 }
694
695 static void do_awh(t_fileio *fio, gmx::AwhParams *awhParams, gmx_bool bRead,
696 int gmx_unused file_version)
697 {
698 gmx_fio_do_int(fio, awhParams->numBias);
699 gmx_fio_do_int(fio, awhParams->nstOut);
700 gmx_fio_do_int64(fio, awhParams->seed);
701 gmx_fio_do_int(fio, awhParams->nstSampleCoord);
702 gmx_fio_do_int(fio, awhParams->numSamplesUpdateFreeEnergy);
703 gmx_fio_do_int(fio, awhParams->ePotential);
704 gmx_fio_do_gmx_bool(fio, awhParams->shareBiasMultisim);
705
706 if (awhParams->numBias > 0)
707 {
708 if (bRead)
709 {
710 snew(awhParams->awhBiasParams, awhParams->numBias);
711 }
712
713 for (int k = 0; k < awhParams->numBias; k++)
714 {
715 do_awhBias(fio, &awhParams->awhBiasParams[k], bRead, file_version);
716 }
717 }
718 }
719
720 static void do_pull(t_fileio *fio, pull_params_t *pull, gmx_bool bRead,
721 int file_version, int ePullOld)
722 {
723 int eGeomOld = -1;
724 ivec dimOld;
725 int g;
726
727 if (file_version >= 95)
728 {
729 gmx_fio_do_int(fio, pull->ngroup);
730 }
731 gmx_fio_do_int(fio, pull->ncoord);
732 if (file_version < 95)
733 {
734 pull->ngroup = pull->ncoord + 1;
735 }
736 if (file_version < tpxv_PullCoordTypeGeom)
737 {
738 real dum;
739
740 gmx_fio_do_int(fio, eGeomOld);
741 gmx_fio_do_ivec(fio, dimOld);
742 /* The inner cylinder radius, now removed */
743 gmx_fio_do_real(fio, dum);
744 }
745 gmx_fio_do_real(fio, pull->cylinder_r);
746 gmx_fio_do_real(fio, pull->constr_tol);
747 if (file_version >= 95)
748 {
749 gmx_fio_do_int(fio, pull->bPrintCOM);
750 /* With file_version < 95 this value is set below */
751 }
752 if (file_version >= tpxv_ReplacePullPrintCOM12)
753 {
754 gmx_fio_do_int(fio, pull->bPrintRefValue);
755 gmx_fio_do_int(fio, pull->bPrintComp);
756 }
757 else if (file_version >= tpxv_PullCoordTypeGeom)
758 {
759 int idum;
760 gmx_fio_do_int(fio, idum); /* used to be bPrintCOM2 */
761 gmx_fio_do_int(fio, pull->bPrintRefValue);
762 gmx_fio_do_int(fio, pull->bPrintComp);
763 }
764 else
765 {
766 pull->bPrintRefValue = FALSE;
767 pull->bPrintComp = TRUE;
768 }
769 gmx_fio_do_int(fio, pull->nstxout);
770 gmx_fio_do_int(fio, pull->nstfout);
771 if (bRead)
772 {
773 snew(pull->group, pull->ngroup);
774 snew(pull->coord, pull->ncoord);
775 }
776 if (file_version < 95)
777 {
778 /* epullgPOS for position pulling, before epullgDIRPBC was removed */
779 if (eGeomOld == epullgDIRPBC)
780 {
781 gmx_fatal(FARGS, "pull-geometry=position is no longer supported");
782 }
783 if (eGeomOld > epullgDIRPBC)
784 {
785 eGeomOld -= 1;
786 }
787
788 for (g = 0; g < pull->ngroup; g++)
789 {
790 /* We read and ignore a pull coordinate for group 0 */
791 do_pullgrp_tpx_pre95(fio, &pull->group[g], &pull->coord[std::max(g-1, 0)],
792 bRead, file_version);
793 if (g > 0)
794 {
795 pull->coord[g-1].group[0] = 0;
796 pull->coord[g-1].group[1] = g;
797 }
798 }
799
800 pull->bPrintCOM = (pull->group[0].nat > 0);
801 }
802 else
803 {
804 for (g = 0; g < pull->ngroup; g++)
805 {
806 do_pull_group(fio, &pull->group[g], bRead);
807 }
808 for (g = 0; g < pull->ncoord; g++)
809 {
810 do_pull_coord(fio, &pull->coord[g],
811 bRead, file_version, ePullOld, eGeomOld, dimOld);
812 }
813 }
814 }
815
816
817 static void do_rotgrp(t_fileio *fio, t_rotgrp *rotg, gmx_bool bRead)
818 {
819 gmx_fio_do_int(fio, rotg->eType);
820 gmx_fio_do_int(fio, rotg->bMassW);
821 gmx_fio_do_int(fio, rotg->nat);
822 if (bRead)
823 {
824 snew(rotg->ind, rotg->nat);
825 }
826 gmx_fio_ndo_int(fio, rotg->ind, rotg->nat);
827 if (bRead)
828 {
829 snew(rotg->x_ref, rotg->nat);
830 }
831 gmx_fio_ndo_rvec(fio, rotg->x_ref, rotg->nat);
832 gmx_fio_do_rvec(fio, rotg->vec);
833 gmx_fio_do_rvec(fio, rotg->pivot);
834 gmx_fio_do_real(fio, rotg->rate);
835 gmx_fio_do_real(fio, rotg->k);
836 gmx_fio_do_real(fio, rotg->slab_dist);
837 gmx_fio_do_real(fio, rotg->min_gaussian);
838 gmx_fio_do_real(fio, rotg->eps);
839 gmx_fio_do_int(fio, rotg->eFittype);
840 gmx_fio_do_int(fio, rotg->PotAngle_nstep);
841 gmx_fio_do_real(fio, rotg->PotAngle_step);
842 }
843
844 static void do_rot(t_fileio *fio, t_rot *rot, gmx_bool bRead)
845 {
846 int g;
847
848 gmx_fio_do_int(fio, rot->ngrp);
849 gmx_fio_do_int(fio, rot->nstrout);
850 gmx_fio_do_int(fio, rot->nstsout);
851 if (bRead)
852 {
853 snew(rot->grp, rot->ngrp);
854 }
855 for (g = 0; g < rot->ngrp; g++)
856 {
857 do_rotgrp(fio, &rot->grp[g], bRead);
858 }
859 }
860
861
862 static void do_swapgroup(t_fileio *fio, t_swapGroup *g, gmx_bool bRead)
863 {
864
865 /* Name of the group or molecule */
866 if (bRead)
867 {
868 char buf[STRLEN];
869
870 gmx_fio_do_string(fio, buf);
871 g->molname = gmx_strdup(buf);
872 }
873 else
874 {
875 gmx_fio_do_string(fio, g->molname);
876 }
877
878 /* Number of atoms in the group */
879 gmx_fio_do_int(fio, g->nat);
880
881 /* The group's atom indices */
882 if (bRead)
883 {
884 snew(g->ind, g->nat);
885 }
886 gmx_fio_ndo_int(fio, g->ind, g->nat);
887
888 /* Requested counts for compartments A and B */
889 gmx_fio_ndo_int(fio, g->nmolReq, eCompNR);
890 }
891
892 static void do_swapcoords_tpx(t_fileio *fio, t_swapcoords *swap, gmx_bool bRead, int file_version)
893 {
894 /* Enums for better readability of the code */
895 enum {
896 eCompA = 0, eCompB
897 };
898 enum {
899 eChannel0 = 0, eChannel1
900 };
901
902
903 if (file_version >= tpxv_CompElPolyatomicIonsAndMultipleIonTypes)
904 {
905 /* The total number of swap groups is the sum of the fixed groups
906 * (split0, split1, solvent), and the user-defined groups (2+ types of ions)
907 */
908 gmx_fio_do_int(fio, swap->ngrp);
909 if (bRead)
910 {
911 snew(swap->grp, swap->ngrp);
912 }
913 for (int ig = 0; ig < swap->ngrp; ig++)
914 {
915 do_swapgroup(fio, &swap->grp[ig], bRead);
916 }
917 gmx_fio_do_int(fio, swap->massw_split[eChannel0]);
918 gmx_fio_do_int(fio, swap->massw_split[eChannel1]);
919 gmx_fio_do_int(fio, swap->nstswap);
920 gmx_fio_do_int(fio, swap->nAverage);
921 gmx_fio_do_real(fio, swap->threshold);
922 gmx_fio_do_real(fio, swap->cyl0r);
923 gmx_fio_do_real(fio, swap->cyl0u);
924 gmx_fio_do_real(fio, swap->cyl0l);
925 gmx_fio_do_real(fio, swap->cyl1r);
926 gmx_fio_do_real(fio, swap->cyl1u);
927 gmx_fio_do_real(fio, swap->cyl1l);
928 }
929 else
930 {
931 /*** Support reading older CompEl .tpr files ***/
932
933 /* In the original CompEl .tpr files, we always have 5 groups: */
934 swap->ngrp = 5;
935 snew(swap->grp, swap->ngrp);
936
937 swap->grp[eGrpSplit0 ].molname = gmx_strdup("split0" ); // group 0: split0
938 swap->grp[eGrpSplit1 ].molname = gmx_strdup("split1" ); // group 1: split1
939 swap->grp[eGrpSolvent].molname = gmx_strdup("solvent"); // group 2: solvent
940 swap->grp[3 ].molname = gmx_strdup("anions" ); // group 3: anions
941 swap->grp[4 ].molname = gmx_strdup("cations"); // group 4: cations
942
943 gmx_fio_do_int(fio, swap->grp[3].nat);
944 gmx_fio_do_int(fio, swap->grp[eGrpSolvent].nat);
945 gmx_fio_do_int(fio, swap->grp[eGrpSplit0].nat);
946 gmx_fio_do_int(fio, swap->massw_split[eChannel0]);
947 gmx_fio_do_int(fio, swap->grp[eGrpSplit1].nat);
948 gmx_fio_do_int(fio, swap->massw_split[eChannel1]);
949 gmx_fio_do_int(fio, swap->nstswap);
950 gmx_fio_do_int(fio, swap->nAverage);
951 gmx_fio_do_real(fio, swap->threshold);
952 gmx_fio_do_real(fio, swap->cyl0r);
953 gmx_fio_do_real(fio, swap->cyl0u);
954 gmx_fio_do_real(fio, swap->cyl0l);
955 gmx_fio_do_real(fio, swap->cyl1r);
956 gmx_fio_do_real(fio, swap->cyl1u);
957 gmx_fio_do_real(fio, swap->cyl1l);
958
959 // The order[] array keeps compatibility with older .tpr files
960 // by reading in the groups in the classic order
961 {
962 const int order[4] = {3, eGrpSolvent, eGrpSplit0, eGrpSplit1 };
963
964 for (int ig = 0; ig < 4; ig++)
965 {
966 int g = order[ig];
967 snew(swap->grp[g].ind, swap->grp[g].nat);
968 gmx_fio_ndo_int(fio, swap->grp[g].ind, swap->grp[g].nat);
969 }
970 }
971
972 for (int j = eCompA; j <= eCompB; j++)
973 {
974 gmx_fio_do_int(fio, swap->grp[3].nmolReq[j]); // group 3 = anions
975 gmx_fio_do_int(fio, swap->grp[4].nmolReq[j]); // group 4 = cations
976 }
977 } /* End support reading older CompEl .tpr files */
978
979 if (file_version >= tpxv_CompElWithSwapLayerOffset)
980 {
981 gmx_fio_do_real(fio, swap->bulkOffset[eCompA]);
982 gmx_fio_do_real(fio, swap->bulkOffset[eCompB]);
983 }
984
985 }
986
987 static void do_legacy_efield(t_fileio *fio, gmx::KeyValueTreeObjectBuilder *root)
988 {
989 const char *const dimName[] = { "x", "y", "z" };
990
991 auto appliedForcesObj = root->addObject("applied-forces");
992 auto efieldObj = appliedForcesObj.addObject("electric-field");
993 // The content of the tpr file for this feature has
994 // been the same since gromacs 4.0 that was used for
995 // developing.
996 for (int j = 0; j < DIM; ++j)
997 {
998 int n, nt;
999 gmx_fio_do_int(fio, n);
1000 gmx_fio_do_int(fio, nt);
1001 std::vector<real> aa(n+1), phi(nt+1), at(nt+1), phit(nt+1);
1002 gmx_fio_ndo_real(fio, aa.data(), n);
1003 gmx_fio_ndo_real(fio, phi.data(), n);
1004 gmx_fio_ndo_real(fio, at.data(), nt);
1005 gmx_fio_ndo_real(fio, phit.data(), nt);
1006 if (n > 0)
1007 {
1008 if (n > 1 || nt > 1)
1009 {
1010 gmx_fatal(FARGS, "Can not handle tpr files with more than one electric field term per direction.");
1011 }
1012 auto dimObj = efieldObj.addObject(dimName[j]);
1013 dimObj.addValue<real>("E0", aa[0]);
1014 dimObj.addValue<real>("omega", at[0]);
1015 dimObj.addValue<real>("t0", phi[0]);
1016 dimObj.addValue<real>("sigma", phit[0]);
1017 }
1018 }
1019 }
1020
1021
1022 static void do_inputrec(t_fileio *fio, t_inputrec *ir, gmx_bool bRead,
1023 int file_version, real *fudgeQQ)
1024 {
1025 int i, j, k, idum = 0;
1026 real rdum, bd_temp;
1027 gmx_bool bdum = 0;
1028
1029 if (file_version != tpx_version)
1030 {
1031 /* Give a warning about features that are not accessible */
1032 fprintf(stderr, "Note: file tpx version %d, software tpx version %d\n",
1033 file_version, tpx_version);
1034 }
1035
1036 if (file_version == 0)
1037 {
1038 return;
1039 }
1040
1041 gmx::KeyValueTreeBuilder paramsBuilder;
1042 gmx::KeyValueTreeObjectBuilder paramsObj = paramsBuilder.rootObject();
1043
1044 /* Basic inputrec stuff */
1045 gmx_fio_do_int(fio, ir->eI);
1046 if (file_version >= 62)
1047 {
1048 gmx_fio_do_int64(fio, ir->nsteps);
1049 }
1050 else
1051 {
1052 gmx_fio_do_int(fio, idum);
1053 ir->nsteps = idum;
1054 }
1055
1056 if (file_version >= 62)
1057 {
1058 gmx_fio_do_int64(fio, ir->init_step);
1059 }
1060 else
1061 {
1062 gmx_fio_do_int(fio, idum);
1063 ir->init_step = idum;
1064 }
1065
1066 if (file_version >= 58)
1067 {
1068 gmx_fio_do_int(fio, ir->simulation_part);
1069 }
1070 else
1071 {
1072 ir->simulation_part = 1;
1073 }
1074
1075 if (file_version >= 67)
1076 {
1077 gmx_fio_do_int(fio, ir->nstcalcenergy);
1078 }
1079 else
1080 {
1081 ir->nstcalcenergy = 1;
1082 }
1083 if (file_version < 53)
1084 {
1085 /* The pbc info has been moved out of do_inputrec,
1086 * since we always want it, also without reading the inputrec.
1087 */
1088 gmx_fio_do_int(fio, ir->ePBC);
1089 if (file_version >= 45)
1090 {
1091 gmx_fio_do_int(fio, ir->bPeriodicMols);
1092 }
1093 else
1094 {
1095 if (ir->ePBC == 2)
1096 {
1097 ir->ePBC = epbcXYZ;
1098 ir->bPeriodicMols = TRUE;
1099 }
1100 else
1101 {
1102 ir->bPeriodicMols = FALSE;
1103 }
1104 }
1105 }
1106 if (file_version >= 81)
1107 {
1108 gmx_fio_do_int(fio, ir->cutoff_scheme);
1109 if (file_version < 94)
1110 {
1111 ir->cutoff_scheme = 1 - ir->cutoff_scheme;
1112 }
1113 }
1114 else
1115 {
1116 ir->cutoff_scheme = ecutsGROUP;
1117 }
1118 gmx_fio_do_int(fio, ir->ns_type);
1119 gmx_fio_do_int(fio, ir->nstlist);
1120 gmx_fio_do_int(fio, idum); /* used to be ndelta; not used anymore */
1121 if (file_version < 41)
1122 {
1123 gmx_fio_do_int(fio, idum);
1124 gmx_fio_do_int(fio, idum);
1125 }
1126 if (file_version >= 45)
1127 {
1128 gmx_fio_do_real(fio, ir->rtpi);
1129 }
1130 else
1131 {
1132 ir->rtpi = 0.05;
1133 }
1134 gmx_fio_do_int(fio, ir->nstcomm);
1135 if (file_version > 34)
1136 {
1137 gmx_fio_do_int(fio, ir->comm_mode);
1138 }
1139 else if (ir->nstcomm < 0)
1140 {
1141 ir->comm_mode = ecmANGULAR;
1142 }
1143 else
1144 {
1145 ir->comm_mode = ecmLINEAR;
1146 }
1147 ir->nstcomm = abs(ir->nstcomm);
1148
1149 /* ignore nstcheckpoint */
1150 if (file_version < tpxv_RemoveObsoleteParameters1)
1151 {
1152 gmx_fio_do_int(fio, idum);
1153 }
1154
1155 gmx_fio_do_int(fio, ir->nstcgsteep);
1156
1157 gmx_fio_do_int(fio, ir->nbfgscorr);
1158
1159 gmx_fio_do_int(fio, ir->nstlog);
1160 gmx_fio_do_int(fio, ir->nstxout);
1161 gmx_fio_do_int(fio, ir->nstvout);
1162 gmx_fio_do_int(fio, ir->nstfout);
1163 gmx_fio_do_int(fio, ir->nstenergy);
1164 gmx_fio_do_int(fio, ir->nstxout_compressed);
1165 if (file_version >= 59)
1166 {
1167 gmx_fio_do_double(fio, ir->init_t);
1168 gmx_fio_do_double(fio, ir->delta_t);
1169 }
1170 else
1171 {
1172 gmx_fio_do_real(fio, rdum);
1173 ir->init_t = rdum;
1174 gmx_fio_do_real(fio, rdum);
1175 ir->delta_t = rdum;
1176 }
1177 gmx_fio_do_real(fio, ir->x_compression_precision);
1178 if (file_version >= 81)
1179 {
1180 gmx_fio_do_real(fio, ir->verletbuf_tol);
1181 }
1182 else
1183 {
1184 ir->verletbuf_tol = 0;
1185 }
1186 gmx_fio_do_real(fio, ir->rlist);
1187 if (file_version >= 67 && file_version < tpxv_RemoveTwinRange)
1188 {
1189 if (bRead)
1190 {
1191 // Reading such a file version could invoke the twin-range
1192 // scheme, about which mdrun should give a fatal error.
1193 real dummy_rlistlong = -1;
1194 gmx_fio_do_real(fio, dummy_rlistlong);
1195
1196 if (ir->rlist > 0 && (dummy_rlistlong == 0 || dummy_rlistlong > ir->rlist))
1197 {
1198 // Get mdrun to issue an error (regardless of
1199 // ir->cutoff_scheme).
1200 ir->useTwinRange = true;
1201 }
1202 else
1203 {
1204 // grompp used to set rlistlong actively. Users were
1205 // probably also confused and set rlistlong == rlist.
1206 // However, in all remaining cases, it is safe to let
1207 // mdrun proceed normally.
1208 ir->useTwinRange = false;
1209 }
1210 }
1211 }
1212 else
1213 {
1214 // No need to read or write anything
1215 ir->useTwinRange = false;
1216 }
1217 if (file_version >= 82 && file_version != 90)
1218 {
1219 // Multiple time-stepping is no longer enabled, but the old
1220 // support required the twin-range scheme, for which mdrun
1221 // already emits a fatal error.
1222 int dummy_nstcalclr = -1;
1223 gmx_fio_do_int(fio, dummy_nstcalclr);
1224 }
1225 gmx_fio_do_int(fio, ir->coulombtype);
1226 if (file_version < 32 && ir->coulombtype == eelRF)
1227 {
1228 ir->coulombtype = eelRF_NEC_UNSUPPORTED;
1229 }
1230 if (file_version >= 81)
1231 {
1232 gmx_fio_do_int(fio, ir->coulomb_modifier);
1233 }
1234 else
1235 {
1236 ir->coulomb_modifier = (ir->cutoff_scheme == ecutsVERLET ? eintmodPOTSHIFT : eintmodNONE);
1237 }
1238 gmx_fio_do_real(fio, ir->rcoulomb_switch);
1239 gmx_fio_do_real(fio, ir->rcoulomb);
1240 gmx_fio_do_int(fio, ir->vdwtype);
1241 if (file_version >= 81)
1242 {
1243 gmx_fio_do_int(fio, ir->vdw_modifier);
1244 }
1245 else
1246 {
1247 ir->vdw_modifier = (ir->cutoff_scheme == ecutsVERLET ? eintmodPOTSHIFT : eintmodNONE);
1248 }
1249 gmx_fio_do_real(fio, ir->rvdw_switch);
1250 gmx_fio_do_real(fio, ir->rvdw);
1251 gmx_fio_do_int(fio, ir->eDispCorr);
1252 gmx_fio_do_real(fio, ir->epsilon_r);
1253 if (file_version >= 37)
1254 {
1255 gmx_fio_do_real(fio, ir->epsilon_rf);
1256 }
1257 else
1258 {
1259 if (EEL_RF(ir->coulombtype))
1260 {
1261 ir->epsilon_rf = ir->epsilon_r;
1262 ir->epsilon_r = 1.0;
1263 }
1264 else
1265 {
1266 ir->epsilon_rf = 1.0;
1267 }
1268 }
1269 gmx_fio_do_real(fio, ir->tabext);
1270
1271 gmx_fio_do_int(fio, ir->gb_algorithm);
1272 gmx_fio_do_int(fio, ir->nstgbradii);
1273 gmx_fio_do_real(fio, ir->rgbradii);
1274 gmx_fio_do_real(fio, ir->gb_saltconc);
1275 gmx_fio_do_int(fio, ir->implicit_solvent);
1276 if (file_version >= 55)
1277 {
1278 gmx_fio_do_real(fio, ir->gb_epsilon_solvent);
1279 gmx_fio_do_real(fio, ir->gb_obc_alpha);
1280 gmx_fio_do_real(fio, ir->gb_obc_beta);
1281 gmx_fio_do_real(fio, ir->gb_obc_gamma);
1282 if (file_version >= 60)
1283 {
1284 gmx_fio_do_real(fio, ir->gb_dielectric_offset);
1285 gmx_fio_do_int(fio, ir->sa_algorithm);
1286 }
1287 else
1288 {
1289 ir->gb_dielectric_offset = 0.009;
1290 ir->sa_algorithm = esaAPPROX;
1291 }
1292 gmx_fio_do_real(fio, ir->sa_surface_tension);
1293
1294 /* Override sa_surface_tension if it is not changed in the mpd-file */
1295 if (ir->sa_surface_tension < 0)
1296 {
1297 if (ir->gb_algorithm == egbSTILL)
1298 {
1299 ir->sa_surface_tension = 0.0049 * 100 * CAL2JOULE;
1300 }
1301 else if (ir->gb_algorithm == egbHCT || ir->gb_algorithm == egbOBC)
1302 {
1303 ir->sa_surface_tension = 0.0054 * 100 * CAL2JOULE;
1304 }
1305 }
1306
1307 }
1308 else
1309 {
1310 /* Better use sensible values than insane (0.0) ones... */
1311 ir->gb_epsilon_solvent = 80;
1312 ir->gb_obc_alpha = 1.0;
1313 ir->gb_obc_beta = 0.8;
1314 ir->gb_obc_gamma = 4.85;
1315 ir->sa_surface_tension = 2.092;
1316 }
1317
1318
1319 if (file_version >= 81)
1320 {
1321 gmx_fio_do_real(fio, ir->fourier_spacing);
1322 }
1323 else
1324 {
1325 ir->fourier_spacing = 0.0;
1326 }
1327 gmx_fio_do_int(fio, ir->nkx);
1328 gmx_fio_do_int(fio, ir->nky);
1329 gmx_fio_do_int(fio, ir->nkz);
1330 gmx_fio_do_int(fio, ir->pme_order);
1331 gmx_fio_do_real(fio, ir->ewald_rtol);
1332
1333 if (file_version >= 93)
1334 {
1335 gmx_fio_do_real(fio, ir->ewald_rtol_lj);
1336 }
1337 else
1338 {
1339 ir->ewald_rtol_lj = ir->ewald_rtol;
1340 }
1341 gmx_fio_do_int(fio, ir->ewald_geometry);
1342 gmx_fio_do_real(fio, ir->epsilon_surface);
1343
1344 /* ignore bOptFFT */
1345 if (file_version < tpxv_RemoveObsoleteParameters1)
1346 {
1347 gmx_fio_do_gmx_bool(fio, bdum);
1348 }
1349
1350 if (file_version >= 93)
1351 {
1352 gmx_fio_do_int(fio, ir->ljpme_combination_rule);
1353 }
1354 gmx_fio_do_gmx_bool(fio, ir->bContinuation);
1355 gmx_fio_do_int(fio, ir->etc);
1356 /* before version 18, ir->etc was a gmx_bool (ir->btc),
1357 * but the values 0 and 1 still mean no and
1358 * berendsen temperature coupling, respectively.
1359 */
1360 if (file_version >= 79)
1361 {
1362 gmx_fio_do_gmx_bool(fio, ir->bPrintNHChains);
1363 }
1364 if (file_version >= 71)
1365 {
1366 gmx_fio_do_int(fio, ir->nsttcouple);
1367 }
1368 else
1369 {
1370 ir->nsttcouple = ir->nstcalcenergy;
1371 }
1372 gmx_fio_do_int(fio, ir->epc);
1373 gmx_fio_do_int(fio, ir->epct);
1374 if (file_version >= 71)
1375 {
1376 gmx_fio_do_int(fio, ir->nstpcouple);
1377 }
1378 else
1379 {
1380 ir->nstpcouple = ir->nstcalcenergy;
1381 }
1382 gmx_fio_do_real(fio, ir->tau_p);
1383 gmx_fio_do_rvec(fio, ir->ref_p[XX]);
1384 gmx_fio_do_rvec(fio, ir->ref_p[YY]);
1385 gmx_fio_do_rvec(fio, ir->ref_p[ZZ]);
1386 gmx_fio_do_rvec(fio, ir->compress[XX]);
1387 gmx_fio_do_rvec(fio, ir->compress[YY]);
1388 gmx_fio_do_rvec(fio, ir->compress[ZZ]);
1389 if (file_version >= 47)
1390 {
1391 gmx_fio_do_int(fio, ir->refcoord_scaling);
1392 gmx_fio_do_rvec(fio, ir->posres_com);
1393 gmx_fio_do_rvec(fio, ir->posres_comB);
1394 }
1395 else
1396 {
1397 ir->refcoord_scaling = erscNO;
1398 clear_rvec(ir->posres_com);
1399 clear_rvec(ir->posres_comB);
1400 }
1401 if (file_version < 79)
1402 {
1403 gmx_fio_do_int(fio, ir->andersen_seed);
1404 }
1405 else
1406 {
1407 ir->andersen_seed = 0;
1408 }
1409
1410 if (file_version < 37)
1411 {
1412 gmx_fio_do_real(fio, rdum);
1413 }
1414
1415 gmx_fio_do_real(fio, ir->shake_tol);
1416 if (file_version < 54)
1417 {
1418 // cppcheck-suppress redundantPointerOp
1419 gmx_fio_do_real(fio, *fudgeQQ);
1420 }
1421
1422 gmx_fio_do_int(fio, ir->efep);
1423 do_fepvals(fio, ir->fepvals, bRead, file_version);
1424
1425 if (file_version >= 79)
1426 {
1427 gmx_fio_do_gmx_bool(fio, ir->bSimTemp);
1428 if (ir->bSimTemp)
1429 {
1430 ir->bSimTemp = TRUE;
1431 }
1432 }
1433 else
1434 {
1435 ir->bSimTemp = FALSE;
1436 }
1437 if (ir->bSimTemp)
1438 {
1439 do_simtempvals(fio, ir->simtempvals, ir->fepvals->n_lambda, bRead, file_version);
1440 }
1441
1442 if (file_version >= 79)
1443 {
1444 gmx_fio_do_gmx_bool(fio, ir->bExpanded);
1445 if (ir->bExpanded)
1446 {
1447 ir->bExpanded = TRUE;
1448 }
1449 else
1450 {
1451 ir->bExpanded = FALSE;
1452 }
1453 }
1454 if (ir->bExpanded)
1455 {
1456 do_expandedvals(fio, ir->expandedvals, ir->fepvals, bRead, file_version);
1457 }
1458 if (file_version >= 57)
1459 {
1460 gmx_fio_do_int(fio, ir->eDisre);
1461 }
1462 gmx_fio_do_int(fio, ir->eDisreWeighting);
1463 gmx_fio_do_gmx_bool(fio, ir->bDisreMixed);
1464 gmx_fio_do_real(fio, ir->dr_fc);
1465 gmx_fio_do_real(fio, ir->dr_tau);
1466 gmx_fio_do_int(fio, ir->nstdisreout);
1467 gmx_fio_do_real(fio, ir->orires_fc);
1468 gmx_fio_do_real(fio, ir->orires_tau);
1469 gmx_fio_do_int(fio, ir->nstorireout);
1470
1471 /* ignore dihre_fc */
1472 if (file_version < 79)
1473 {
1474 gmx_fio_do_real(fio, rdum);
1475 if (file_version < 56)
1476 {
1477 gmx_fio_do_real(fio, rdum);
1478 gmx_fio_do_int(fio, idum);
1479 }
1480 }
1481
1482 gmx_fio_do_real(fio, ir->em_stepsize);
1483 gmx_fio_do_real(fio, ir->em_tol);
1484 gmx_fio_do_gmx_bool(fio, ir->bShakeSOR);
1485 gmx_fio_do_int(fio, ir->niter);
1486 gmx_fio_do_real(fio, ir->fc_stepsize);
1487 gmx_fio_do_int(fio, ir->eConstrAlg);
1488 gmx_fio_do_int(fio, ir->nProjOrder);
1489 gmx_fio_do_real(fio, ir->LincsWarnAngle);
1490 gmx_fio_do_int(fio, ir->nLincsIter);
1491 if (file_version < 33)
1492 {
1493 gmx_fio_do_real(fio, bd_temp);
1494 }
1495 gmx_fio_do_real(fio, ir->bd_fric);
1496 if (file_version >= tpxv_Use64BitRandomSeed)
1497 {
1498 gmx_fio_do_int64(fio, ir->ld_seed);
1499 }
1500 else
1501 {
1502 gmx_fio_do_int(fio, idum);
1503 ir->ld_seed = idum;
1504 }
1505 if (file_version >= 33)
1506 {
1507 for (i = 0; i < DIM; i++)
1508 {
1509 gmx_fio_do_rvec(fio, ir->deform[i]);
1510 }
1511 }
1512 else
1513 {
1514 for (i = 0; i < DIM; i++)
1515 {
1516 clear_rvec(ir->deform[i]);
1517 }
1518 }
1519 gmx_fio_do_real(fio, ir->cos_accel);
1520
1521 gmx_fio_do_int(fio, ir->userint1);
1522 gmx_fio_do_int(fio, ir->userint2);
1523 gmx_fio_do_int(fio, ir->userint3);
1524 gmx_fio_do_int(fio, ir->userint4);
1525 gmx_fio_do_real(fio, ir->userreal1);
1526 gmx_fio_do_real(fio, ir->userreal2);
1527 gmx_fio_do_real(fio, ir->userreal3);
1528 gmx_fio_do_real(fio, ir->userreal4);
1529
1530 /* AdResS is removed, but we need to be able to read old files,
1531 and let mdrun refuse to run them */
1532 if (file_version >= 77 && file_version < tpxv_RemoveAdress)
1533 {
1534 gmx_fio_do_gmx_bool(fio, ir->bAdress);
1535 if (ir->bAdress)
1536 {
1537 int idum, numThermoForceGroups, numEnergyGroups;
1538 real rdum;
1539 rvec rvecdum;
1540 gmx_fio_do_int(fio, idum);
1541 gmx_fio_do_real(fio, rdum);
1542 gmx_fio_do_real(fio, rdum);
1543 gmx_fio_do_real(fio, rdum);
1544 gmx_fio_do_int(fio, idum);
1545 gmx_fio_do_int(fio, idum);
1546 gmx_fio_do_rvec(fio, rvecdum);
1547 gmx_fio_do_int(fio, numThermoForceGroups);
1548 gmx_fio_do_real(fio, rdum);
1549 gmx_fio_do_int(fio, numEnergyGroups);
1550 gmx_fio_do_int(fio, idum);
1551
1552 if (numThermoForceGroups > 0)
1553 {
1554 std::vector<int> idumn(numThermoForceGroups);
1555 gmx_fio_ndo_int(fio, idumn.data(), idumn.size());
1556 }
1557 if (numEnergyGroups > 0)
1558 {
1559 std::vector<int> idumn(numEnergyGroups);
1560 gmx_fio_ndo_int(fio, idumn.data(), idumn.size());
1561 }
1562 }
1563 }
1564 else
1565 {
1566 ir->bAdress = FALSE;
1567 }
1568
1569 /* pull stuff */
1570 if (file_version >= 48)
1571 {
1572 int ePullOld = 0;
1573
1574 if (file_version >= tpxv_PullCoordTypeGeom)
1575 {
1576 gmx_fio_do_gmx_bool(fio, ir->bPull);
1577 }
1578 else
1579 {
1580 gmx_fio_do_int(fio, ePullOld);
1581 ir->bPull = (ePullOld > 0);
1582 /* We removed the first ePull=ePullNo for the enum */
1583 ePullOld -= 1;
1584 }
1585 if (ir->bPull)
1586 {
1587 if (bRead)
1588 {
1589 snew(ir->pull, 1);
1590 }
1591 do_pull(fio, ir->pull, bRead, file_version, ePullOld);
1592 }
1593 }
1594 else
1595 {
1596 ir->bPull = FALSE;
1597 }
1598
1599 if (file_version >= tpxv_AcceleratedWeightHistogram)
1600 {
1601 gmx_fio_do_gmx_bool(fio, ir->bDoAwh);
1602
1603 if (ir->bDoAwh)
1604 {
1605 if (bRead)
1606 {
1607 snew(ir->awhParams, 1);
1608 }
1609 do_awh(fio, ir->awhParams, bRead, file_version);
1610 }
1611 }
1612 else
1613 {
1614 ir->bDoAwh = FALSE;
1615 }
1616
1617 /* Enforced rotation */
1618 if (file_version >= 74)
1619 {
1620 gmx_fio_do_int(fio, ir->bRot);
1621 if (ir->bRot == TRUE)
1622 {
1623 if (bRead)
1624 {
1625 snew(ir->rot, 1);
1626 }
1627 do_rot(fio, ir->rot, bRead);
1628 }
1629 }
1630 else
1631 {
1632 ir->bRot = FALSE;
1633 }
1634
1635 /* Interactive molecular dynamics */
1636 if (file_version >= tpxv_InteractiveMolecularDynamics)
1637 {
1638 gmx_fio_do_int(fio, ir->bIMD);
1639 if (TRUE == ir->bIMD)
1640 {
1641 if (bRead)
1642 {
1643 snew(ir->imd, 1);
1644 }
1645 do_imd(fio, ir->imd, bRead);
1646 }
1647 }
1648 else
1649 {
1650 /* We don't support IMD sessions for old .tpr files */
1651 ir->bIMD = FALSE;
1652 }
1653
1654 /* grpopts stuff */
1655 gmx_fio_do_int(fio, ir->opts.ngtc);
1656 if (file_version >= 69)
1657 {
1658 gmx_fio_do_int(fio, ir->opts.nhchainlength);
1659 }
1660 else
1661 {
1662 ir->opts.nhchainlength = 1;
1663 }
1664 gmx_fio_do_int(fio, ir->opts.ngacc);
1665 gmx_fio_do_int(fio, ir->opts.ngfrz);
1666 gmx_fio_do_int(fio, ir->opts.ngener);
1667
1668 if (bRead)
1669 {
1670 snew(ir->opts.nrdf, ir->opts.ngtc);
1671 snew(ir->opts.ref_t, ir->opts.ngtc);
1672 snew(ir->opts.annealing, ir->opts.ngtc);
1673 snew(ir->opts.anneal_npoints, ir->opts.ngtc);
1674 snew(ir->opts.anneal_time, ir->opts.ngtc);
1675 snew(ir->opts.anneal_temp, ir->opts.ngtc);
1676 snew(ir->opts.tau_t, ir->opts.ngtc);
1677 snew(ir->opts.nFreeze, ir->opts.ngfrz);
1678 snew(ir->opts.acc, ir->opts.ngacc);
1679 snew(ir->opts.egp_flags, ir->opts.ngener*ir->opts.ngener);
1680 }
1681 if (ir->opts.ngtc > 0)
1682 {
1683 gmx_fio_ndo_real(fio, ir->opts.nrdf, ir->opts.ngtc);
1684 gmx_fio_ndo_real(fio, ir->opts.ref_t, ir->opts.ngtc);
1685 gmx_fio_ndo_real(fio, ir->opts.tau_t, ir->opts.ngtc);
1686 if (file_version < 33 && ir->eI == eiBD)
1687 {
1688 for (i = 0; i < ir->opts.ngtc; i++)
1689 {
1690 ir->opts.tau_t[i] = bd_temp;
1691 }
1692 }
1693 }
1694 if (ir->opts.ngfrz > 0)
1695 {
1696 gmx_fio_ndo_ivec(fio, ir->opts.nFreeze, ir->opts.ngfrz);
1697 }
1698 if (ir->opts.ngacc > 0)
1699 {
1700 gmx_fio_ndo_rvec(fio, ir->opts.acc, ir->opts.ngacc);
1701 }
1702 gmx_fio_ndo_int(fio, ir->opts.egp_flags,
1703 ir->opts.ngener*ir->opts.ngener);
1704
1705 /* First read the lists with annealing and npoints for each group */
1706 gmx_fio_ndo_int(fio, ir->opts.annealing, ir->opts.ngtc);
1707 gmx_fio_ndo_int(fio, ir->opts.anneal_npoints, ir->opts.ngtc);
1708 for (j = 0; j < (ir->opts.ngtc); j++)
1709 {
1710 k = ir->opts.anneal_npoints[j];
1711 if (bRead)
1712 {
1713 snew(ir->opts.anneal_time[j], k);
1714 snew(ir->opts.anneal_temp[j], k);
1715 }
1716 gmx_fio_ndo_real(fio, ir->opts.anneal_time[j], k);
1717 gmx_fio_ndo_real(fio, ir->opts.anneal_temp[j], k);
1718 }
1719 /* Walls */
1720 if (file_version >= 45)
1721 {
1722 gmx_fio_do_int(fio, ir->nwall);
1723 gmx_fio_do_int(fio, ir->wall_type);
1724 if (file_version >= 50)
1725 {
1726 gmx_fio_do_real(fio, ir->wall_r_linpot);
1727 }
1728 else
1729 {
1730 ir->wall_r_linpot = -1;
1731 }
1732 gmx_fio_do_int(fio, ir->wall_atomtype[0]);
1733 gmx_fio_do_int(fio, ir->wall_atomtype[1]);
1734 gmx_fio_do_real(fio, ir->wall_density[0]);
1735 gmx_fio_do_real(fio, ir->wall_density[1]);
1736 gmx_fio_do_real(fio, ir->wall_ewald_zfac);
1737 }
1738 else
1739 {
1740 ir->nwall = 0;
1741 ir->wall_type = 0;
1742 ir->wall_atomtype[0] = -1;
1743 ir->wall_atomtype[1] = -1;
1744 ir->wall_density[0] = 0;
1745 ir->wall_density[1] = 0;
1746 ir->wall_ewald_zfac = 3;
1747 }
1748 /* Cosine stuff for electric fields */
1749 if (file_version < tpxv_GenericParamsForElectricField)
1750 {
1751 do_legacy_efield(fio, &paramsObj);
1752 }
1753
1754 /* Swap ions */
1755 if (file_version >= tpxv_ComputationalElectrophysiology)
1756 {
1757 gmx_fio_do_int(fio, ir->eSwapCoords);
1758 if (ir->eSwapCoords != eswapNO)
1759 {
1760 if (bRead)
1761 {
1762 snew(ir->swap, 1);
1763 }
1764 do_swapcoords_tpx(fio, ir->swap, bRead, file_version);
1765 }
1766 }
1767
1768 /* QMMM stuff */
1769 if (file_version >= 39)
1770 {
1771 gmx_fio_do_gmx_bool(fio, ir->bQMMM);
1772 gmx_fio_do_int(fio, ir->QMMMscheme);
1773 gmx_fio_do_real(fio, ir->scalefactor);
1774 gmx_fio_do_int(fio, ir->opts.ngQM);
1775 if (bRead)
1776 {
1777 snew(ir->opts.QMmethod, ir->opts.ngQM);
1778 snew(ir->opts.QMbasis, ir->opts.ngQM);
1779 snew(ir->opts.QMcharge, ir->opts.ngQM);
1780 snew(ir->opts.QMmult, ir->opts.ngQM);
1781 snew(ir->opts.bSH, ir->opts.ngQM);
1782 snew(ir->opts.CASorbitals, ir->opts.ngQM);
1783 snew(ir->opts.CASelectrons, ir->opts.ngQM);
1784 snew(ir->opts.SAon, ir->opts.ngQM);
1785 snew(ir->opts.SAoff, ir->opts.ngQM);
1786 snew(ir->opts.SAsteps, ir->opts.ngQM);
1787 }
1788 if (ir->opts.ngQM > 0)
1789 {
1790 gmx_fio_ndo_int(fio, ir->opts.QMmethod, ir->opts.ngQM);
1791 gmx_fio_ndo_int(fio, ir->opts.QMbasis, ir->opts.ngQM);
1792 gmx_fio_ndo_int(fio, ir->opts.QMcharge, ir->opts.ngQM);
1793 gmx_fio_ndo_int(fio, ir->opts.QMmult, ir->opts.ngQM);
1794 gmx_fio_ndo_gmx_bool(fio, ir->opts.bSH, ir->opts.ngQM);
1795 gmx_fio_ndo_int(fio, ir->opts.CASorbitals, ir->opts.ngQM);
1796 gmx_fio_ndo_int(fio, ir->opts.CASelectrons, ir->opts.ngQM);
1797 gmx_fio_ndo_real(fio, ir->opts.SAon, ir->opts.ngQM);
1798 gmx_fio_ndo_real(fio, ir->opts.SAoff, ir->opts.ngQM);
1799 gmx_fio_ndo_int(fio, ir->opts.SAsteps, ir->opts.ngQM);
1800 /* We leave in dummy i/o for removed parameters to avoid
1801 * changing the tpr format for every QMMM change.
1802 */
1803 std::vector<int> dummy(ir->opts.ngQM, 0);
1804 gmx_fio_ndo_gmx_bool(fio, dummy.data(), ir->opts.ngQM);
1805 gmx_fio_ndo_gmx_bool(fio, dummy.data(), ir->opts.ngQM);
1806 }
1807 /* end of QMMM stuff */
1808 }
1809
1810 if (file_version >= tpxv_GenericParamsForElectricField)
1811 {
1812 gmx::FileIOXdrSerializer serializer(fio);
1813 if (bRead)
1814 {
1815 paramsObj.mergeObject(
1816 gmx::deserializeKeyValueTree(&serializer));
1817 }
1818 else
1819 {
1820 GMX_RELEASE_ASSERT(ir->params != nullptr,
1821 "Parameters should be present when writing inputrec");
1822 gmx::serializeKeyValueTree(*ir->params, &serializer);
1823 }
1824 }
1825 if (bRead)
1826 {
1827 ir->params = new gmx::KeyValueTreeObject(paramsBuilder.build());
1828 }
1829 }
1830
1831
1832 static void do_harm(t_fileio *fio, t_iparams *iparams)
1833 {
1834 gmx_fio_do_real(fio, iparams->harmonic.rA);
1835 gmx_fio_do_real(fio, iparams->harmonic.krA);
1836 gmx_fio_do_real(fio, iparams->harmonic.rB);
1837 gmx_fio_do_real(fio, iparams->harmonic.krB);
1838 }
1839
1840 static void do_iparams(t_fileio *fio, t_functype ftype, t_iparams *iparams,
1841 gmx_bool bRead, int file_version)
1842 {
1843 int idum;
1844 real rdum;
1845
1846 switch (ftype)
1847 {
1848 case F_ANGLES:
1849 case F_G96ANGLES:
1850 case F_BONDS:
1851 case F_G96BONDS:
1852 case F_HARMONIC:
1853 case F_IDIHS:
1854 do_harm(fio, iparams);
1855 if ((ftype == F_ANGRES || ftype == F_ANGRESZ) && bRead)
1856 {
1857 /* Correct incorrect storage of parameters */
1858 iparams->pdihs.phiB = iparams->pdihs.phiA;
1859 iparams->pdihs.cpB = iparams->pdihs.cpA;
1860 }
1861 break;
1862 case F_RESTRANGLES:
1863 gmx_fio_do_real(fio, iparams->harmonic.rA);
1864 gmx_fio_do_real(fio, iparams->harmonic.krA);
1865 break;
1866 case F_LINEAR_ANGLES:
1867 gmx_fio_do_real(fio, iparams->linangle.klinA);
1868 gmx_fio_do_real(fio, iparams->linangle.aA);
1869 gmx_fio_do_real(fio, iparams->linangle.klinB);
1870 gmx_fio_do_real(fio, iparams->linangle.aB);
1871 break;
1872 case F_FENEBONDS:
1873 gmx_fio_do_real(fio, iparams->fene.bm);
1874 gmx_fio_do_real(fio, iparams->fene.kb);
1875 break;
1876
1877 case F_RESTRBONDS:
1878 gmx_fio_do_real(fio, iparams->restraint.lowA);
1879 gmx_fio_do_real(fio, iparams->restraint.up1A);
1880 gmx_fio_do_real(fio, iparams->restraint.up2A);
1881 gmx_fio_do_real(fio, iparams->restraint.kA);
1882 gmx_fio_do_real(fio, iparams->restraint.lowB);
1883 gmx_fio_do_real(fio, iparams->restraint.up1B);
1884 gmx_fio_do_real(fio, iparams->restraint.up2B);
1885 gmx_fio_do_real(fio, iparams->restraint.kB);
1886 break;
1887 case F_TABBONDS:
1888 case F_TABBONDSNC:
1889 case F_TABANGLES:
1890 case F_TABDIHS:
1891 gmx_fio_do_real(fio, iparams->tab.kA);
1892 gmx_fio_do_int(fio, iparams->tab.table);
1893 gmx_fio_do_real(fio, iparams->tab.kB);
1894 break;
1895 case F_CROSS_BOND_BONDS:
1896 gmx_fio_do_real(fio, iparams->cross_bb.r1e);
1897 gmx_fio_do_real(fio, iparams->cross_bb.r2e);
1898 gmx_fio_do_real(fio, iparams->cross_bb.krr);
1899 break;
1900 case F_CROSS_BOND_ANGLES:
1901 gmx_fio_do_real(fio, iparams->cross_ba.r1e);
1902 gmx_fio_do_real(fio, iparams->cross_ba.r2e);
1903 gmx_fio_do_real(fio, iparams->cross_ba.r3e);
1904 gmx_fio_do_real(fio, iparams->cross_ba.krt);
1905 break;
1906 case F_UREY_BRADLEY:
1907 gmx_fio_do_real(fio, iparams->u_b.thetaA);
1908 gmx_fio_do_real(fio, iparams->u_b.kthetaA);
1909 gmx_fio_do_real(fio, iparams->u_b.r13A);
1910 gmx_fio_do_real(fio, iparams->u_b.kUBA);
1911 if (file_version >= 79)
1912 {
1913 gmx_fio_do_real(fio, iparams->u_b.thetaB);
1914 gmx_fio_do_real(fio, iparams->u_b.kthetaB);
1915 gmx_fio_do_real(fio, iparams->u_b.r13B);
1916 gmx_fio_do_real(fio, iparams->u_b.kUBB);
1917 }
1918 else
1919 {
1920 iparams->u_b.thetaB = iparams->u_b.thetaA;
1921 iparams->u_b.kthetaB = iparams->u_b.kthetaA;
1922 iparams->u_b.r13B = iparams->u_b.r13A;
1923 iparams->u_b.kUBB = iparams->u_b.kUBA;
1924 }
1925 break;
1926 case F_QUARTIC_ANGLES:
1927 gmx_fio_do_real(fio, iparams->qangle.theta);
1928 gmx_fio_ndo_real(fio, iparams->qangle.c, 5);
1929 break;
1930 case F_BHAM:
1931 gmx_fio_do_real(fio, iparams->bham.a);
1932 gmx_fio_do_real(fio, iparams->bham.b);
1933 gmx_fio_do_real(fio, iparams->bham.c);
1934 break;
1935 case F_MORSE:
1936 gmx_fio_do_real(fio, iparams->morse.b0A);
1937 gmx_fio_do_real(fio, iparams->morse.cbA);
1938 gmx_fio_do_real(fio, iparams->morse.betaA);
1939 if (file_version >= 79)
1940 {
1941 gmx_fio_do_real(fio, iparams->morse.b0B);
1942 gmx_fio_do_real(fio, iparams->morse.cbB);
1943 gmx_fio_do_real(fio, iparams->morse.betaB);
1944 }
1945 else
1946 {
1947 iparams->morse.b0B = iparams->morse.b0A;
1948 iparams->morse.cbB = iparams->morse.cbA;
1949 iparams->morse.betaB = iparams->morse.betaA;
1950 }
1951 break;
1952 case F_CUBICBONDS:
1953 gmx_fio_do_real(fio, iparams->cubic.b0);
1954 gmx_fio_do_real(fio, iparams->cubic.kb);
1955 gmx_fio_do_real(fio, iparams->cubic.kcub);
1956 break;
1957 case F_CONNBONDS:
1958 break;
1959 case F_POLARIZATION:
1960 gmx_fio_do_real(fio, iparams->polarize.alpha);
1961 break;
1962 case F_ANHARM_POL:
1963 gmx_fio_do_real(fio, iparams->anharm_polarize.alpha);
1964 gmx_fio_do_real(fio, iparams->anharm_polarize.drcut);
1965 gmx_fio_do_real(fio, iparams->anharm_polarize.khyp);
1966 break;
1967 case F_WATER_POL:
1968 gmx_fio_do_real(fio, iparams->wpol.al_x);
1969 gmx_fio_do_real(fio, iparams->wpol.al_y);
1970 gmx_fio_do_real(fio, iparams->wpol.al_z);
1971 gmx_fio_do_real(fio, iparams->wpol.rOH);
1972 gmx_fio_do_real(fio, iparams->wpol.rHH);
1973 gmx_fio_do_real(fio, iparams->wpol.rOD);
1974 break;
1975 case F_THOLE_POL:
1976 gmx_fio_do_real(fio, iparams->thole.a);
1977 gmx_fio_do_real(fio, iparams->thole.alpha1);
1978 gmx_fio_do_real(fio, iparams->thole.alpha2);
1979 gmx_fio_do_real(fio, iparams->thole.rfac);
1980 break;
1981 case F_LJ:
1982 gmx_fio_do_real(fio, iparams->lj.c6);
1983 gmx_fio_do_real(fio, iparams->lj.c12);
1984 break;
1985 case F_LJ14:
1986 gmx_fio_do_real(fio, iparams->lj14.c6A);
1987 gmx_fio_do_real(fio, iparams->lj14.c12A);
1988 gmx_fio_do_real(fio, iparams->lj14.c6B);
1989 gmx_fio_do_real(fio, iparams->lj14.c12B);
1990 break;
1991 case F_LJC14_Q:
1992 gmx_fio_do_real(fio, iparams->ljc14.fqq);
1993 gmx_fio_do_real(fio, iparams->ljc14.qi);
1994 gmx_fio_do_real(fio, iparams->ljc14.qj);
1995 gmx_fio_do_real(fio, iparams->ljc14.c6);
1996 gmx_fio_do_real(fio, iparams->ljc14.c12);
1997 break;
1998 case F_LJC_PAIRS_NB:
1999 gmx_fio_do_real(fio, iparams->ljcnb.qi);
2000 gmx_fio_do_real(fio, iparams->ljcnb.qj);
2001 gmx_fio_do_real(fio, iparams->ljcnb.c6);
2002 gmx_fio_do_real(fio, iparams->ljcnb.c12);
2003 break;
2004 case F_PDIHS:
2005 case F_PIDIHS:
2006 case F_ANGRES:
2007 case F_ANGRESZ:
2008 gmx_fio_do_real(fio, iparams->pdihs.phiA);
2009 gmx_fio_do_real(fio, iparams->pdihs.cpA);
2010 if ((ftype == F_ANGRES || ftype == F_ANGRESZ) && file_version < 42)
2011 {
2012 /* Read the incorrectly stored multiplicity */
2013 gmx_fio_do_real(fio, iparams->harmonic.rB);
2014 gmx_fio_do_real(fio, iparams->harmonic.krB);
2015 iparams->pdihs.phiB = iparams->pdihs.phiA;
2016 iparams->pdihs.cpB = iparams->pdihs.cpA;
2017 }
2018 else
2019 {
2020 gmx_fio_do_real(fio, iparams->pdihs.phiB);
2021 gmx_fio_do_real(fio, iparams->pdihs.cpB);
2022 gmx_fio_do_int(fio, iparams->pdihs.mult);
2023 }
2024 break;
2025 case F_RESTRDIHS:
2026 gmx_fio_do_real(fio, iparams->pdihs.phiA);
2027 gmx_fio_do_real(fio, iparams->pdihs.cpA);
2028 break;
2029 case F_DISRES:
2030 gmx_fio_do_int(fio, iparams->disres.label);
2031 gmx_fio_do_int(fio, iparams->disres.type);
2032 gmx_fio_do_real(fio, iparams->disres.low);
2033 gmx_fio_do_real(fio, iparams->disres.up1);
2034 gmx_fio_do_real(fio, iparams->disres.up2);
2035 gmx_fio_do_real(fio, iparams->disres.kfac);
2036 break;
2037 case F_ORIRES:
2038 gmx_fio_do_int(fio, iparams->orires.ex);
2039 gmx_fio_do_int(fio, iparams->orires.label);
2040 gmx_fio_do_int(fio, iparams->orires.power);
2041 gmx_fio_do_real(fio, iparams->orires.c);
2042 gmx_fio_do_real(fio, iparams->orires.obs);
2043 gmx_fio_do_real(fio, iparams->orires.kfac);
2044 break;
2045 case F_DIHRES:
2046 if (file_version < 82)
2047 {
2048 gmx_fio_do_int(fio, idum);
2049 gmx_fio_do_int(fio, idum);
2050 }
2051 gmx_fio_do_real(fio, iparams->dihres.phiA);
2052 gmx_fio_do_real(fio, iparams->dihres.dphiA);
2053 gmx_fio_do_real(fio, iparams->dihres.kfacA);
2054 if (file_version >= 82)
2055 {
2056 gmx_fio_do_real(fio, iparams->dihres.phiB);
2057 gmx_fio_do_real(fio, iparams->dihres.dphiB);
2058 gmx_fio_do_real(fio, iparams->dihres.kfacB);
2059 }
2060 else
2061 {
2062 iparams->dihres.phiB = iparams->dihres.phiA;
2063 iparams->dihres.dphiB = iparams->dihres.dphiA;
2064 iparams->dihres.kfacB = iparams->dihres.kfacA;
2065 }
2066 break;
2067 case F_POSRES:
2068 gmx_fio_do_rvec(fio, iparams->posres.pos0A);
2069 gmx_fio_do_rvec(fio, iparams->posres.fcA);
2070 gmx_fio_do_rvec(fio, iparams->posres.pos0B);
2071 gmx_fio_do_rvec(fio, iparams->posres.fcB);
2072 break;
2073 case F_FBPOSRES:
2074 gmx_fio_do_int(fio, iparams->fbposres.geom);
2075 gmx_fio_do_rvec(fio, iparams->fbposres.pos0);
2076 gmx_fio_do_real(fio, iparams->fbposres.r);
2077 gmx_fio_do_real(fio, iparams->fbposres.k);
2078 break;
2079 case F_CBTDIHS:
2080 gmx_fio_ndo_real(fio, iparams->cbtdihs.cbtcA, NR_CBTDIHS);
2081 break;
2082 case F_RBDIHS:
2083 gmx_fio_ndo_real(fio, iparams->rbdihs.rbcA, NR_RBDIHS);
2084 gmx_fio_ndo_real(fio, iparams->rbdihs.rbcB, NR_RBDIHS);
2085 break;
2086 case F_FOURDIHS:
2087 /* Fourier dihedrals are internally represented
2088 * as Ryckaert-Bellemans since those are faster to compute.
2089 */
2090 gmx_fio_ndo_real(fio, iparams->rbdihs.rbcA, NR_RBDIHS);
2091 gmx_fio_ndo_real(fio, iparams->rbdihs.rbcB, NR_RBDIHS);
2092 break;
2093 case F_CONSTR:
2094 case F_CONSTRNC:
2095 gmx_fio_do_real(fio, iparams->constr.dA);
2096 gmx_fio_do_real(fio, iparams->constr.dB);
2097 break;
2098 case F_SETTLE:
2099 gmx_fio_do_real(fio, iparams->settle.doh);
2100 gmx_fio_do_real(fio, iparams->settle.dhh);
2101 break;
2102 case F_VSITE2:
2103 gmx_fio_do_real(fio, iparams->vsite.a);
2104 break;
2105 case F_VSITE3:
2106 case F_VSITE3FD:
2107 case F_VSITE3FAD:
2108 gmx_fio_do_real(fio, iparams->vsite.a);
2109 gmx_fio_do_real(fio, iparams->vsite.b);
2110 break;
2111 case F_VSITE3OUT:
2112 case F_VSITE4FD:
2113 case F_VSITE4FDN:
2114 gmx_fio_do_real(fio, iparams->vsite.a);
2115 gmx_fio_do_real(fio, iparams->vsite.b);
2116 gmx_fio_do_real(fio, iparams->vsite.c);
2117 break;
2118 case F_VSITEN:
2119 gmx_fio_do_int(fio, iparams->vsiten.n);
2120 gmx_fio_do_real(fio, iparams->vsiten.a);
2121 break;
2122 case F_GB12:
2123 case F_GB13:
2124 case F_GB14:
2125 /* We got rid of some parameters in version 68 */
2126 if (bRead && file_version < 68)
2127 {
2128 gmx_fio_do_real(fio, rdum);
2129 gmx_fio_do_real(fio, rdum);
2130 gmx_fio_do_real(fio, rdum);
2131 gmx_fio_do_real(fio, rdum);
2132 }
2133 gmx_fio_do_real(fio, iparams->gb.sar);
2134 gmx_fio_do_real(fio, iparams->gb.st);
2135 gmx_fio_do_real(fio, iparams->gb.pi);
2136 gmx_fio_do_real(fio, iparams->gb.gbr);
2137 gmx_fio_do_real(fio, iparams->gb.bmlt);
2138 break;
2139 case F_CMAP:
2140 gmx_fio_do_int(fio, iparams->cmap.cmapA);
2141 gmx_fio_do_int(fio, iparams->cmap.cmapB);
2142 break;
2143 default:
2144 gmx_fatal(FARGS, "unknown function type %d (%s) in %s line %d",
2145 ftype, interaction_function[ftype].name, __FILE__, __LINE__);
2146 }
2147 }
2148
2149 static void do_ilist(t_fileio *fio, t_ilist *ilist, gmx_bool bRead, int file_version)
2150 {
2151 int i, idum;
2152
2153 if (file_version < 44)
2154 {
2155 for (i = 0; i < MAXNODES; i++)
2156 {
2157 gmx_fio_do_int(fio, idum);
2158 }
2159 }
2160 gmx_fio_do_int(fio, ilist->nr);
2161 if (bRead)
2162 {
2163 snew(ilist->iatoms, ilist->nr);
2164 }
2165 gmx_fio_ndo_int(fio, ilist->iatoms, ilist->nr);
2166 }
2167
2168 static void do_ffparams(t_fileio *fio, gmx_ffparams_t *ffparams,
2169 gmx_bool bRead, int file_version)
2170 {
2171 int idum, i;
2172 unsigned int k;
2173
2174 gmx_fio_do_int(fio, ffparams->atnr);
2175 if (file_version < 57)
2176 {
2177 gmx_fio_do_int(fio, idum);
2178 }
2179 gmx_fio_do_int(fio, ffparams->ntypes);
2180 if (bRead)
2181 {
2182 snew(ffparams->functype, ffparams->ntypes);
2183 snew(ffparams->iparams, ffparams->ntypes);
2184 }
2185 /* Read/write all the function types */
2186 gmx_fio_ndo_int(fio, ffparams->functype, ffparams->ntypes);
2187
2188 if (file_version >= 66)
2189 {
2190 gmx_fio_do_double(fio, ffparams->reppow);
2191 }
2192 else
2193 {
2194 ffparams->reppow = 12.0;
2195 }
2196
2197 if (file_version >= 57)
2198 {
2199 gmx_fio_do_real(fio, ffparams->fudgeQQ);
2200 }
2201
2202 /* Check whether all these function types are supported by the code.
2203 * In practice the code is backwards compatible, which means that the
2204 * numbering may have to be altered from old numbering to new numbering
2205 */
2206 for (i = 0; (i < ffparams->ntypes); i++)
2207 {
2208 if (bRead)
2209 {
2210 /* Loop over file versions */
2211 for (k = 0; (k < NFTUPD); k++)
2212 {
2213 /* Compare the read file_version to the update table */
2214 if ((file_version < ftupd[k].fvnr) &&
2215 (ffparams->functype[i] >= ftupd[k].ftype))
2216 {
2217 ffparams->functype[i] += 1;
2218 }
2219 }
2220 }
2221
2222 do_iparams(fio, ffparams->functype[i], &ffparams->iparams[i], bRead,
2223 file_version);
2224 }
2225 }
2226
2227 static void add_settle_atoms(t_ilist *ilist)
2228 {
2229 int i;
2230
2231 /* Settle used to only store the first atom: add the other two */
2232 srenew(ilist->iatoms, 2*ilist->nr);
2233 for (i = ilist->nr/2-1; i >= 0; i--)
2234 {
2235 ilist->iatoms[4*i+0] = ilist->iatoms[2*i+0];
2236 ilist->iatoms[4*i+1] = ilist->iatoms[2*i+1];
2237 ilist->iatoms[4*i+2] = ilist->iatoms[2*i+1] + 1;
2238 ilist->iatoms[4*i+3] = ilist->iatoms[2*i+1] + 2;
2239 }
2240 ilist->nr = 2*ilist->nr;
2241 }
2242
2243 static void do_ilists(t_fileio *fio, t_ilist *ilist, gmx_bool bRead,
2244 int file_version)
2245 {
2246 int j;
2247 gmx_bool bClear;
2248 unsigned int k;
2249
2250 for (j = 0; (j < F_NRE); j++)
2251 {
2252 bClear = FALSE;
2253 if (bRead)
2254 {
2255 for (k = 0; k < NFTUPD; k++)
2256 {
2257 if ((file_version < ftupd[k].fvnr) && (j == ftupd[k].ftype))
2258 {
2259 bClear = TRUE;
2260 }
2261 }
2262 }
2263 if (bClear)
2264 {
2265 ilist[j].nr = 0;
2266 ilist[j].iatoms = nullptr;
2267 }
2268 else
2269 {
2270 do_ilist(fio, &ilist[j], bRead, file_version);
2271 if (file_version < 78 && j == F_SETTLE && ilist[j].nr > 0)
2272 {
2273 add_settle_atoms(&ilist[j]);
2274 }
2275 }
2276 }
2277 }
2278
2279 static void do_idef(t_fileio *fio, gmx_ffparams_t *ffparams, gmx_moltype_t *molt,
2280 gmx_bool bRead, int file_version)
2281 {
2282 do_ffparams(fio, ffparams, bRead, file_version);
2283
2284 if (file_version >= 54)
2285 {
2286 gmx_fio_do_real(fio, ffparams->fudgeQQ);
2287 }
2288
2289 do_ilists(fio, molt->ilist, bRead, file_version);
2290 }
2291
2292 static void do_block(t_fileio *fio, t_block *block, gmx_bool bRead, int file_version)
2293 {
2294 int i, idum, dum_nra, *dum_a;
2295
2296 if (file_version < 44)
2297 {
2298 for (i = 0; i < MAXNODES; i++)
2299 {
2300 gmx_fio_do_int(fio, idum);
2301 }
2302 }
2303 gmx_fio_do_int(fio, block->nr);
2304 if (file_version < 51)
2305 {
2306 gmx_fio_do_int(fio, dum_nra);
2307 }
2308 if (bRead)
2309 {
2310 if ((block->nalloc_index > 0) && (nullptr != block->index))
2311 {
2312 sfree(block->index);
2313 }
2314 block->nalloc_index = block->nr+1;
2315 snew(block->index, block->nalloc_index);
2316 }
2317 gmx_fio_ndo_int(fio, block->index, block->nr+1);
2318
2319 if (file_version < 51 && dum_nra > 0)
2320 {
2321 snew(dum_a, dum_nra);
2322 gmx_fio_ndo_int(fio, dum_a, dum_nra);
2323 sfree(dum_a);
2324 }
2325 }
2326
2327 static void do_blocka(t_fileio *fio, t_blocka *block, gmx_bool bRead,
2328 int file_version)
2329 {
2330 int i, idum;
2331
2332 if (file_version < 44)
2333 {
2334 for (i = 0; i < MAXNODES; i++)
2335 {
2336 gmx_fio_do_int(fio, idum);
2337 }
2338 }
2339 gmx_fio_do_int(fio, block->nr);
2340 gmx_fio_do_int(fio, block->nra);
2341 if (bRead)
2342 {
2343 block->nalloc_index = block->nr+1;
2344 snew(block->index, block->nalloc_index);
2345 block->nalloc_a = block->nra;
2346 snew(block->a, block->nalloc_a);
2347 }
2348 gmx_fio_ndo_int(fio, block->index, block->nr+1);
2349 gmx_fio_ndo_int(fio, block->a, block->nra);
2350 }
2351
2352 /* This is a primitive routine to make it possible to translate atomic numbers
2353 * to element names when reading TPR files, without making the Gromacs library
2354 * directory a dependency on mdrun (which is the case if we need elements.dat).
2355 */
2356 static const char *
2357 atomicnumber_to_element(int atomicnumber)
2358 {
2359 const char * p;
2360
2361 /* This does not have to be complete, so we only include elements likely
2362 * to occur in PDB files.
2363 */
2364 switch (atomicnumber)
2365 {
2366 case 1: p = "H"; break;
2367 case 5: p = "B"; break;
2368 case 6: p = "C"; break;
2369 case 7: p = "N"; break;
2370 case 8: p = "O"; break;
2371 case 9: p = "F"; break;
2372 case 11: p = "Na"; break;
2373 case 12: p = "Mg"; break;
2374 case 15: p = "P"; break;
2375 case 16: p = "S"; break;
2376 case 17: p = "Cl"; break;
2377 case 18: p = "Ar"; break;
2378 case 19: p = "K"; break;
2379 case 20: p = "Ca"; break;
2380 case 25: p = "Mn"; break;
2381 case 26: p = "Fe"; break;
2382 case 28: p = "Ni"; break;
2383 case 29: p = "Cu"; break;
2384 case 30: p = "Zn"; break;
2385 case 35: p = "Br"; break;
2386 case 47: p = "Ag"; break;
2387 default: p = ""; break;
2388 }
2389 return p;
2390 }
2391
2392
2393 static void do_atom(t_fileio *fio, t_atom *atom, int ngrp, gmx_bool bRead,
2394 int file_version, gmx_groups_t *groups, int atnr)
2395 {
2396 int i, myngrp;
2397
2398 gmx_fio_do_real(fio, atom->m);
2399 gmx_fio_do_real(fio, atom->q);
2400 gmx_fio_do_real(fio, atom->mB);
2401 gmx_fio_do_real(fio, atom->qB);
2402 gmx_fio_do_ushort(fio, atom->type);
2403 gmx_fio_do_ushort(fio, atom->typeB);
2404 gmx_fio_do_int(fio, atom->ptype);
2405 gmx_fio_do_int(fio, atom->resind);
2406 if (file_version >= 52)
2407 {
2408 gmx_fio_do_int(fio, atom->atomnumber);
2409 if (bRead)
2410 {
2411 /* Set element string from atomic number if present.
2412 * This routine returns an empty string if the name is not found.
2413 */
2414 std::strncpy(atom->elem, atomicnumber_to_element(atom->atomnumber), 4);
2415 /* avoid warnings about potentially unterminated string */
2416 atom->elem[3] = '\0';
2417 }
2418 }
2419 else if (bRead)
2420 {
2421 atom->atomnumber = 0;
2422 }
2423 if (file_version < 39)
2424 {
2425 myngrp = 9;
2426 }
2427 else
2428 {
2429 myngrp = ngrp;
2430 }
2431
2432 if (file_version < 57)
2433 {
2434 unsigned char uchar[egcNR];
2435 gmx_fio_ndo_uchar(fio, uchar, myngrp);
2436 for (i = myngrp; (i < ngrp); i++)
2437 {
2438 uchar[i] = 0;
2439 }
2440 /* Copy the old data format to the groups struct */
2441 for (i = 0; i < ngrp; i++)
2442 {
2443 groups->grpnr[i][atnr] = uchar[i];
2444 }
2445 }
2446 }
2447
2448 static void do_grps(t_fileio *fio, int ngrp, t_grps grps[], gmx_bool bRead,
2449 int file_version)
2450 {
2451 int j, myngrp;
2452
2453 if (file_version < 39)
2454 {
2455 myngrp = 9;
2456 }
2457 else
2458 {
2459 myngrp = ngrp;
2460 }
2461
2462 for (j = 0; (j < ngrp); j++)
2463 {
2464 if (j < myngrp)
2465 {
2466 gmx_fio_do_int(fio, grps[j].nr);
2467 if (bRead)
2468 {
2469 snew(grps[j].nm_ind, grps[j].nr);
2470 }
2471 gmx_fio_ndo_int(fio, grps[j].nm_ind, grps[j].nr);
2472 }
2473 else
2474 {
2475 grps[j].nr = 1;
2476 snew(grps[j].nm_ind, grps[j].nr);
2477 }
2478 }
2479 }
2480
2481 static void do_symstr(t_fileio *fio, char ***nm, gmx_bool bRead, t_symtab *symtab)
2482 {
2483 int ls;
2484
2485 if (bRead)
2486 {
2487 gmx_fio_do_int(fio, ls);
2488 *nm = get_symtab_handle(symtab, ls);
2489 }
2490 else
2491 {
2492 ls = lookup_symtab(symtab, *nm);
2493 gmx_fio_do_int(fio, ls);
2494 }
2495 }
2496
2497 static void do_strstr(t_fileio *fio, int nstr, char ***nm, gmx_bool bRead,
2498 t_symtab *symtab)
2499 {
2500 int j;
2501
2502 for (j = 0; (j < nstr); j++)
2503 {
2504 do_symstr(fio, &(nm[j]), bRead, symtab);
2505 }
2506 }
2507
2508 static void do_resinfo(t_fileio *fio, int n, t_resinfo *ri, gmx_bool bRead,
2509 t_symtab *symtab, int file_version)
2510 {
2511 int j;
2512
2513 for (j = 0; (j < n); j++)
2514 {
2515 do_symstr(fio, &(ri[j].name), bRead, symtab);
2516 if (file_version >= 63)
2517 {
2518 gmx_fio_do_int(fio, ri[j].nr);
2519 gmx_fio_do_uchar(fio, ri[j].ic);
2520 }
2521 else
2522 {
2523 ri[j].nr = j + 1;
2524 ri[j].ic = ' ';
2525 }
2526 }
2527 }
2528
2529 static void do_atoms(t_fileio *fio, t_atoms *atoms, gmx_bool bRead, t_symtab *symtab,
2530 int file_version,
2531 gmx_groups_t *groups)
2532 {
2533 int i;
2534
2535 gmx_fio_do_int(fio, atoms->nr);
2536 gmx_fio_do_int(fio, atoms->nres);
2537 if (file_version < 57)
2538 {
2539 gmx_fio_do_int(fio, groups->ngrpname);
2540 for (i = 0; i < egcNR; i++)
2541 {
2542 groups->ngrpnr[i] = atoms->nr;
2543 snew(groups->grpnr[i], groups->ngrpnr[i]);
2544 }
2545 }
2546 if (bRead)
2547 {
2548 /* Since we have always written all t_atom properties in the tpr file
2549 * (at least for all backward compatible versions), we don't store
2550 * but simple set the booleans here.
2551 */
2552 atoms->haveMass = TRUE;
2553 atoms->haveCharge = TRUE;
2554 atoms->haveType = TRUE;
2555 atoms->haveBState = TRUE;
2556 atoms->havePdbInfo = FALSE;
2557
2558 snew(atoms->atom, atoms->nr);
2559 snew(atoms->atomname, atoms->nr);
2560 snew(atoms->atomtype, atoms->nr);
2561 snew(atoms->atomtypeB, atoms->nr);
2562 snew(atoms->resinfo, atoms->nres);
2563 if (file_version < 57)
2564 {
2565 snew(groups->grpname, groups->ngrpname);
2566 }
2567 atoms->pdbinfo = nullptr;
2568 }
2569 else
2570 {
2571 GMX_RELEASE_ASSERT(atoms->haveMass && atoms->haveCharge && atoms->haveType && atoms->haveBState, "Mass, charge, atomtype and B-state parameters should be present in t_atoms when writing a tpr file");
2572 }
2573 for (i = 0; (i < atoms->nr); i++)
2574 {
2575 do_atom(fio, &atoms->atom[i], egcNR, bRead, file_version, groups, i);
2576 }
2577 do_strstr(fio, atoms->nr, atoms->atomname, bRead, symtab);
2578 do_strstr(fio, atoms->nr, atoms->atomtype, bRead, symtab);
2579 do_strstr(fio, atoms->nr, atoms->atomtypeB, bRead, symtab);
2580
2581 do_resinfo(fio, atoms->nres, atoms->resinfo, bRead, symtab, file_version);
2582
2583 if (file_version < 57)
2584 {
2585 do_strstr(fio, groups->ngrpname, groups->grpname, bRead, symtab);
2586
2587 do_grps(fio, egcNR, groups->grps, bRead, file_version);
2588 }
2589 }
2590
2591 static void do_groups(t_fileio *fio, gmx_groups_t *groups,
2592 gmx_bool bRead, t_symtab *symtab,
2593 int file_version)
2594 {
2595 int g;
2596
2597 do_grps(fio, egcNR, groups->grps, bRead, file_version);
2598 gmx_fio_do_int(fio, groups->ngrpname);
2599 if (bRead)
2600 {
2601 snew(groups->grpname, groups->ngrpname);
2602 }
2603 do_strstr(fio, groups->ngrpname, groups->grpname, bRead, symtab);
2604 for (g = 0; g < egcNR; g++)
2605 {
2606 gmx_fio_do_int(fio, groups->ngrpnr[g]);
2607 if (groups->ngrpnr[g] == 0)
2608 {
2609 if (bRead)
2610 {
2611 groups->grpnr[g] = nullptr;
2612 }
2613 }
2614 else
2615 {
2616 if (bRead)
2617 {
2618 snew(groups->grpnr[g], groups->ngrpnr[g]);
2619 }
2620 gmx_fio_ndo_uchar(fio, groups->grpnr[g], groups->ngrpnr[g]);
2621 }
2622 }
2623 }
2624
2625 static void do_atomtypes(t_fileio *fio, t_atomtypes *atomtypes, gmx_bool bRead,
2626 int file_version)
2627 {
2628 int j;
2629
2630 gmx_fio_do_int(fio, atomtypes->nr);
2631 j = atomtypes->nr;
2632 if (bRead)
2633 {
2634 snew(atomtypes->radius, j);
2635 snew(atomtypes->vol, j);
2636 snew(atomtypes->surftens, j);
2637 snew(atomtypes->atomnumber, j);
2638 snew(atomtypes->gb_radius, j);
2639 snew(atomtypes->S_hct, j);
2640 }
2641 gmx_fio_ndo_real(fio, atomtypes->radius, j);
2642 gmx_fio_ndo_real(fio, atomtypes->vol, j);
2643 gmx_fio_ndo_real(fio, atomtypes->surftens, j);
2644 if (file_version >= 40)
2645 {
2646 gmx_fio_ndo_int(fio, atomtypes->atomnumber, j);
2647 }
2648 if (file_version >= 60)
2649 {
2650 gmx_fio_ndo_real(fio, atomtypes->gb_radius, j);
2651 gmx_fio_ndo_real(fio, atomtypes->S_hct, j);
2652 }
2653 }
2654
2655 static void do_symtab(t_fileio *fio, t_symtab *symtab, gmx_bool bRead)
2656 {
2657 int i, nr;
2658 t_symbuf *symbuf;
2659 char buf[STRLEN];
2660
2661 gmx_fio_do_int(fio, symtab->nr);
2662 nr = symtab->nr;
2663 if (bRead)
2664 {
2665 snew(symtab->symbuf, 1);
2666 symbuf = symtab->symbuf;
2667 symbuf->bufsize = nr;
2668 snew(symbuf->buf, nr);
2669 for (i = 0; (i < nr); i++)
2670 {
2671 gmx_fio_do_string(fio, buf);
2672 symbuf->buf[i] = gmx_strdup(buf);
2673 }
2674 }
2675 else
2676 {
2677 symbuf = symtab->symbuf;
2678 while (symbuf != nullptr)
2679 {
2680 for (i = 0; (i < symbuf->bufsize) && (i < nr); i++)
2681 {
2682 gmx_fio_do_string(fio, symbuf->buf[i]);
2683 }
2684 nr -= i;
2685 symbuf = symbuf->next;
2686 }
2687 if (nr != 0)
2688 {
2689 gmx_fatal(FARGS, "nr of symtab strings left: %d", nr);
2690 }
2691 }
2692 }
2693
2694 static void do_cmap(t_fileio *fio, gmx_cmap_t *cmap_grid, gmx_bool bRead)
2695 {
2696 int i, j, ngrid, gs, nelem;
2697
2698 gmx_fio_do_int(fio, cmap_grid->ngrid);
2699 gmx_fio_do_int(fio, cmap_grid->grid_spacing);
2700
2701 ngrid = cmap_grid->ngrid;
2702 gs = cmap_grid->grid_spacing;
2703 nelem = gs * gs;
2704
2705 if (bRead)
2706 {
2707 snew(cmap_grid->cmapdata, ngrid);
2708
2709 for (i = 0; i < cmap_grid->ngrid; i++)
2710 {
2711 snew(cmap_grid->cmapdata[i].cmap, 4*nelem);
2712 }
2713 }
2714
2715 for (i = 0; i < cmap_grid->ngrid; i++)
2716 {
2717 for (j = 0; j < nelem; j++)
2718 {
2719 gmx_fio_do_real(fio, cmap_grid->cmapdata[i].cmap[j*4]);
2720 gmx_fio_do_real(fio, cmap_grid->cmapdata[i].cmap[j*4+1]);
2721 gmx_fio_do_real(fio, cmap_grid->cmapdata[i].cmap[j*4+2]);
2722 gmx_fio_do_real(fio, cmap_grid->cmapdata[i].cmap[j*4+3]);
2723 }
2724 }
2725 }
2726
2727
2728 static void do_moltype(t_fileio *fio, gmx_moltype_t *molt, gmx_bool bRead,
2729 t_symtab *symtab, int file_version,
2730 gmx_groups_t *groups)
2731 {
2732 if (file_version >= 57)
2733 {
2734 do_symstr(fio, &(molt->name), bRead, symtab);
2735 }
2736
2737 do_atoms(fio, &molt->atoms, bRead, symtab, file_version, groups);
2738
2739 if (file_version >= 57)
2740 {
2741 do_ilists(fio, molt->ilist, bRead, file_version);
2742
2743 do_block(fio, &molt->cgs, bRead, file_version);
2744 }
2745
2746 /* This used to be in the atoms struct */
2747 do_blocka(fio, &molt->excls, bRead, file_version);
2748 }
2749
2750 static void do_molblock(t_fileio *fio, gmx_molblock_t *molb, gmx_bool bRead)
2751 {
2752 gmx_fio_do_int(fio, molb->type);
2753 gmx_fio_do_int(fio, molb->nmol);
2754 gmx_fio_do_int(fio, molb->natoms_mol);
2755 /* Position restraint coordinates */
2756 gmx_fio_do_int(fio, molb->nposres_xA);
2757 if (molb->nposres_xA > 0)
2758 {
2759 if (bRead)
2760 {
2761 snew(molb->posres_xA, molb->nposres_xA);
2762 }
2763 gmx_fio_ndo_rvec(fio, molb->posres_xA, molb->nposres_xA);
2764 }
2765 gmx_fio_do_int(fio, molb->nposres_xB);
2766 if (molb->nposres_xB > 0)
2767 {
2768 if (bRead)
2769 {
2770 snew(molb->posres_xB, molb->nposres_xB);
2771 }
2772 gmx_fio_ndo_rvec(fio, molb->posres_xB, molb->nposres_xB);
2773 }
2774
2775 }
2776
2777 static t_block mtop_mols(gmx_mtop_t *mtop)
2778 {
2779 int mb, m, a, mol;
2780 t_block mols;
2781
2782 mols.nr = 0;
2783 for (mb = 0; mb < mtop->nmolblock; mb++)
2784 {
2785 mols.nr += mtop->molblock[mb].nmol;
2786 }
2787 mols.nalloc_index = mols.nr + 1;
2788 snew(mols.index, mols.nalloc_index);
2789
2790 a = 0;
2791 m = 0;
2792 mols.index[m] = a;
2793 for (mb = 0; mb < mtop->nmolblock; mb++)
2794 {
2795 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
2796 {
2797 a += mtop->molblock[mb].natoms_mol;
2798 m++;
2799 mols.index[m] = a;
2800 }
2801 }
2802
2803 return mols;
2804 }
2805
2806 static void add_posres_molblock(gmx_mtop_t *mtop)
2807 {
2808 t_ilist *il, *ilfb;
2809 int am, i, mol, a;
2810 gmx_bool bFE;
2811 gmx_molblock_t *molb;
2812 t_iparams *ip;
2813
2814 /* posres reference positions are stored in ip->posres (if present) and
2815 in ip->fbposres (if present). If normal and flat-bottomed posres are present,
2816 posres.pos0A are identical to fbposres.pos0. */
2817 il = &mtop->moltype[0].ilist[F_POSRES];
2818 ilfb = &mtop->moltype[0].ilist[F_FBPOSRES];
2819 if (il->nr == 0 && ilfb->nr == 0)
2820 {
2821 return;
2822 }
2823 am = 0;
2824 bFE = FALSE;
2825 for (i = 0; i < il->nr; i += 2)
2826 {
2827 ip = &mtop->ffparams.iparams[il->iatoms[i]];
2828 am = std::max(am, il->iatoms[i+1]);
2829 if (ip->posres.pos0B[XX] != ip->posres.pos0A[XX] ||
2830 ip->posres.pos0B[YY] != ip->posres.pos0A[YY] ||
2831 ip->posres.pos0B[ZZ] != ip->posres.pos0A[ZZ])
2832 {
2833 bFE = TRUE;
2834 }
2835 }
2836 /* This loop is required if we have only flat-bottomed posres:
2837 - set am
2838 - bFE == FALSE (no B-state for flat-bottomed posres) */
2839 if (il->nr == 0)
2840 {
2841 for (i = 0; i < ilfb->nr; i += 2)
2842 {
2843 am = std::max(am, ilfb->iatoms[i+1]);
2844 }
2845 }
2846 /* Make the posres coordinate block end at a molecule end */
2847 mol = 0;
2848 while (am >= mtop->mols.index[mol+1])
2849 {
2850 mol++;
2851 }
2852 molb = &mtop->molblock[0];
2853 molb->nposres_xA = mtop->mols.index[mol+1];
2854 snew(molb->posres_xA, molb->nposres_xA);
2855 if (bFE)
2856 {
2857 molb->nposres_xB = molb->nposres_xA;
2858 snew(molb->posres_xB, molb->nposres_xB);
2859 }
2860 else
2861 {
2862 molb->nposres_xB = 0;
2863 }
2864 for (i = 0; i < il->nr; i += 2)
2865 {
2866 ip = &mtop->ffparams.iparams[il->iatoms[i]];
2867 a = il->iatoms[i+1];
2868 molb->posres_xA[a][XX] = ip->posres.pos0A[XX];
2869 molb->posres_xA[a][YY] = ip->posres.pos0A[YY];
2870 molb->posres_xA[a][ZZ] = ip->posres.pos0A[ZZ];
2871 if (bFE)
2872 {
2873 molb->posres_xB[a][XX] = ip->posres.pos0B[XX];
2874 molb->posres_xB[a][YY] = ip->posres.pos0B[YY];
2875 molb->posres_xB[a][ZZ] = ip->posres.pos0B[ZZ];
2876 }
2877 }
2878 if (il->nr == 0)
2879 {
2880 /* If only flat-bottomed posres are present, take reference pos from them.
2881 Here: bFE == FALSE */
2882 for (i = 0; i < ilfb->nr; i += 2)
2883 {
2884 ip = &mtop->ffparams.iparams[ilfb->iatoms[i]];
2885 a = ilfb->iatoms[i+1];
2886 molb->posres_xA[a][XX] = ip->fbposres.pos0[XX];
2887 molb->posres_xA[a][YY] = ip->fbposres.pos0[YY];
2888 molb->posres_xA[a][ZZ] = ip->fbposres.pos0[ZZ];
2889 }
2890 }
2891 }
2892
2893 static void set_disres_npair(gmx_mtop_t *mtop)
2894 {
2895 t_iparams *ip;
2896 gmx_mtop_ilistloop_t iloop;
2897 t_ilist *ilist, *il;
2898 int nmol, i, npair;
2899 t_iatom *a;
2900
2901 ip = mtop->ffparams.iparams;
2902
2903 iloop = gmx_mtop_ilistloop_init(mtop);
2904 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
2905 {
2906 il = &ilist[F_DISRES];
2907
2908 if (il->nr > 0)
2909 {
2910 a = il->iatoms;
2911 npair = 0;
2912 for (i = 0; i < il->nr; i += 3)
2913 {
2914 npair++;
2915 if (i+3 == il->nr || ip[a[i]].disres.label != ip[a[i+3]].disres.label)
2916 {
2917 ip[a[i]].disres.npair = npair;
2918 npair = 0;
2919 }
2920 }
2921 }
2922 }
2923 }
2924
2925 static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop, gmx_bool bRead,
2926 int file_version)
2927 {
2928 int mt, mb;
2929 t_blocka dumb;
2930
2931 if (bRead)
2932 {
2933 init_mtop(mtop);
2934 }
2935 do_symtab(fio, &(mtop->symtab), bRead);
2936
2937 do_symstr(fio, &(mtop->name), bRead, &(mtop->symtab));
2938
2939 if (file_version >= 57)
2940 {
2941 do_ffparams(fio, &mtop->ffparams, bRead, file_version);
2942
2943 gmx_fio_do_int(fio, mtop->nmoltype);
2944 }
2945 else
2946 {
2947 mtop->nmoltype = 1;
2948 }
2949 if (bRead)
2950 {
2951 snew(mtop->moltype, mtop->nmoltype);
2952 if (file_version < 57)
2953 {
2954 mtop->moltype[0].name = mtop->name;
2955 }
2956 }
2957 for (mt = 0; mt < mtop->nmoltype; mt++)
2958 {
2959 do_moltype(fio, &mtop->moltype[mt], bRead, &mtop->symtab, file_version,
2960 &mtop->groups);
2961 }
2962
2963 if (file_version >= 57)
2964 {
2965 gmx_fio_do_int(fio, mtop->nmolblock);
2966 }
2967 else
2968 {
2969 mtop->nmolblock = 1;
2970 }
2971 if (bRead)
2972 {
2973 snew(mtop->molblock, mtop->nmolblock);
2974 }
2975 if (file_version >= 57)
2976 {
2977 for (mb = 0; mb < mtop->nmolblock; mb++)
2978 {
2979 do_molblock(fio, &mtop->molblock[mb], bRead);
2980 }
2981 gmx_fio_do_int(fio, mtop->natoms);
2982 }
2983 else
2984 {
2985 mtop->molblock[0].type = 0;
2986 mtop->molblock[0].nmol = 1;
2987 mtop->molblock[0].natoms_mol = mtop->moltype[0].atoms.nr;
2988 mtop->molblock[0].nposres_xA = 0;
2989 mtop->molblock[0].nposres_xB = 0;
2990 }
2991
2992 if (file_version >= tpxv_IntermolecularBondeds)
2993 {
2994 gmx_fio_do_gmx_bool(fio, mtop->bIntermolecularInteractions);
2995 if (mtop->bIntermolecularInteractions)
2996 {
2997 if (bRead)
2998 {
2999 snew(mtop->intermolecular_ilist, F_NRE);
3000 }
3001 do_ilists(fio, mtop->intermolecular_ilist, bRead, file_version);
3002 }
3003 }
3004 else
3005 {
3006 mtop->bIntermolecularInteractions = FALSE;
3007 }
3008
3009 do_atomtypes (fio, &(mtop->atomtypes), bRead, file_version);
3010
3011 if (file_version < 57)
3012 {
3013 do_idef (fio, &mtop->ffparams, &mtop->moltype[0], bRead, file_version);
3014 mtop->natoms = mtop->moltype[0].atoms.nr;
3015 }
3016
3017 if (file_version >= 65)
3018 {
3019 do_cmap(fio, &mtop->ffparams.cmap_grid, bRead);
3020 }
3021 else
3022 {
3023 mtop->ffparams.cmap_grid.ngrid = 0;
3024 mtop->ffparams.cmap_grid.grid_spacing = 0;
3025 mtop->ffparams.cmap_grid.cmapdata = nullptr;
3026 }
3027
3028 if (file_version >= 57)
3029 {
3030 do_groups(fio, &mtop->groups, bRead, &(mtop->symtab), file_version);
3031 }
3032
3033 if (file_version < 57)
3034 {
3035 do_block(fio, &mtop->moltype[0].cgs, bRead, file_version);
3036 do_block(fio, &mtop->mols, bRead, file_version);
3037 /* Add the posres coordinates to the molblock */
3038 add_posres_molblock(mtop);
3039 }
3040 if (bRead)
3041 {
3042 if (file_version >= 57)
3043 {
3044 done_block(&mtop->mols);
3045 mtop->mols = mtop_mols(mtop);
3046 }
3047 }
3048
3049 if (file_version < 51)
3050 {
3051 /* Here used to be the shake blocks */
3052 do_blocka(fio, &dumb, bRead, file_version);
3053 if (dumb.nr > 0)
3054 {
3055 sfree(dumb.index);
3056 }
3057 if (dumb.nra > 0)
3058 {
3059 sfree(dumb.a);
3060 }
3061 }
3062
3063 if (bRead)
3064 {
3065 close_symtab(&(mtop->symtab));
3066 }
3067 }
3068
3069 /* If TopOnlyOK is TRUE then we can read even future versions
3070 * of tpx files, provided the fileGeneration hasn't changed.
3071 * If it is FALSE, we need the inputrecord too, and bail out
3072 * if the file is newer than the program.
3073 *
3074 * The version and generation of the topology (see top of this file)
3075 * are returned in the two last arguments, if those arguments are non-NULL.
3076 *
3077 * If possible, we will read the inputrec even when TopOnlyOK is TRUE.
3078 */
3079 static void do_tpxheader(t_fileio *fio, gmx_bool bRead, t_tpxheader *tpx,
3080 gmx_bool TopOnlyOK, int *fileVersionPointer, int *fileGenerationPointer)
3081 {
3082 char buf[STRLEN];
3083 char file_tag[STRLEN];
3084 gmx_bool bDouble;
3085 int precision;
3086 int idum = 0;
3087 real rdum = 0;
3088 int fileVersion; /* Version number of the code that wrote the file */
3089 int fileGeneration; /* Generation version number of the code that wrote the file */
3090
3091 /* XDR binary topology file */
3092 precision = sizeof(real);
3093 if (bRead)
3094 {
3095 gmx_fio_do_string(fio, buf);
3096 if (std::strncmp(buf, "VERSION", 7))
3097 {
3098 gmx_fatal(FARGS, "Can not read file %s,\n"
3099 " this file is from a GROMACS version which is older than 2.0\n"
3100 " Make a new one with grompp or use a gro or pdb file, if possible",
3101 gmx_fio_getname(fio));
3102 }
3103 gmx_fio_do_int(fio, precision);
3104 bDouble = (precision == sizeof(double));
3105 if ((precision != sizeof(float)) && !bDouble)
3106 {
3107 gmx_fatal(FARGS, "Unknown precision in file %s: real is %d bytes "
3108 "instead of %d or %d",
3109 gmx_fio_getname(fio), precision, sizeof(float), sizeof(double));
3110 }
3111 gmx_fio_setprecision(fio, bDouble);
3112 fprintf(stderr, "Reading file %s, %s (%s precision)\n",
3113 gmx_fio_getname(fio), buf, bDouble ? "double" : "single");
3114 }
3115 else
3116 {
3117 snprintf(buf, STRLEN, "VERSION %s", gmx_version());
3118 gmx_fio_write_string(fio, buf);
3119 bDouble = (precision == sizeof(double));
3120 gmx_fio_setprecision(fio, bDouble);
3121 gmx_fio_do_int(fio, precision);
3122 fileVersion = tpx_version;
3123 sprintf(file_tag, "%s", tpx_tag);
3124 fileGeneration = tpx_generation;
3125 }
3126
3127 /* Check versions! */
3128 gmx_fio_do_int(fio, fileVersion);
3129
3130 /* This is for backward compatibility with development versions 77-79
3131 * where the tag was, mistakenly, placed before the generation,
3132 * which would cause a segv instead of a proper error message
3133 * when reading the topology only from tpx with <77 code.
3134 */
3135 if (fileVersion >= 77 && fileVersion <= 79)
3136 {
3137 gmx_fio_do_string(fio, file_tag);
3138 }
3139
3140 gmx_fio_do_int(fio, fileGeneration);
3141
3142 if (fileVersion >= 81)
3143 {
3144 gmx_fio_do_string(fio, file_tag);
3145 }
3146 if (bRead)
3147 {
3148 if (fileVersion < 77)
3149 {
3150 /* Versions before 77 don't have the tag, set it to release */
3151 sprintf(file_tag, "%s", TPX_TAG_RELEASE);
3152 }
3153
3154 if (std::strcmp(file_tag, tpx_tag) != 0)
3155 {
3156 fprintf(stderr, "Note: file tpx tag '%s', software tpx tag '%s'\n",
3157 file_tag, tpx_tag);
3158
3159 /* We only support reading tpx files with the same tag as the code
3160 * or tpx files with the release tag and with lower version number.
3161 */
3162 if (std::strcmp(file_tag, TPX_TAG_RELEASE) != 0 && fileVersion < tpx_version)
3163 {
3164 gmx_fatal(FARGS, "tpx tag/version mismatch: reading tpx file (%s) version %d, tag '%s' with program for tpx version %d, tag '%s'",
3165 gmx_fio_getname(fio), fileVersion, file_tag,
3166 tpx_version, tpx_tag);
3167 }
3168 }
3169 }
3170
3171 if (fileVersionPointer)
3172 {
3173 *fileVersionPointer = fileVersion;
3174 }
3175 if (fileGenerationPointer)
3176 {
3177 *fileGenerationPointer = fileGeneration;
3178 }
3179
3180 if ((fileVersion <= tpx_incompatible_version) ||
3181 ((fileVersion > tpx_version) && !TopOnlyOK) ||
3182 (fileGeneration > tpx_generation) ||
3183 tpx_version == 80) /*80 was used by both 5.0-dev and 4.6-dev*/
3184 {
3185 gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program",
3186 gmx_fio_getname(fio), fileVersion, tpx_version);
3187 }
3188
3189 gmx_fio_do_int(fio, tpx->natoms);
3190 gmx_fio_do_int(fio, tpx->ngtc);
3191
3192 if (fileVersion < 62)
3193 {
3194 gmx_fio_do_int(fio, idum);
3195 gmx_fio_do_real(fio, rdum);
3196 }
3197 if (fileVersion >= 79)
3198 {
3199 gmx_fio_do_int(fio, tpx->fep_state);
3200 }
3201 gmx_fio_do_real(fio, tpx->lambda);
3202 gmx_fio_do_int(fio, tpx->bIr);
3203 gmx_fio_do_int(fio, tpx->bTop);
3204 gmx_fio_do_int(fio, tpx->bX);
3205 gmx_fio_do_int(fio, tpx->bV);
3206 gmx_fio_do_int(fio, tpx->bF);
3207 gmx_fio_do_int(fio, tpx->bBox);
3208
3209 if ((fileGeneration > tpx_generation))
3210 {
3211 /* This can only happen if TopOnlyOK=TRUE */
3212 tpx->bIr = FALSE;
3213 }
3214 }
3215
3216 static int do_tpx(t_fileio *fio, gmx_bool bRead,
3217 t_inputrec *ir, t_state *state, rvec *x, rvec *v,
3218 gmx_mtop_t *mtop)
3219 {
3220 t_tpxheader tpx;
3221 gmx_mtop_t dum_top;
3222 gmx_bool TopOnlyOK;
3223 int ePBC;
3224 gmx_bool bPeriodicMols;
3225
3226 if (!bRead)
3227 {
3228 GMX_RELEASE_ASSERT(state != nullptr, "Cannot write a null state");
3229 GMX_RELEASE_ASSERT(x == nullptr && v == nullptr, "Passing separate x and v pointers to do_tpx() is not supported when writing");
3230
3231 tpx.natoms = state->natoms;
3232 tpx.ngtc = state->ngtc;
3233 tpx.fep_state = state->fep_state;
3234 tpx.lambda = state->lambda[efptFEP];
3235 tpx.bIr = (ir != nullptr);
3236 tpx.bTop = (mtop != nullptr);
3237 tpx.bX = (state->flags & (1 << estX));
3238 tpx.bV = (state->flags & (1 << estV));
3239 tpx.bF = FALSE;
3240 tpx.bBox = TRUE;
3241 }
3242 else
3243 {
3244 GMX_RELEASE_ASSERT(!(x == nullptr && v != nullptr), "Passing x==NULL and v!=NULL is not supported");
3245 }
3246
3247 TopOnlyOK = (ir == nullptr);
3248
3249 int fileVersion; /* Version number of the code that wrote the file */
3250 int fileGeneration; /* Generation version number of the code that wrote the file */
3251 do_tpxheader(fio, bRead, &tpx, TopOnlyOK, &fileVersion, &fileGeneration);
3252
3253 if (bRead)
3254 {
3255 state->flags = 0;
3256 init_gtc_state(state, tpx.ngtc, 0, 0);
3257 if (x == nullptr)
3258 {
3259 // v is also nullptr by the above assertion, so we may
3260 // need to make memory in state for storing the contents
3261 // of the tpx file.
3262 if (tpx.bX)
3263 {
3264 state->flags |= (1 << estX);
3265 }
3266 if (tpx.bV)
3267 {
3268 state->flags |= (1 << estV);
3269 }
3270 state_change_natoms(state, tpx.natoms);
3271 }
3272 }
3273
3274 if (x == nullptr)
3275 {
3276 x = as_rvec_array(state->x.data());
3277 v = as_rvec_array(state->v.data());
3278 }
3279
3280 #define do_test(fio, b, p) if (bRead && (p != NULL) && !b) gmx_fatal(FARGS, "No %s in %s",#p, gmx_fio_getname(fio))
3281
3282 do_test(fio, tpx.bBox, state->box);
3283 if (tpx.bBox)
3284 {
3285 gmx_fio_ndo_rvec(fio, state->box, DIM);
3286 if (fileVersion >= 51)
3287 {
3288 gmx_fio_ndo_rvec(fio, state->box_rel, DIM);
3289 }
3290 else
3291 {
3292 /* We initialize box_rel after reading the inputrec */
3293 clear_mat(state->box_rel);
3294 }
3295 gmx_fio_ndo_rvec(fio, state->boxv, DIM);
3296 if (fileVersion < 56)
3297 {
3298 matrix mdum;
3299 gmx_fio_ndo_rvec(fio, mdum, DIM);
3300 }
3301 }
3302
3303 if (state->ngtc > 0)
3304 {
3305 real *dumv;
3306 snew(dumv, state->ngtc);
3307 if (fileVersion < 69)
3308 {
3309 gmx_fio_ndo_real(fio, dumv, state->ngtc);
3310 }
3311 /* These used to be the Berendsen tcoupl_lambda's */
3312 gmx_fio_ndo_real(fio, dumv, state->ngtc);
3313 sfree(dumv);
3314 }
3315
3316 do_test(fio, tpx.bTop, mtop);
3317 if (tpx.bTop)
3318 {
3319 if (mtop)
3320 {
3321 do_mtop(fio, mtop, bRead, fileVersion);
3322 }
3323 else
3324 {
3325 do_mtop(fio, &dum_top, bRead, fileVersion);
3326 done_mtop(&dum_top);
3327 }
3328 }
3329 do_test(fio, tpx.bX, x);
3330 if (tpx.bX)
3331 {
3332 if (bRead)
3333 {
3334 state->flags |= (1<<estX);
3335 }
3336 gmx_fio_ndo_rvec(fio, x, tpx.natoms);
3337 }
3338
3339 do_test(fio, tpx.bV, v);
3340 if (tpx.bV)
3341 {
3342 if (bRead)
3343 {
3344 state->flags |= (1<<estV);
3345 }
3346 gmx_fio_ndo_rvec(fio, v, tpx.natoms);
3347 }
3348
3349 // No need to run do_test when the last argument is NULL
3350 if (tpx.bF)
3351 {
3352 rvec *dummyForces;
3353 snew(dummyForces, state->natoms);
3354 gmx_fio_ndo_rvec(fio, dummyForces, tpx.natoms);
3355 sfree(dummyForces);
3356 }
3357
3358 /* Starting with tpx version 26, we have the inputrec
3359 * at the end of the file, so we can ignore it
3360 * if the file is never than the software (but still the
3361 * same generation - see comments at the top of this file.
3362 *
3363 *
3364 */
3365 ePBC = -1;
3366 bPeriodicMols = FALSE;
3367
3368 do_test(fio, tpx.bIr, ir);
3369 if (tpx.bIr)
3370 {
3371 if (fileVersion >= 53)
3372 {
3373 /* Removed the pbc info from do_inputrec, since we always want it */
3374 if (!bRead)
3375 {
3376 ePBC = ir->ePBC;
3377 bPeriodicMols = ir->bPeriodicMols;
3378 }
3379 gmx_fio_do_int(fio, ePBC);
3380 gmx_fio_do_gmx_bool(fio, bPeriodicMols);
3381 }
3382 if (fileGeneration <= tpx_generation && ir)
3383 {
3384 do_inputrec(fio, ir, bRead, fileVersion, mtop ? &mtop->ffparams.fudgeQQ : nullptr);
3385 if (fileVersion < 51)
3386 {
3387 set_box_rel(ir, state);
3388 }
3389 if (fileVersion < 53)
3390 {
3391 ePBC = ir->ePBC;
3392 bPeriodicMols = ir->bPeriodicMols;
3393 }
3394 }
3395 if (bRead && ir && fileVersion >= 53)
3396 {
3397 /* We need to do this after do_inputrec, since that initializes ir */
3398 ir->ePBC = ePBC;
3399 ir->bPeriodicMols = bPeriodicMols;
3400 }
3401 }
3402
3403 if (bRead)
3404 {
3405 if (tpx.bIr && ir)
3406 {
3407 if (state->ngtc == 0)
3408 {
3409 /* Reading old version without tcoupl state data: set it */
3410 init_gtc_state(state, ir->opts.ngtc, 0, ir->opts.nhchainlength);
3411 }
3412 if (tpx.bTop && mtop)
3413 {
3414 if (fileVersion < 57)
3415 {
3416 if (mtop->moltype[0].ilist[F_DISRES].nr > 0)
3417 {
3418 ir->eDisre = edrSimple;
3419 }
3420 else
3421 {
3422 ir->eDisre = edrNone;
3423 }
3424 }
3425 set_disres_npair(mtop);
3426 }
3427 }
3428
3429 if (tpx.bTop && mtop)
3430 {
3431 gmx_mtop_finalize(mtop);
3432 }
3433 }
3434
3435 return ePBC;
3436 }
3437
3438 static t_fileio *open_tpx(const char *fn, const char *mode)
3439 {
3440 return gmx_fio_open(fn, mode);
3441 }
3442
3443 static void close_tpx(t_fileio *fio)
3444 {
3445 gmx_fio_close(fio);
3446 }
3447
3448 /************************************************************
3449 *
3450 * The following routines are the exported ones
3451 *
3452 ************************************************************/
3453
3454 void read_tpxheader(const char *fn, t_tpxheader *tpx, gmx_bool TopOnlyOK)
3455 {
3456 t_fileio *fio;
3457
3458 fio = open_tpx(fn, "r");
3459 do_tpxheader(fio, TRUE, tpx, TopOnlyOK, nullptr, nullptr);
3460 close_tpx(fio);
3461 }
3462
3463 void write_tpx_state(const char *fn,
3464 const t_inputrec *ir, const t_state *state, const gmx_mtop_t *mtop)
3465 {
3466 t_fileio *fio;
3467
3468 fio = open_tpx(fn, "w");
3469 do_tpx(fio, FALSE,
3470 const_cast<t_inputrec *>(ir),
3471 const_cast<t_state *>(state), nullptr, nullptr,
3472 const_cast<gmx_mtop_t *>(mtop));
3473 close_tpx(fio);
3474 }
3475
3476 void read_tpx_state(const char *fn,
3477 t_inputrec *ir, t_state *state, gmx_mtop_t *mtop)
3478 {
3479 t_fileio *fio;
3480
3481 fio = open_tpx(fn, "r");
3482 do_tpx(fio, TRUE, ir, state, nullptr, nullptr, mtop);
3483 close_tpx(fio);
3484 }
3485
3486 int read_tpx(const char *fn,
3487 t_inputrec *ir, matrix box, int *natoms,
3488 rvec *x, rvec *v, gmx_mtop_t *mtop)
3489 {
3490 t_fileio *fio;
3491 t_state state;
3492 int ePBC;
3493
3494 fio = open_tpx(fn, "r");
3495 ePBC = do_tpx(fio, TRUE, ir, &state, x, v, mtop);
3496 close_tpx(fio);
3497 if (mtop != nullptr && natoms != nullptr)
3498 {
3499 *natoms = mtop->natoms;
3500 }
3501 if (box)
3502 {
3503 copy_mat(state.box, box);
3504 }
3505
3506 return ePBC;
3507 }
3508
3509 int read_tpx_top(const char *fn,
3510 t_inputrec *ir, matrix box, int *natoms,
3511 rvec *x, rvec *v, t_topology *top)
3512 {
3513 gmx_mtop_t mtop;
3514 int ePBC;
3515
3516 ePBC = read_tpx(fn, ir, box, natoms, x, v, &mtop);
3517
3518 *top = gmx_mtop_t_to_t_topology(&mtop, true);
3519
3520 return ePBC;
3521 }
3522
3523 gmx_bool fn2bTPX(const char *file)
3524 {
3525 return (efTPR == fn2ftp(file));
3526 }
3527
3528 void pr_tpxheader(FILE *fp, int indent, const char *title, const t_tpxheader *sh)
3529 {
3530 if (available(fp, sh, indent, title))
3531 {
3532 indent = pr_title(fp, indent, title);
3533 pr_indent(fp, indent);
3534 fprintf(fp, "bIr = %spresent\n", sh->bIr ? "" : "not ");
3535 pr_indent(fp, indent);
3536 fprintf(fp, "bBox = %spresent\n", sh->bBox ? "" : "not ");
3537 pr_indent(fp, indent);
3538 fprintf(fp, "bTop = %spresent\n", sh->bTop ? "" : "not ");
3539 pr_indent(fp, indent);
3540 fprintf(fp, "bX = %spresent\n", sh->bX ? "" : "not ");
3541 pr_indent(fp, indent);
3542 fprintf(fp, "bV = %spresent\n", sh->bV ? "" : "not ");
3543 pr_indent(fp, indent);
3544 fprintf(fp, "bF = %spresent\n", sh->bF ? "" : "not ");
3545
3546 pr_indent(fp, indent);
3547 fprintf(fp, "natoms = %d\n", sh->natoms);
3548 pr_indent(fp, indent);
3549 fprintf(fp, "lambda = %e\n", sh->lambda);
3550 }
3551 }
The ultimate molecular dynamics simulation package