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