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