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