?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merging in free energy, exp. ensemble, & andersen t-control to 4.6
[gromacs.git] / src / gmxlib / tpxio.c
blob68bcadbc28e6aab97dd1c98051bb3719e395e0c7
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROningen Mixture of Alchemy and Childrens' Stories
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 /* This file is completely threadsafe - keep it that way! */
41 #ifdef GMX_THREAD_MPI
42 #include <thread_mpi.h>
43 #endif
44
45
46 #include <ctype.h>
47 #include "sysstuff.h"
48 #include "smalloc.h"
49 #include "string2.h"
50 #include "gmx_fatal.h"
51 #include "macros.h"
52 #include "names.h"
53 #include "symtab.h"
54 #include "futil.h"
55 #include "filenm.h"
56 #include "gmxfio.h"
57 #include "topsort.h"
58 #include "tpxio.h"
59 #include "txtdump.h"
60 #include "confio.h"
61 #include "atomprop.h"
62 #include "copyrite.h"
63 #include "vec.h"
64 #include "mtop_util.h"
65
66 #define TPX_TAG_RELEASE "release"
67
68 /* This is the tag string which is stored in the tpx file.
69 * Change this if you want to change the tpx format in a feature branch.
70 * This ensures that there will not be different tpx formats around which
71 * can not be distinguished.
72 */
73 static const char *tpx_tag = TPX_TAG_RELEASE;
74
75 /* This number should be increased whenever the file format changes! */
76 static const int tpx_version = 79;
77
78 /* This number should only be increased when you edit the TOPOLOGY section
79 * of the tpx format. This way we can maintain forward compatibility too
80 * for all analysis tools and/or external programs that only need to
81 * know the atom/residue names, charges, and bond connectivity.
82 *
83 * It first appeared in tpx version 26, when I also moved the inputrecord
84 * to the end of the tpx file, so we can just skip it if we only
85 * want the topology.
86 */
87 static const int tpx_generation = 24;
88
89 /* This number should be the most recent backwards incompatible version
90 * I.e., if this number is 9, we cannot read tpx version 9 with this code.
91 */
92 static const int tpx_incompatible_version = 9;
93
94
95
96 /* Struct used to maintain tpx compatibility when function types are added */
97 typedef struct {
98 int fvnr; /* file version number in which the function type first appeared */
99 int ftype; /* function type */
100 } t_ftupd;
101
102 /*
103 *The entries should be ordered in:
104 * 1. ascending file version number
105 * 2. ascending function type number
106 */
107 /*static const t_ftupd ftupd[] = {
108 { 20, F_CUBICBONDS },
109 { 20, F_CONNBONDS },
110 { 20, F_HARMONIC },
111 { 20, F_EQM, },
112 { 22, F_DISRESVIOL },
113 { 22, F_ORIRES },
114 { 22, F_ORIRESDEV },
115 { 26, F_FOURDIHS },
116 { 26, F_PIDIHS },
117 { 26, F_DIHRES },
118 { 26, F_DIHRESVIOL },
119 { 30, F_CROSS_BOND_BONDS },
120 { 30, F_CROSS_BOND_ANGLES },
121 { 30, F_UREY_BRADLEY },
122 { 30, F_POLARIZATION },
123 { 54, F_DHDL_CON },
124 };*/
125 /*
126 *The entries should be ordered in:
127 * 1. ascending function type number
128 * 2. ascending file version number
129 */
130 /* question; what is the purpose of the commented code above? */
131 static const t_ftupd ftupd[] = {
132 { 20, F_CUBICBONDS },
133 { 20, F_CONNBONDS },
134 { 20, F_HARMONIC },
135 { 34, F_FENEBONDS },
136 { 43, F_TABBONDS },
137 { 43, F_TABBONDSNC },
138 { 70, F_RESTRBONDS },
139 { 76, F_LINEAR_ANGLES },
140 { 30, F_CROSS_BOND_BONDS },
141 { 30, F_CROSS_BOND_ANGLES },
142 { 30, F_UREY_BRADLEY },
143 { 34, F_QUARTIC_ANGLES },
144 { 43, F_TABANGLES },
145 { 26, F_FOURDIHS },
146 { 26, F_PIDIHS },
147 { 43, F_TABDIHS },
148 { 65, F_CMAP },
149 { 60, F_GB12 },
150 { 61, F_GB13 },
151 { 61, F_GB14 },
152 { 72, F_GBPOL },
153 { 72, F_NPSOLVATION },
154 { 41, F_LJC14_Q },
155 { 41, F_LJC_PAIRS_NB },
156 { 32, F_BHAM_LR },
157 { 32, F_RF_EXCL },
158 { 32, F_COUL_RECIP },
159 { 46, F_DPD },
160 { 30, F_POLARIZATION },
161 { 36, F_THOLE_POL },
162 { 22, F_DISRESVIOL },
163 { 22, F_ORIRES },
164 { 22, F_ORIRESDEV },
165 { 26, F_DIHRES },
166 { 26, F_DIHRESVIOL },
167 { 49, F_VSITE4FDN },
168 { 50, F_VSITEN },
169 { 46, F_COM_PULL },
170 { 20, F_EQM },
171 { 46, F_ECONSERVED },
172 { 69, F_VTEMP },
173 { 66, F_PDISPCORR },
174 { 54, F_DHDL_CON },
175 { 76, F_ANHARM_POL },
176 { 79, F_DVDL_COUL },
177 { 79, F_DVDL_VDW, },
178 { 79, F_DVDL_BONDED, },
179 { 79, F_DVDL_RESTRAINT },
180 { 79, F_DVDL_TEMPERATURE },
181 { 54, F_DHDL_CON }
182 };
183 #define NFTUPD asize(ftupd)
184
185 /* Needed for backward compatibility */
186 #define MAXNODES 256
187
188 static void _do_section(t_fileio *fio,int key,gmx_bool bRead,const char *src,
189 int line)
190 {
191 char buf[STRLEN];
192 gmx_bool bDbg;
193
194 if (gmx_fio_getftp(fio) == efTPA) {
195 if (!bRead) {
196 gmx_fio_write_string(fio,itemstr[key]);
197 bDbg = gmx_fio_getdebug(fio);
198 gmx_fio_setdebug(fio,FALSE);
199 gmx_fio_write_string(fio,comment_str[key]);
200 gmx_fio_setdebug(fio,bDbg);
201 }
202 else {
203 if (gmx_fio_getdebug(fio))
204 fprintf(stderr,"Looking for section %s (%s, %d)",
205 itemstr[key],src,line);
206
207 do {
208 gmx_fio_do_string(fio,buf);
209 } while ((gmx_strcasecmp(buf,itemstr[key]) != 0));
210
211 if (gmx_strcasecmp(buf,itemstr[key]) != 0)
212 gmx_fatal(FARGS,"\nCould not find section heading %s",itemstr[key]);
213 else if (gmx_fio_getdebug(fio))
214 fprintf(stderr," and found it\n");
215 }
216 }
217 }
218
219 #define do_section(fio,key,bRead) _do_section(fio,key,bRead,__FILE__,__LINE__)
220
221 /**************************************************************
222 *
223 * Now the higer level routines that do io of the structures and arrays
224 *
225 **************************************************************/
226 static void do_pullgrp(t_fileio *fio, t_pullgrp *pgrp, gmx_bool bRead,
227 int file_version)
228 {
229 gmx_bool bDum=TRUE;
230 int i;
231
232 gmx_fio_do_int(fio,pgrp->nat);
233 if (bRead)
234 snew(pgrp->ind,pgrp->nat);
235 bDum=gmx_fio_ndo_int(fio,pgrp->ind,pgrp->nat);
236 gmx_fio_do_int(fio,pgrp->nweight);
237 if (bRead)
238 snew(pgrp->weight,pgrp->nweight);
239 bDum=gmx_fio_ndo_real(fio,pgrp->weight,pgrp->nweight);
240 gmx_fio_do_int(fio,pgrp->pbcatom);
241 gmx_fio_do_rvec(fio,pgrp->vec);
242 gmx_fio_do_rvec(fio,pgrp->init);
243 gmx_fio_do_real(fio,pgrp->rate);
244 gmx_fio_do_real(fio,pgrp->k);
245 if (file_version >= 56) {
246 gmx_fio_do_real(fio,pgrp->kB);
247 } else {
248 pgrp->kB = pgrp->k;
249 }
250 }
251
252 static void do_expandedvals(t_fileio *fio,t_expanded *expand,int n_lambda, gmx_bool bRead, int file_version)
253 {
254 /* i is used in the ndo_double macro*/
255 int i;
256 real fv;
257 gmx_bool bDum=TRUE;
258 real rdum;
259
260 if (file_version >= 79)
261 {
262 if (n_lambda>0)
263 {
264 if (bRead)
265 {
266 snew(expand->init_lambda_weights,n_lambda);
267 }
268 bDum=gmx_fio_ndo_real(fio,expand->init_lambda_weights,n_lambda);
269 gmx_fio_do_gmx_bool(fio,expand->bInit_weights);
270 }
271
272 gmx_fio_do_int(fio,expand->nstexpanded);
273 gmx_fio_do_int(fio,expand->elmcmove);
274 gmx_fio_do_int(fio,expand->elamstats);
275 gmx_fio_do_int(fio,expand->lmc_repeats);
276 gmx_fio_do_int(fio,expand->gibbsdeltalam);
277 gmx_fio_do_int(fio,expand->lmc_forced_nstart);
278 gmx_fio_do_int(fio,expand->lmc_seed);
279 gmx_fio_do_real(fio,expand->mc_temp);
280 gmx_fio_do_int(fio,expand->bSymmetrizedTMatrix);
281 gmx_fio_do_int(fio,expand->nstTij);
282 gmx_fio_do_int(fio,expand->minvarmin);
283 gmx_fio_do_int(fio,expand->c_range);
284 gmx_fio_do_real(fio,expand->wl_scale);
285 gmx_fio_do_real(fio,expand->wl_ratio);
286 gmx_fio_do_real(fio,expand->init_wl_delta);
287 gmx_fio_do_gmx_bool(fio,expand->bWLoneovert);
288 gmx_fio_do_int(fio,expand->elmceq);
289 gmx_fio_do_int(fio,expand->equil_steps);
290 gmx_fio_do_int(fio,expand->equil_samples);
291 gmx_fio_do_int(fio,expand->equil_n_at_lam);
292 gmx_fio_do_real(fio,expand->equil_wl_delta);
293 gmx_fio_do_real(fio,expand->equil_ratio);
294 }
295 }
296
297 static void do_simtempvals(t_fileio *fio,t_simtemp *simtemp, int n_lambda, gmx_bool bRead,
298 int file_version)
299 {
300 gmx_bool bDum=TRUE;
301
302 if (file_version >= 79)
303 {
304 gmx_fio_do_int(fio,simtemp->eSimTempScale);
305 gmx_fio_do_real(fio,simtemp->simtemp_high);
306 gmx_fio_do_real(fio,simtemp->simtemp_low);
307 if (n_lambda>0)
308 {
309 if (bRead)
310 {
311 snew(simtemp->temperatures,n_lambda);
312 }
313 bDum=gmx_fio_ndo_real(fio,simtemp->temperatures,n_lambda);
314 }
315 }
316 }
317
318 static void do_fepvals(t_fileio *fio,t_lambda *fepvals,gmx_bool bRead, int file_version)
319 {
320 /* i is defined in the ndo_double macro; use g to iterate. */
321 int i,g;
322 real fv;
323 gmx_bool bDum=TRUE;
324 real rdum;
325
326 /* free energy values */
327 if (file_version >= 79)
328 {
329 gmx_fio_do_int(fio,fepvals->init_fep_state);
330 gmx_fio_do_double(fio,fepvals->init_lambda);
331 gmx_fio_do_double(fio,fepvals->delta_lambda);
332 }
333 else if (file_version >= 59) {
334 gmx_fio_do_double(fio,fepvals->init_lambda);
335 gmx_fio_do_double(fio,fepvals->delta_lambda);
336 } else {
337 gmx_fio_do_real(fio,rdum);
338 fepvals->init_lambda = rdum;
339 gmx_fio_do_real(fio,rdum);
340 fepvals->delta_lambda = rdum;
341 }
342 if (file_version >= 79)
343 {
344 gmx_fio_do_int(fio,fepvals->n_lambda);
345 if (bRead)
346 {
347 snew(fepvals->all_lambda,efptNR);
348 }
349 for (g=0;g<efptNR;g++)
350 {
351 if (fepvals->n_lambda > 0) {
352 if (bRead)
353 {
354 snew(fepvals->all_lambda[g],fepvals->n_lambda);
355 }
356 bDum=gmx_fio_ndo_double(fio,fepvals->all_lambda[g],fepvals->n_lambda);
357 bDum=gmx_fio_ndo_int(fio,fepvals->separate_dvdl,efptNR);
358 }
359 else if (fepvals->init_lambda >= 0)
360 {
361 fepvals->separate_dvdl[efptFEP] = TRUE;
362 }
363 }
364 }
365 else if (file_version >= 64)
366 {
367 gmx_fio_do_int(fio,fepvals->n_lambda);
368 snew(fepvals->all_lambda,efptNR);
369 if (bRead)
370 {
371 snew(fepvals->all_lambda[efptFEP],fepvals->n_lambda);
372 }
373 bDum=gmx_fio_ndo_double(fio,fepvals->all_lambda[efptFEP],fepvals->n_lambda);
374 if (fepvals->init_lambda >= 0)
375 {
376 fepvals->separate_dvdl[efptFEP] = TRUE;
377 }
378 }
379 else
380 {
381 fepvals->n_lambda = 0;
382 fepvals->all_lambda = NULL;
383 if (fepvals->init_lambda >= 0)
384 {
385 fepvals->separate_dvdl[efptFEP] = TRUE;
386 }
387 }
388 if (file_version >= 13)
389 {
390 gmx_fio_do_real(fio,fepvals->sc_alpha);
391 }
392 else
393 {
394 fepvals->sc_alpha = 0;
395 }
396 if (file_version >= 38)
397 {
398 gmx_fio_do_int(fio,fepvals->sc_power);
399 }
400 else
401 {
402 fepvals->sc_power = 2;
403 }
404 if (file_version >= 79)
405 {
406 gmx_fio_do_real(fio,fepvals->sc_r_power);
407 }
408 else
409 {
410 fepvals->sc_r_power = 6.0;
411 }
412 if (file_version >= 15)
413 {
414 gmx_fio_do_real(fio,fepvals->sc_sigma);
415 }
416 else
417 {
418 fepvals->sc_sigma = 0.3;
419 }
420 if (bRead)
421 {
422 if (file_version >= 71)
423 {
424 fepvals->sc_sigma_min = fepvals->sc_sigma;
425 }
426 else
427 {
428 fepvals->sc_sigma_min = 0;
429 }
430 }
431 if (file_version >= 79)
432 {
433 gmx_fio_do_int(fio,fepvals->bScCoul);
434 }
435 else
436 {
437 fepvals->bScCoul = TRUE;
438 }
439 if (file_version >= 64) {
440 gmx_fio_do_int(fio,fepvals->nstdhdl);
441 } else {
442 fepvals->nstdhdl = 1;
443 }
444
445 if (file_version >= 73)
446 {
447 gmx_fio_do_int(fio, fepvals->separate_dhdl_file);
448 gmx_fio_do_int(fio, fepvals->dhdl_derivatives);
449 }
450 else
451 {
452 fepvals->separate_dhdl_file = esepdhdlfileYES;
453 fepvals->dhdl_derivatives = edhdlderivativesYES;
454 }
455 if (file_version >= 71)
456 {
457 gmx_fio_do_int(fio,fepvals->dh_hist_size);
458 gmx_fio_do_double(fio,fepvals->dh_hist_spacing);
459 }
460 else
461 {
462 fepvals->dh_hist_size = 0;
463 fepvals->dh_hist_spacing = 0.1;
464 }
465 if (file_version >= 79)
466 {
467 gmx_fio_do_int(fio,fepvals->bPrintEnergy);
468 }
469 else
470 {
471 fepvals->bPrintEnergy = FALSE;
472 }
473 }
474
475 static void do_pull(t_fileio *fio, t_pull *pull,gmx_bool bRead, int file_version)
476 {
477 int g;
478
479 gmx_fio_do_int(fio,pull->ngrp);
480 gmx_fio_do_int(fio,pull->eGeom);
481 gmx_fio_do_ivec(fio,pull->dim);
482 gmx_fio_do_real(fio,pull->cyl_r1);
483 gmx_fio_do_real(fio,pull->cyl_r0);
484 gmx_fio_do_real(fio,pull->constr_tol);
485 gmx_fio_do_int(fio,pull->nstxout);
486 gmx_fio_do_int(fio,pull->nstfout);
487 if (bRead)
488 snew(pull->grp,pull->ngrp+1);
489 for(g=0; g<pull->ngrp+1; g++)
490 do_pullgrp(fio,&pull->grp[g],bRead,file_version);
491 }
492
493
494 static void do_rotgrp(t_fileio *fio, t_rotgrp *rotg,gmx_bool bRead, int file_version)
495 {
496 gmx_bool bDum=TRUE;
497 int i;
498
499 gmx_fio_do_int(fio,rotg->eType);
500 gmx_fio_do_int(fio,rotg->bMassW);
501 gmx_fio_do_int(fio,rotg->nat);
502 if (bRead)
503 snew(rotg->ind,rotg->nat);
504 gmx_fio_ndo_int(fio,rotg->ind,rotg->nat);
505 if (bRead)
506 snew(rotg->x_ref,rotg->nat);
507 gmx_fio_ndo_rvec(fio,rotg->x_ref,rotg->nat);
508 gmx_fio_do_rvec(fio,rotg->vec);
509 gmx_fio_do_rvec(fio,rotg->pivot);
510 gmx_fio_do_real(fio,rotg->rate);
511 gmx_fio_do_real(fio,rotg->k);
512 gmx_fio_do_real(fio,rotg->slab_dist);
513 gmx_fio_do_real(fio,rotg->min_gaussian);
514 gmx_fio_do_real(fio,rotg->eps);
515 gmx_fio_do_int(fio,rotg->eFittype);
516 gmx_fio_do_int(fio,rotg->PotAngle_nstep);
517 gmx_fio_do_real(fio,rotg->PotAngle_step);
518 }
519
520 static void do_rot(t_fileio *fio, t_rot *rot,gmx_bool bRead, int file_version)
521 {
522 int g;
523
524 gmx_fio_do_int(fio,rot->ngrp);
525 gmx_fio_do_int(fio,rot->nstrout);
526 gmx_fio_do_int(fio,rot->nstsout);
527 if (bRead)
528 snew(rot->grp,rot->ngrp);
529 for(g=0; g<rot->ngrp; g++)
530 do_rotgrp(fio, &rot->grp[g],bRead,file_version);
531 }
532
533
534 static void do_inputrec(t_fileio *fio, t_inputrec *ir,gmx_bool bRead,
535 int file_version, real *fudgeQQ)
536 {
537 int i,j,k,*tmp,idum=0;
538 gmx_bool bDum=TRUE;
539 real rdum,bd_temp;
540 rvec vdum;
541 gmx_bool bSimAnn;
542 real zerotemptime,finish_t,init_temp,finish_temp;
543
544 if (file_version != tpx_version)
545 {
546 /* Give a warning about features that are not accessible */
547 fprintf(stderr,"Note: file tpx version %d, software tpx version %d\n",
548 file_version,tpx_version);
549 }
550
551 if (bRead)
552 {
553 init_inputrec(ir);
554 }
555
556 if (file_version == 0)
557 {
558 return;
559 }
560
561 /* Basic inputrec stuff */
562 gmx_fio_do_int(fio,ir->eI);
563 if (file_version >= 62) {
564 gmx_fio_do_gmx_large_int(fio, ir->nsteps);
565 } else {
566 gmx_fio_do_int(fio,idum);
567 ir->nsteps = idum;
568 }
569 if(file_version > 25) {
570 if (file_version >= 62) {
571 gmx_fio_do_gmx_large_int(fio, ir->init_step);
572 } else {
573 gmx_fio_do_int(fio,idum);
574 ir->init_step = idum;
575 }
576 } else {
577 ir->init_step=0;
578 }
579
580 if(file_version >= 58)
581 gmx_fio_do_int(fio,ir->simulation_part);
582 else
583 ir->simulation_part=1;
584
585 if (file_version >= 67) {
586 gmx_fio_do_int(fio,ir->nstcalcenergy);
587 } else {
588 ir->nstcalcenergy = 1;
589 }
590 if (file_version < 53) {
591 /* The pbc info has been moved out of do_inputrec,
592 * since we always want it, also without reading the inputrec.
593 */
594 gmx_fio_do_int(fio,ir->ePBC);
595 if ((file_version <= 15) && (ir->ePBC == 2))
596 ir->ePBC = epbcNONE;
597 if (file_version >= 45) {
598 gmx_fio_do_int(fio,ir->bPeriodicMols);
599 } else {
600 if (ir->ePBC == 2) {
601 ir->ePBC = epbcXYZ;
602 ir->bPeriodicMols = TRUE;
603 } else {
604 ir->bPeriodicMols = FALSE;
605 }
606 }
607 }
608 gmx_fio_do_int(fio,ir->ns_type);
609 gmx_fio_do_int(fio,ir->nstlist);
610 gmx_fio_do_int(fio,ir->ndelta);
611 if (file_version < 41) {
612 gmx_fio_do_int(fio,idum);
613 gmx_fio_do_int(fio,idum);
614 }
615 if (file_version >= 45)
616 gmx_fio_do_real(fio,ir->rtpi);
617 else
618 ir->rtpi = 0.05;
619 gmx_fio_do_int(fio,ir->nstcomm);
620 if (file_version > 34)
621 gmx_fio_do_int(fio,ir->comm_mode);
622 else if (ir->nstcomm < 0)
623 ir->comm_mode = ecmANGULAR;
624 else
625 ir->comm_mode = ecmLINEAR;
626 ir->nstcomm = abs(ir->nstcomm);
627
628 if(file_version > 25)
629 gmx_fio_do_int(fio,ir->nstcheckpoint);
630 else
631 ir->nstcheckpoint=0;
632
633 gmx_fio_do_int(fio,ir->nstcgsteep);
634
635 if(file_version>=30)
636 gmx_fio_do_int(fio,ir->nbfgscorr);
637 else if (bRead)
638 ir->nbfgscorr = 10;
639
640 gmx_fio_do_int(fio,ir->nstlog);
641 gmx_fio_do_int(fio,ir->nstxout);
642 gmx_fio_do_int(fio,ir->nstvout);
643 gmx_fio_do_int(fio,ir->nstfout);
644 gmx_fio_do_int(fio,ir->nstenergy);
645 gmx_fio_do_int(fio,ir->nstxtcout);
646 if (file_version >= 59) {
647 gmx_fio_do_double(fio,ir->init_t);
648 gmx_fio_do_double(fio,ir->delta_t);
649 } else {
650 gmx_fio_do_real(fio,rdum);
651 ir->init_t = rdum;
652 gmx_fio_do_real(fio,rdum);
653 ir->delta_t = rdum;
654 }
655 gmx_fio_do_real(fio,ir->xtcprec);
656 if (file_version < 19) {
657 gmx_fio_do_int(fio,idum);
658 gmx_fio_do_int(fio,idum);
659 }
660 if(file_version < 18)
661 gmx_fio_do_int(fio,idum);
662 gmx_fio_do_real(fio,ir->rlist);
663 if (file_version >= 67) {
664 gmx_fio_do_real(fio,ir->rlistlong);
665 }
666 gmx_fio_do_int(fio,ir->coulombtype);
667 if (file_version < 32 && ir->coulombtype == eelRF)
668 ir->coulombtype = eelRF_NEC;
669 gmx_fio_do_real(fio,ir->rcoulomb_switch);
670 gmx_fio_do_real(fio,ir->rcoulomb);
671 gmx_fio_do_int(fio,ir->vdwtype);
672 gmx_fio_do_real(fio,ir->rvdw_switch);
673 gmx_fio_do_real(fio,ir->rvdw);
674 if (file_version < 67) {
675 ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
676 }
677 gmx_fio_do_int(fio,ir->eDispCorr);
678 gmx_fio_do_real(fio,ir->epsilon_r);
679 if (file_version >= 37) {
680 gmx_fio_do_real(fio,ir->epsilon_rf);
681 } else {
682 if (EEL_RF(ir->coulombtype)) {
683 ir->epsilon_rf = ir->epsilon_r;
684 ir->epsilon_r = 1.0;
685 } else {
686 ir->epsilon_rf = 1.0;
687 }
688 }
689 if (file_version >= 29)
690 gmx_fio_do_real(fio,ir->tabext);
691 else
692 ir->tabext=1.0;
693
694 if(file_version > 25) {
695 gmx_fio_do_int(fio,ir->gb_algorithm);
696 gmx_fio_do_int(fio,ir->nstgbradii);
697 gmx_fio_do_real(fio,ir->rgbradii);
698 gmx_fio_do_real(fio,ir->gb_saltconc);
699 gmx_fio_do_int(fio,ir->implicit_solvent);
700 } else {
701 ir->gb_algorithm=egbSTILL;
702 ir->nstgbradii=1;
703 ir->rgbradii=1.0;
704 ir->gb_saltconc=0;
705 ir->implicit_solvent=eisNO;
706 }
707 if(file_version>=55)
708 {
709 gmx_fio_do_real(fio,ir->gb_epsilon_solvent);
710 gmx_fio_do_real(fio,ir->gb_obc_alpha);
711 gmx_fio_do_real(fio,ir->gb_obc_beta);
712 gmx_fio_do_real(fio,ir->gb_obc_gamma);
713 if(file_version>=60)
714 {
715 gmx_fio_do_real(fio,ir->gb_dielectric_offset);
716 gmx_fio_do_int(fio,ir->sa_algorithm);
717 }
718 else
719 {
720 ir->gb_dielectric_offset = 0.009;
721 ir->sa_algorithm = esaAPPROX;
722 }
723 gmx_fio_do_real(fio,ir->sa_surface_tension);
724
725 /* Override sa_surface_tension if it is not changed in the mpd-file */
726 if(ir->sa_surface_tension<0)
727 {
728 if(ir->gb_algorithm==egbSTILL)
729 {
730 ir->sa_surface_tension = 0.0049 * 100 * CAL2JOULE;
731 }
732 else if(ir->gb_algorithm==egbHCT || ir->gb_algorithm==egbOBC)
733 {
734 ir->sa_surface_tension = 0.0054 * 100 * CAL2JOULE;
735 }
736 }
737
738 }
739 else
740 {
741 /* Better use sensible values than insane (0.0) ones... */
742 ir->gb_epsilon_solvent = 80;
743 ir->gb_obc_alpha = 1.0;
744 ir->gb_obc_beta = 0.8;
745 ir->gb_obc_gamma = 4.85;
746 ir->sa_surface_tension = 2.092;
747 }
748
749
750 gmx_fio_do_int(fio,ir->nkx);
751 gmx_fio_do_int(fio,ir->nky);
752 gmx_fio_do_int(fio,ir->nkz);
753 gmx_fio_do_int(fio,ir->pme_order);
754 gmx_fio_do_real(fio,ir->ewald_rtol);
755
756 if (file_version >=24)
757 gmx_fio_do_int(fio,ir->ewald_geometry);
758 else
759 ir->ewald_geometry=eewg3D;
760
761 if (file_version <=17) {
762 ir->epsilon_surface=0;
763 if (file_version==17)
764 gmx_fio_do_int(fio,idum);
765 }
766 else
767 gmx_fio_do_real(fio,ir->epsilon_surface);
768
769 gmx_fio_do_gmx_bool(fio,ir->bOptFFT);
770
771 gmx_fio_do_gmx_bool(fio,ir->bContinuation);
772 gmx_fio_do_int(fio,ir->etc);
773 /* before version 18, ir->etc was a gmx_bool (ir->btc),
774 * but the values 0 and 1 still mean no and
775 * berendsen temperature coupling, respectively.
776 */
777 if (file_version >= 79) {
778 gmx_fio_do_gmx_bool(fio,ir->bPrintNHChains);
779 }
780 if (file_version >= 71)
781 {
782 gmx_fio_do_int(fio,ir->nsttcouple);
783 }
784 else
785 {
786 ir->nsttcouple = ir->nstcalcenergy;
787 }
788 if (file_version <= 15)
789 {
790 gmx_fio_do_int(fio,idum);
791 }
792 if (file_version <=17)
793 {
794 gmx_fio_do_int(fio,ir->epct);
795 if (file_version <= 15)
796 {
797 if (ir->epct == 5)
798 {
799 ir->epct = epctSURFACETENSION;
800 }
801 gmx_fio_do_int(fio,idum);
802 }
803 ir->epct -= 1;
804 /* we have removed the NO alternative at the beginning */
805 if(ir->epct==-1)
806 {
807 ir->epc=epcNO;
808 ir->epct=epctISOTROPIC;
809 }
810 else
811 {
812 ir->epc=epcBERENDSEN;
813 }
814 }
815 else
816 {
817 gmx_fio_do_int(fio,ir->epc);
818 gmx_fio_do_int(fio,ir->epct);
819 }
820 if (file_version >= 71)
821 {
822 gmx_fio_do_int(fio,ir->nstpcouple);
823 }
824 else
825 {
826 ir->nstpcouple = ir->nstcalcenergy;
827 }
828 gmx_fio_do_real(fio,ir->tau_p);
829 if (file_version <= 15) {
830 gmx_fio_do_rvec(fio,vdum);
831 clear_mat(ir->ref_p);
832 for(i=0; i<DIM; i++)
833 ir->ref_p[i][i] = vdum[i];
834 } else {
835 gmx_fio_do_rvec(fio,ir->ref_p[XX]);
836 gmx_fio_do_rvec(fio,ir->ref_p[YY]);
837 gmx_fio_do_rvec(fio,ir->ref_p[ZZ]);
838 }
839 if (file_version <= 15) {
840 gmx_fio_do_rvec(fio,vdum);
841 clear_mat(ir->compress);
842 for(i=0; i<DIM; i++)
843 ir->compress[i][i] = vdum[i];
844 }
845 else {
846 gmx_fio_do_rvec(fio,ir->compress[XX]);
847 gmx_fio_do_rvec(fio,ir->compress[YY]);
848 gmx_fio_do_rvec(fio,ir->compress[ZZ]);
849 }
850 if (file_version >= 47) {
851 gmx_fio_do_int(fio,ir->refcoord_scaling);
852 gmx_fio_do_rvec(fio,ir->posres_com);
853 gmx_fio_do_rvec(fio,ir->posres_comB);
854 } else {
855 ir->refcoord_scaling = erscNO;
856 clear_rvec(ir->posres_com);
857 clear_rvec(ir->posres_comB);
858 }
859 if((file_version > 25) && (file_version < 79))
860 gmx_fio_do_int(fio,ir->andersen_seed);
861 else
862 ir->andersen_seed=0;
863 if(file_version < 26) {
864 gmx_fio_do_gmx_bool(fio,bSimAnn);
865 gmx_fio_do_real(fio,zerotemptime);
866 }
867
868 if (file_version < 37)
869 gmx_fio_do_real(fio,rdum);
870
871 gmx_fio_do_real(fio,ir->shake_tol);
872 if (file_version < 54)
873 gmx_fio_do_real(fio,*fudgeQQ);
874
875 gmx_fio_do_int(fio,ir->efep);
876 if (file_version <= 14 && ir->efep != efepNO)
877 {
878 ir->efep = efepYES;
879 }
880 do_fepvals(fio,ir->fepvals,bRead,file_version);
881
882 if (file_version >= 79)
883 {
884 gmx_fio_do_gmx_bool(fio,ir->bSimTemp);
885 if (ir->bSimTemp)
886 {
887 ir->bSimTemp = TRUE;
888 }
889 }
890 else
891 {
892 ir->bSimTemp = FALSE;
893 }
894 if (ir->bSimTemp)
895 {
896 do_simtempvals(fio,ir->simtempvals,ir->fepvals->n_lambda,bRead,file_version);
897 }
898
899 if (file_version >= 79)
900 {
901 gmx_fio_do_gmx_bool(fio,ir->bExpanded);
902 if (ir->bExpanded)
903 {
904 ir->bExpanded = TRUE;
905 }
906 else
907 {
908 ir->bExpanded = FALSE;
909 }
910 }
911 if (ir->bExpanded)
912 {
913 do_expandedvals(fio,ir->expandedvals,ir->fepvals->n_lambda,bRead,file_version);
914 }
915 if (file_version >= 57) {
916 gmx_fio_do_int(fio,ir->eDisre);
917 }
918 gmx_fio_do_int(fio,ir->eDisreWeighting);
919 if (file_version < 22) {
920 if (ir->eDisreWeighting == 0)
921 ir->eDisreWeighting = edrwEqual;
922 else
923 ir->eDisreWeighting = edrwConservative;
924 }
925 gmx_fio_do_gmx_bool(fio,ir->bDisreMixed);
926 gmx_fio_do_real(fio,ir->dr_fc);
927 gmx_fio_do_real(fio,ir->dr_tau);
928 gmx_fio_do_int(fio,ir->nstdisreout);
929 if (file_version >= 22) {
930 gmx_fio_do_real(fio,ir->orires_fc);
931 gmx_fio_do_real(fio,ir->orires_tau);
932 gmx_fio_do_int(fio,ir->nstorireout);
933 } else {
934 ir->orires_fc = 0;
935 ir->orires_tau = 0;
936 ir->nstorireout = 0;
937 }
938 if(file_version >= 26 && file_version < 79) {
939 gmx_fio_do_real(fio,ir->dihre_fc);
940 if (file_version < 56)
941 {
942 gmx_fio_do_real(fio,rdum);
943 gmx_fio_do_int(fio,idum);
944 }
945 } else {
946 ir->dihre_fc=0;
947 }
948
949 gmx_fio_do_real(fio,ir->em_stepsize);
950 gmx_fio_do_real(fio,ir->em_tol);
951 if (file_version >= 22)
952 gmx_fio_do_gmx_bool(fio,ir->bShakeSOR);
953 else if (bRead)
954 ir->bShakeSOR = TRUE;
955 if (file_version >= 11)
956 gmx_fio_do_int(fio,ir->niter);
957 else if (bRead) {
958 ir->niter = 25;
959 fprintf(stderr,"Note: niter not in run input file, setting it to %d\n",
960 ir->niter);
961 }
962 if (file_version >= 21)
963 gmx_fio_do_real(fio,ir->fc_stepsize);
964 else
965 ir->fc_stepsize = 0;
966 gmx_fio_do_int(fio,ir->eConstrAlg);
967 gmx_fio_do_int(fio,ir->nProjOrder);
968 gmx_fio_do_real(fio,ir->LincsWarnAngle);
969 if (file_version <= 14)
970 gmx_fio_do_int(fio,idum);
971 if (file_version >=26)
972 gmx_fio_do_int(fio,ir->nLincsIter);
973 else if (bRead) {
974 ir->nLincsIter = 1;
975 fprintf(stderr,"Note: nLincsIter not in run input file, setting it to %d\n",
976 ir->nLincsIter);
977 }
978 if (file_version < 33)
979 gmx_fio_do_real(fio,bd_temp);
980 gmx_fio_do_real(fio,ir->bd_fric);
981 gmx_fio_do_int(fio,ir->ld_seed);
982 if (file_version >= 33) {
983 for(i=0; i<DIM; i++)
984 gmx_fio_do_rvec(fio,ir->deform[i]);
985 } else {
986 for(i=0; i<DIM; i++)
987 clear_rvec(ir->deform[i]);
988 }
989 if (file_version >= 14)
990 gmx_fio_do_real(fio,ir->cos_accel);
991 else if (bRead)
992 ir->cos_accel = 0;
993 gmx_fio_do_int(fio,ir->userint1);
994 gmx_fio_do_int(fio,ir->userint2);
995 gmx_fio_do_int(fio,ir->userint3);
996 gmx_fio_do_int(fio,ir->userint4);
997 gmx_fio_do_real(fio,ir->userreal1);
998 gmx_fio_do_real(fio,ir->userreal2);
999 gmx_fio_do_real(fio,ir->userreal3);
1000 gmx_fio_do_real(fio,ir->userreal4);
1001
1002 /* AdResS stuff */
1003 if (file_version >= 77) {
1004 gmx_fio_do_gmx_bool(fio,ir->bAdress);
1005 if(ir->bAdress){
1006 if (bRead) snew(ir->adress, 1);
1007 gmx_fio_do_int(fio,ir->adress->type);
1008 gmx_fio_do_real(fio,ir->adress->const_wf);
1009 gmx_fio_do_real(fio,ir->adress->ex_width);
1010 gmx_fio_do_real(fio,ir->adress->hy_width);
1011 gmx_fio_do_int(fio,ir->adress->icor);
1012 gmx_fio_do_int(fio,ir->adress->site);
1013 gmx_fio_do_rvec(fio,ir->adress->refs);
1014 gmx_fio_do_int(fio,ir->adress->n_tf_grps);
1015 gmx_fio_do_real(fio, ir->adress->ex_forcecap);
1016 gmx_fio_do_int(fio, ir->adress->n_energy_grps);
1017 gmx_fio_do_int(fio,ir->adress->do_hybridpairs);
1018
1019 if (bRead)snew(ir->adress->tf_table_index,ir->adress->n_tf_grps);
1020 if (ir->adress->n_tf_grps > 0) {
1021 bDum=gmx_fio_ndo_int(fio,ir->adress->tf_table_index,ir->adress->n_tf_grps);
1022 }
1023 if (bRead)snew(ir->adress->group_explicit,ir->adress->n_energy_grps);
1024 if (ir->adress->n_energy_grps > 0) {
1025 bDum=gmx_fio_ndo_int(fio, ir->adress->group_explicit,ir->adress->n_energy_grps);
1026 }
1027 }
1028 } else {
1029 ir->bAdress = FALSE;
1030 }
1031
1032 /* pull stuff */
1033 if (file_version >= 48) {
1034 gmx_fio_do_int(fio,ir->ePull);
1035 if (ir->ePull != epullNO) {
1036 if (bRead)
1037 snew(ir->pull,1);
1038 do_pull(fio, ir->pull,bRead,file_version);
1039 }
1040 } else {
1041 ir->ePull = epullNO;
1042 }
1043
1044 /* Enforced rotation */
1045 if (file_version >= 74) {
1046 gmx_fio_do_int(fio,ir->bRot);
1047 if (ir->bRot == TRUE) {
1048 if (bRead)
1049 snew(ir->rot,1);
1050 do_rot(fio, ir->rot,bRead,file_version);
1051 }
1052 } else {
1053 ir->bRot = FALSE;
1054 }
1055
1056 /* grpopts stuff */
1057 gmx_fio_do_int(fio,ir->opts.ngtc);
1058 if (file_version >= 69) {
1059 gmx_fio_do_int(fio,ir->opts.nhchainlength);
1060 } else {
1061 ir->opts.nhchainlength = 1;
1062 }
1063 gmx_fio_do_int(fio,ir->opts.ngacc);
1064 gmx_fio_do_int(fio,ir->opts.ngfrz);
1065 gmx_fio_do_int(fio,ir->opts.ngener);
1066
1067 if (bRead) {
1068 snew(ir->opts.nrdf, ir->opts.ngtc);
1069 snew(ir->opts.ref_t, ir->opts.ngtc);
1070 snew(ir->opts.annealing, ir->opts.ngtc);
1071 snew(ir->opts.anneal_npoints, ir->opts.ngtc);
1072 snew(ir->opts.anneal_time, ir->opts.ngtc);
1073 snew(ir->opts.anneal_temp, ir->opts.ngtc);
1074 snew(ir->opts.tau_t, ir->opts.ngtc);
1075 snew(ir->opts.nFreeze,ir->opts.ngfrz);
1076 snew(ir->opts.acc, ir->opts.ngacc);
1077 snew(ir->opts.egp_flags,ir->opts.ngener*ir->opts.ngener);
1078 }
1079 if (ir->opts.ngtc > 0) {
1080 if (bRead && file_version<13) {
1081 snew(tmp,ir->opts.ngtc);
1082 bDum=gmx_fio_ndo_int(fio,tmp, ir->opts.ngtc);
1083 for(i=0; i<ir->opts.ngtc; i++)
1084 ir->opts.nrdf[i] = tmp[i];
1085 sfree(tmp);
1086 } else {
1087 bDum=gmx_fio_ndo_real(fio,ir->opts.nrdf, ir->opts.ngtc);
1088 }
1089 bDum=gmx_fio_ndo_real(fio,ir->opts.ref_t,ir->opts.ngtc);
1090 bDum=gmx_fio_ndo_real(fio,ir->opts.tau_t,ir->opts.ngtc);
1091 if (file_version<33 && ir->eI==eiBD) {
1092 for(i=0; i<ir->opts.ngtc; i++)
1093 ir->opts.tau_t[i] = bd_temp;
1094 }
1095 }
1096 if (ir->opts.ngfrz > 0)
1097 bDum=gmx_fio_ndo_ivec(fio,ir->opts.nFreeze,ir->opts.ngfrz);
1098 if (ir->opts.ngacc > 0)
1099 gmx_fio_ndo_rvec(fio,ir->opts.acc,ir->opts.ngacc);
1100 if (file_version >= 12)
1101 bDum=gmx_fio_ndo_int(fio,ir->opts.egp_flags,
1102 ir->opts.ngener*ir->opts.ngener);
1103
1104 if(bRead && file_version < 26) {
1105 for(i=0;i<ir->opts.ngtc;i++) {
1106 if(bSimAnn) {
1107 ir->opts.annealing[i] = eannSINGLE;
1108 ir->opts.anneal_npoints[i] = 2;
1109 snew(ir->opts.anneal_time[i],2);
1110 snew(ir->opts.anneal_temp[i],2);
1111 /* calculate the starting/ending temperatures from reft, zerotemptime, and nsteps */
1112 finish_t = ir->init_t + ir->nsteps * ir->delta_t;
1113 init_temp = ir->opts.ref_t[i]*(1-ir->init_t/zerotemptime);
1114 finish_temp = ir->opts.ref_t[i]*(1-finish_t/zerotemptime);
1115 ir->opts.anneal_time[i][0] = ir->init_t;
1116 ir->opts.anneal_time[i][1] = finish_t;
1117 ir->opts.anneal_temp[i][0] = init_temp;
1118 ir->opts.anneal_temp[i][1] = finish_temp;
1119 } else {
1120 ir->opts.annealing[i] = eannNO;
1121 ir->opts.anneal_npoints[i] = 0;
1122 }
1123 }
1124 } else {
1125 /* file version 26 or later */
1126 /* First read the lists with annealing and npoints for each group */
1127 bDum=gmx_fio_ndo_int(fio,ir->opts.annealing,ir->opts.ngtc);
1128 bDum=gmx_fio_ndo_int(fio,ir->opts.anneal_npoints,ir->opts.ngtc);
1129 for(j=0;j<(ir->opts.ngtc);j++) {
1130 k=ir->opts.anneal_npoints[j];
1131 if(bRead) {
1132 snew(ir->opts.anneal_time[j],k);
1133 snew(ir->opts.anneal_temp[j],k);
1134 }
1135 bDum=gmx_fio_ndo_real(fio,ir->opts.anneal_time[j],k);
1136 bDum=gmx_fio_ndo_real(fio,ir->opts.anneal_temp[j],k);
1137 }
1138 }
1139 /* Walls */
1140 if (file_version >= 45) {
1141 gmx_fio_do_int(fio,ir->nwall);
1142 gmx_fio_do_int(fio,ir->wall_type);
1143 if (file_version >= 50)
1144 gmx_fio_do_real(fio,ir->wall_r_linpot);
1145 else
1146 ir->wall_r_linpot = -1;
1147 gmx_fio_do_int(fio,ir->wall_atomtype[0]);
1148 gmx_fio_do_int(fio,ir->wall_atomtype[1]);
1149 gmx_fio_do_real(fio,ir->wall_density[0]);
1150 gmx_fio_do_real(fio,ir->wall_density[1]);
1151 gmx_fio_do_real(fio,ir->wall_ewald_zfac);
1152 } else {
1153 ir->nwall = 0;
1154 ir->wall_type = 0;
1155 ir->wall_atomtype[0] = -1;
1156 ir->wall_atomtype[1] = -1;
1157 ir->wall_density[0] = 0;
1158 ir->wall_density[1] = 0;
1159 ir->wall_ewald_zfac = 3;
1160 }
1161 /* Cosine stuff for electric fields */
1162 for(j=0; (j<DIM); j++) {
1163 gmx_fio_do_int(fio,ir->ex[j].n);
1164 gmx_fio_do_int(fio,ir->et[j].n);
1165 if (bRead) {
1166 snew(ir->ex[j].a, ir->ex[j].n);
1167 snew(ir->ex[j].phi,ir->ex[j].n);
1168 snew(ir->et[j].a, ir->et[j].n);
1169 snew(ir->et[j].phi,ir->et[j].n);
1170 }
1171 bDum=gmx_fio_ndo_real(fio,ir->ex[j].a, ir->ex[j].n);
1172 bDum=gmx_fio_ndo_real(fio,ir->ex[j].phi,ir->ex[j].n);
1173 bDum=gmx_fio_ndo_real(fio,ir->et[j].a, ir->et[j].n);
1174 bDum=gmx_fio_ndo_real(fio,ir->et[j].phi,ir->et[j].n);
1175 }
1176
1177 /* QMMM stuff */
1178 if(file_version>=39){
1179 gmx_fio_do_gmx_bool(fio,ir->bQMMM);
1180 gmx_fio_do_int(fio,ir->QMMMscheme);
1181 gmx_fio_do_real(fio,ir->scalefactor);
1182 gmx_fio_do_int(fio,ir->opts.ngQM);
1183 if (bRead) {
1184 snew(ir->opts.QMmethod, ir->opts.ngQM);
1185 snew(ir->opts.QMbasis, ir->opts.ngQM);
1186 snew(ir->opts.QMcharge, ir->opts.ngQM);
1187 snew(ir->opts.QMmult, ir->opts.ngQM);
1188 snew(ir->opts.bSH, ir->opts.ngQM);
1189 snew(ir->opts.CASorbitals, ir->opts.ngQM);
1190 snew(ir->opts.CASelectrons,ir->opts.ngQM);
1191 snew(ir->opts.SAon, ir->opts.ngQM);
1192 snew(ir->opts.SAoff, ir->opts.ngQM);
1193 snew(ir->opts.SAsteps, ir->opts.ngQM);
1194 snew(ir->opts.bOPT, ir->opts.ngQM);
1195 snew(ir->opts.bTS, ir->opts.ngQM);
1196 }
1197 if (ir->opts.ngQM > 0) {
1198 bDum=gmx_fio_ndo_int(fio,ir->opts.QMmethod,ir->opts.ngQM);
1199 bDum=gmx_fio_ndo_int(fio,ir->opts.QMbasis,ir->opts.ngQM);
1200 bDum=gmx_fio_ndo_int(fio,ir->opts.QMcharge,ir->opts.ngQM);
1201 bDum=gmx_fio_ndo_int(fio,ir->opts.QMmult,ir->opts.ngQM);
1202 bDum=gmx_fio_ndo_gmx_bool(fio,ir->opts.bSH,ir->opts.ngQM);
1203 bDum=gmx_fio_ndo_int(fio,ir->opts.CASorbitals,ir->opts.ngQM);
1204 bDum=gmx_fio_ndo_int(fio,ir->opts.CASelectrons,ir->opts.ngQM);
1205 bDum=gmx_fio_ndo_real(fio,ir->opts.SAon,ir->opts.ngQM);
1206 bDum=gmx_fio_ndo_real(fio,ir->opts.SAoff,ir->opts.ngQM);
1207 bDum=gmx_fio_ndo_int(fio,ir->opts.SAsteps,ir->opts.ngQM);
1208 bDum=gmx_fio_ndo_gmx_bool(fio,ir->opts.bOPT,ir->opts.ngQM);
1209 bDum=gmx_fio_ndo_gmx_bool(fio,ir->opts.bTS,ir->opts.ngQM);
1210 }
1211 /* end of QMMM stuff */
1212 }
1213 }
1214
1215
1216 static void do_harm(t_fileio *fio, t_iparams *iparams,gmx_bool bRead)
1217 {
1218 gmx_fio_do_real(fio,iparams->harmonic.rA);
1219 gmx_fio_do_real(fio,iparams->harmonic.krA);
1220 gmx_fio_do_real(fio,iparams->harmonic.rB);
1221 gmx_fio_do_real(fio,iparams->harmonic.krB);
1222 }
1223
1224 void do_iparams(t_fileio *fio, t_functype ftype,t_iparams *iparams,
1225 gmx_bool bRead, int file_version)
1226 {
1227 int i;
1228 gmx_bool bDum;
1229 real rdum;
1230
1231 if (!bRead)
1232 gmx_fio_set_comment(fio, interaction_function[ftype].name);
1233 switch (ftype) {
1234 case F_ANGLES:
1235 case F_G96ANGLES:
1236 case F_BONDS:
1237 case F_G96BONDS:
1238 case F_HARMONIC:
1239 case F_IDIHS:
1240 do_harm(fio, iparams,bRead);
1241 if ((ftype == F_ANGRES || ftype == F_ANGRESZ) && bRead) {
1242 /* Correct incorrect storage of parameters */
1243 iparams->pdihs.phiB = iparams->pdihs.phiA;
1244 iparams->pdihs.cpB = iparams->pdihs.cpA;
1245 }
1246 break;
1247 case F_LINEAR_ANGLES:
1248 gmx_fio_do_real(fio,iparams->linangle.klinA);
1249 gmx_fio_do_real(fio,iparams->linangle.aA);
1250 gmx_fio_do_real(fio,iparams->linangle.klinB);
1251 gmx_fio_do_real(fio,iparams->linangle.aB);
1252 break;
1253 case F_FENEBONDS:
1254 gmx_fio_do_real(fio,iparams->fene.bm);
1255 gmx_fio_do_real(fio,iparams->fene.kb);
1256 break;
1257 case F_RESTRBONDS:
1258 gmx_fio_do_real(fio,iparams->restraint.lowA);
1259 gmx_fio_do_real(fio,iparams->restraint.up1A);
1260 gmx_fio_do_real(fio,iparams->restraint.up2A);
1261 gmx_fio_do_real(fio,iparams->restraint.kA);
1262 gmx_fio_do_real(fio,iparams->restraint.lowB);
1263 gmx_fio_do_real(fio,iparams->restraint.up1B);
1264 gmx_fio_do_real(fio,iparams->restraint.up2B);
1265 gmx_fio_do_real(fio,iparams->restraint.kB);
1266 break;
1267 case F_TABBONDS:
1268 case F_TABBONDSNC:
1269 case F_TABANGLES:
1270 case F_TABDIHS:
1271 gmx_fio_do_real(fio,iparams->tab.kA);
1272 gmx_fio_do_int(fio,iparams->tab.table);
1273 gmx_fio_do_real(fio,iparams->tab.kB);
1274 break;
1275 case F_CROSS_BOND_BONDS:
1276 gmx_fio_do_real(fio,iparams->cross_bb.r1e);
1277 gmx_fio_do_real(fio,iparams->cross_bb.r2e);
1278 gmx_fio_do_real(fio,iparams->cross_bb.krr);
1279 break;
1280 case F_CROSS_BOND_ANGLES:
1281 gmx_fio_do_real(fio,iparams->cross_ba.r1e);
1282 gmx_fio_do_real(fio,iparams->cross_ba.r2e);
1283 gmx_fio_do_real(fio,iparams->cross_ba.r3e);
1284 gmx_fio_do_real(fio,iparams->cross_ba.krt);
1285 break;
1286 case F_UREY_BRADLEY:
1287 gmx_fio_do_real(fio,iparams->u_b.thetaA);
1288 gmx_fio_do_real(fio,iparams->u_b.kthetaA);
1289 gmx_fio_do_real(fio,iparams->u_b.r13A);
1290 gmx_fio_do_real(fio,iparams->u_b.kUBA);
1291 if (file_version >= 79) {
1292 gmx_fio_do_real(fio,iparams->u_b.thetaB);
1293 gmx_fio_do_real(fio,iparams->u_b.kthetaB);
1294 gmx_fio_do_real(fio,iparams->u_b.r13B);
1295 gmx_fio_do_real(fio,iparams->u_b.kUBB);
1296 } else {
1297 iparams->u_b.thetaB=iparams->u_b.thetaA;
1298 iparams->u_b.kthetaB=iparams->u_b.kthetaA;
1299 iparams->u_b.r13B=iparams->u_b.r13A;
1300 iparams->u_b.kUBB=iparams->u_b.kUBA;
1301 }
1302 break;
1303 case F_QUARTIC_ANGLES:
1304 gmx_fio_do_real(fio,iparams->qangle.theta);
1305 bDum=gmx_fio_ndo_real(fio,iparams->qangle.c,5);
1306 break;
1307 case F_BHAM:
1308 gmx_fio_do_real(fio,iparams->bham.a);
1309 gmx_fio_do_real(fio,iparams->bham.b);
1310 gmx_fio_do_real(fio,iparams->bham.c);
1311 break;
1312 case F_MORSE:
1313 gmx_fio_do_real(fio,iparams->morse.b0A);
1314 gmx_fio_do_real(fio,iparams->morse.cbA);
1315 gmx_fio_do_real(fio,iparams->morse.betaA);
1316 if (file_version >= 79) {
1317 gmx_fio_do_real(fio,iparams->morse.b0B);
1318 gmx_fio_do_real(fio,iparams->morse.cbB);
1319 gmx_fio_do_real(fio,iparams->morse.betaB);
1320 } else {
1321 iparams->morse.b0B = iparams->morse.b0A;
1322 iparams->morse.cbB = iparams->morse.cbA;
1323 iparams->morse.betaB = iparams->morse.betaA;
1324 }
1325 break;
1326 case F_CUBICBONDS:
1327 gmx_fio_do_real(fio,iparams->cubic.b0);
1328 gmx_fio_do_real(fio,iparams->cubic.kb);
1329 gmx_fio_do_real(fio,iparams->cubic.kcub);
1330 break;
1331 case F_CONNBONDS:
1332 break;
1333 case F_POLARIZATION:
1334 gmx_fio_do_real(fio,iparams->polarize.alpha);
1335 break;
1336 case F_ANHARM_POL:
1337 gmx_fio_do_real(fio,iparams->anharm_polarize.alpha);
1338 gmx_fio_do_real(fio,iparams->anharm_polarize.drcut);
1339 gmx_fio_do_real(fio,iparams->anharm_polarize.khyp);
1340 break;
1341 case F_WATER_POL:
1342 if (file_version < 31)
1343 gmx_fatal(FARGS,"Old tpr files with water_polarization not supported. Make a new.");
1344 gmx_fio_do_real(fio,iparams->wpol.al_x);
1345 gmx_fio_do_real(fio,iparams->wpol.al_y);
1346 gmx_fio_do_real(fio,iparams->wpol.al_z);
1347 gmx_fio_do_real(fio,iparams->wpol.rOH);
1348 gmx_fio_do_real(fio,iparams->wpol.rHH);
1349 gmx_fio_do_real(fio,iparams->wpol.rOD);
1350 break;
1351 case F_THOLE_POL:
1352 gmx_fio_do_real(fio,iparams->thole.a);
1353 gmx_fio_do_real(fio,iparams->thole.alpha1);
1354 gmx_fio_do_real(fio,iparams->thole.alpha2);
1355 gmx_fio_do_real(fio,iparams->thole.rfac);
1356 break;
1357 case F_LJ:
1358 gmx_fio_do_real(fio,iparams->lj.c6);
1359 gmx_fio_do_real(fio,iparams->lj.c12);
1360 break;
1361 case F_LJ14:
1362 gmx_fio_do_real(fio,iparams->lj14.c6A);
1363 gmx_fio_do_real(fio,iparams->lj14.c12A);
1364 gmx_fio_do_real(fio,iparams->lj14.c6B);
1365 gmx_fio_do_real(fio,iparams->lj14.c12B);
1366 break;
1367 case F_LJC14_Q:
1368 gmx_fio_do_real(fio,iparams->ljc14.fqq);
1369 gmx_fio_do_real(fio,iparams->ljc14.qi);
1370 gmx_fio_do_real(fio,iparams->ljc14.qj);
1371 gmx_fio_do_real(fio,iparams->ljc14.c6);
1372 gmx_fio_do_real(fio,iparams->ljc14.c12);
1373 break;
1374 case F_LJC_PAIRS_NB:
1375 gmx_fio_do_real(fio,iparams->ljcnb.qi);
1376 gmx_fio_do_real(fio,iparams->ljcnb.qj);
1377 gmx_fio_do_real(fio,iparams->ljcnb.c6);
1378 gmx_fio_do_real(fio,iparams->ljcnb.c12);
1379 break;
1380 case F_PDIHS:
1381 case F_PIDIHS:
1382 case F_ANGRES:
1383 case F_ANGRESZ:
1384 gmx_fio_do_real(fio,iparams->pdihs.phiA);
1385 gmx_fio_do_real(fio,iparams->pdihs.cpA);
1386 if ((ftype == F_ANGRES || ftype == F_ANGRESZ) && file_version < 42) {
1387 /* Read the incorrectly stored multiplicity */
1388 gmx_fio_do_real(fio,iparams->harmonic.rB);
1389 gmx_fio_do_real(fio,iparams->harmonic.krB);
1390 iparams->pdihs.phiB = iparams->pdihs.phiA;
1391 iparams->pdihs.cpB = iparams->pdihs.cpA;
1392 } else {
1393 gmx_fio_do_real(fio,iparams->pdihs.phiB);
1394 gmx_fio_do_real(fio,iparams->pdihs.cpB);
1395 gmx_fio_do_int(fio,iparams->pdihs.mult);
1396 }
1397 break;
1398 case F_DISRES:
1399 gmx_fio_do_int(fio,iparams->disres.label);
1400 gmx_fio_do_int(fio,iparams->disres.type);
1401 gmx_fio_do_real(fio,iparams->disres.low);
1402 gmx_fio_do_real(fio,iparams->disres.up1);
1403 gmx_fio_do_real(fio,iparams->disres.up2);
1404 gmx_fio_do_real(fio,iparams->disres.kfac);
1405 break;
1406 case F_ORIRES:
1407 gmx_fio_do_int(fio,iparams->orires.ex);
1408 gmx_fio_do_int(fio,iparams->orires.label);
1409 gmx_fio_do_int(fio,iparams->orires.power);
1410 gmx_fio_do_real(fio,iparams->orires.c);
1411 gmx_fio_do_real(fio,iparams->orires.obs);
1412 gmx_fio_do_real(fio,iparams->orires.kfac);
1413 break;
1414 case F_DIHRES:
1415 gmx_fio_do_real(fio,iparams->dihres.phiA);
1416 gmx_fio_do_real(fio,iparams->dihres.dphiA);
1417 gmx_fio_do_real(fio,iparams->dihres.kfacA);
1418 if (file_version >= 72) {
1419 gmx_fio_do_real(fio,iparams->dihres.phiB);
1420 gmx_fio_do_real(fio,iparams->dihres.dphiB);
1421 gmx_fio_do_real(fio,iparams->dihres.kfacB);
1422 } else {
1423 iparams->dihres.phiB=iparams->dihres.phiA;
1424 iparams->dihres.dphiB=iparams->dihres.dphiA;
1425 iparams->dihres.kfacB=iparams->dihres.kfacA;
1426 }
1427 break;
1428 case F_POSRES:
1429 gmx_fio_do_rvec(fio,iparams->posres.pos0A);
1430 gmx_fio_do_rvec(fio,iparams->posres.fcA);
1431 if (bRead && file_version < 27) {
1432 copy_rvec(iparams->posres.pos0A,iparams->posres.pos0B);
1433 copy_rvec(iparams->posres.fcA,iparams->posres.fcB);
1434 } else {
1435 gmx_fio_do_rvec(fio,iparams->posres.pos0B);
1436 gmx_fio_do_rvec(fio,iparams->posres.fcB);
1437 }
1438 break;
1439 case F_RBDIHS:
1440 bDum=gmx_fio_ndo_real(fio,iparams->rbdihs.rbcA,NR_RBDIHS);
1441 if(file_version>=25)
1442 bDum=gmx_fio_ndo_real(fio,iparams->rbdihs.rbcB,NR_RBDIHS);
1443 break;
1444 case F_FOURDIHS:
1445 /* Fourier dihedrals are internally represented
1446 * as Ryckaert-Bellemans since those are faster to compute.
1447 */
1448 bDum=gmx_fio_ndo_real(fio,iparams->rbdihs.rbcA, NR_RBDIHS);
1449 bDum=gmx_fio_ndo_real(fio,iparams->rbdihs.rbcB, NR_RBDIHS);
1450 break;
1451 case F_CONSTR:
1452 case F_CONSTRNC:
1453 gmx_fio_do_real(fio,iparams->constr.dA);
1454 gmx_fio_do_real(fio,iparams->constr.dB);
1455 break;
1456 case F_SETTLE:
1457 gmx_fio_do_real(fio,iparams->settle.doh);
1458 gmx_fio_do_real(fio,iparams->settle.dhh);
1459 break;
1460 case F_VSITE2:
1461 gmx_fio_do_real(fio,iparams->vsite.a);
1462 break;
1463 case F_VSITE3:
1464 case F_VSITE3FD:
1465 case F_VSITE3FAD:
1466 gmx_fio_do_real(fio,iparams->vsite.a);
1467 gmx_fio_do_real(fio,iparams->vsite.b);
1468 break;
1469 case F_VSITE3OUT:
1470 case F_VSITE4FD:
1471 case F_VSITE4FDN:
1472 gmx_fio_do_real(fio,iparams->vsite.a);
1473 gmx_fio_do_real(fio,iparams->vsite.b);
1474 gmx_fio_do_real(fio,iparams->vsite.c);
1475 break;
1476 case F_VSITEN:
1477 gmx_fio_do_int(fio,iparams->vsiten.n);
1478 gmx_fio_do_real(fio,iparams->vsiten.a);
1479 break;
1480 case F_GB12:
1481 case F_GB13:
1482 case F_GB14:
1483 /* We got rid of some parameters in version 68 */
1484 if(bRead && file_version<68)
1485 {
1486 gmx_fio_do_real(fio,rdum);
1487 gmx_fio_do_real(fio,rdum);
1488 gmx_fio_do_real(fio,rdum);
1489 gmx_fio_do_real(fio,rdum);
1490 }
1491 gmx_fio_do_real(fio,iparams->gb.sar);
1492 gmx_fio_do_real(fio,iparams->gb.st);
1493 gmx_fio_do_real(fio,iparams->gb.pi);
1494 gmx_fio_do_real(fio,iparams->gb.gbr);
1495 gmx_fio_do_real(fio,iparams->gb.bmlt);
1496 break;
1497 case F_CMAP:
1498 gmx_fio_do_int(fio,iparams->cmap.cmapA);
1499 gmx_fio_do_int(fio,iparams->cmap.cmapB);
1500 break;
1501 default:
1502 gmx_fatal(FARGS,"unknown function type %d (%s) in %s line %d",
1503 ftype,interaction_function[ftype].name,__FILE__,__LINE__);
1504 }
1505 if (!bRead)
1506 gmx_fio_unset_comment(fio);
1507 }
1508
1509 static void do_ilist(t_fileio *fio, t_ilist *ilist,gmx_bool bRead,int file_version,
1510 int ftype)
1511 {
1512 int i,k,idum;
1513 gmx_bool bDum=TRUE;
1514
1515 if (!bRead) {
1516 gmx_fio_set_comment(fio, interaction_function[ftype].name);
1517 }
1518 if (file_version < 44) {
1519 for(i=0; i<MAXNODES; i++)
1520 gmx_fio_do_int(fio,idum);
1521 }
1522 gmx_fio_do_int(fio,ilist->nr);
1523 if (bRead)
1524 snew(ilist->iatoms,ilist->nr);
1525 bDum=gmx_fio_ndo_int(fio,ilist->iatoms,ilist->nr);
1526 if (!bRead)
1527 gmx_fio_unset_comment(fio);
1528 }
1529
1530 static void do_ffparams(t_fileio *fio, gmx_ffparams_t *ffparams,
1531 gmx_bool bRead, int file_version)
1532 {
1533 int idum,i,j;
1534 gmx_bool bDum=TRUE;
1535 unsigned int k;
1536
1537 gmx_fio_do_int(fio,ffparams->atnr);
1538 if (file_version < 57) {
1539 gmx_fio_do_int(fio,idum);
1540 }
1541 gmx_fio_do_int(fio,ffparams->ntypes);
1542 if (bRead && debug)
1543 fprintf(debug,"ffparams->atnr = %d, ntypes = %d\n",
1544 ffparams->atnr,ffparams->ntypes);
1545 if (bRead) {
1546 snew(ffparams->functype,ffparams->ntypes);
1547 snew(ffparams->iparams,ffparams->ntypes);
1548 }
1549 /* Read/write all the function types */
1550 bDum=gmx_fio_ndo_int(fio,ffparams->functype,ffparams->ntypes);
1551 if (bRead && debug)
1552 pr_ivec(debug,0,"functype",ffparams->functype,ffparams->ntypes,TRUE);
1553
1554 if (file_version >= 66) {
1555 gmx_fio_do_double(fio,ffparams->reppow);
1556 } else {
1557 ffparams->reppow = 12.0;
1558 }
1559
1560 if (file_version >= 57) {
1561 gmx_fio_do_real(fio,ffparams->fudgeQQ);
1562 }
1563
1564 /* Check whether all these function types are supported by the code.
1565 * In practice the code is backwards compatible, which means that the
1566 * numbering may have to be altered from old numbering to new numbering
1567 */
1568 for (i=0; (i<ffparams->ntypes); i++) {
1569 if (bRead)
1570 /* Loop over file versions */
1571 for (k=0; (k<NFTUPD); k++)
1572 /* Compare the read file_version to the update table */
1573 if ((file_version < ftupd[k].fvnr) &&
1574 (ffparams->functype[i] >= ftupd[k].ftype)) {
1575 ffparams->functype[i] += 1;
1576 if (debug) {
1577 fprintf(debug,"Incrementing function type %d to %d (due to %s)\n",
1578 i,ffparams->functype[i],
1579 interaction_function[ftupd[k].ftype].longname);
1580 fflush(debug);
1581 }
1582 }
1583
1584 do_iparams(fio, ffparams->functype[i],&ffparams->iparams[i],bRead,
1585 file_version);
1586 if (bRead && debug)
1587 pr_iparams(debug,ffparams->functype[i],&ffparams->iparams[i]);
1588 }
1589 }
1590
1591 static void add_settle_atoms(t_ilist *ilist)
1592 {
1593 int i;
1594
1595 /* Settle used to only store the first atom: add the other two */
1596 srenew(ilist->iatoms,2*ilist->nr);
1597 for(i=ilist->nr/2-1; i>=0; i--)
1598 {
1599 ilist->iatoms[4*i+0] = ilist->iatoms[2*i+0];
1600 ilist->iatoms[4*i+1] = ilist->iatoms[2*i+1];
1601 ilist->iatoms[4*i+2] = ilist->iatoms[2*i+1] + 1;
1602 ilist->iatoms[4*i+3] = ilist->iatoms[2*i+1] + 2;
1603 }
1604 ilist->nr = 2*ilist->nr;
1605 }
1606
1607 static void do_ilists(t_fileio *fio, t_ilist *ilist,gmx_bool bRead,
1608 int file_version)
1609 {
1610 int i,j,renum[F_NRE];
1611 gmx_bool bDum=TRUE,bClear;
1612 unsigned int k;
1613
1614 for(j=0; (j<F_NRE); j++) {
1615 bClear = FALSE;
1616 if (bRead)
1617 for (k=0; k<NFTUPD; k++)
1618 if ((file_version < ftupd[k].fvnr) && (j == ftupd[k].ftype))
1619 bClear = TRUE;
1620 if (bClear) {
1621 ilist[j].nr = 0;
1622 ilist[j].iatoms = NULL;
1623 } else {
1624 do_ilist(fio, &ilist[j],bRead,file_version,j);
1625 if (file_version < 78 && j == F_SETTLE && ilist[j].nr > 0)
1626 {
1627 add_settle_atoms(&ilist[j]);
1628 }
1629 }
1630 /*
1631 if (bRead && gmx_debug_at)
1632 pr_ilist(debug,0,interaction_function[j].longname,
1633 functype,&ilist[j],TRUE);
1634 */
1635 }
1636 }
1637
1638 static void do_idef(t_fileio *fio, gmx_ffparams_t *ffparams,gmx_moltype_t *molt,
1639 gmx_bool bRead, int file_version)
1640 {
1641 do_ffparams(fio, ffparams,bRead,file_version);
1642
1643 if (file_version >= 54) {
1644 gmx_fio_do_real(fio,ffparams->fudgeQQ);
1645 }
1646
1647 do_ilists(fio, molt->ilist,bRead,file_version);
1648 }
1649
1650 static void do_block(t_fileio *fio, t_block *block,gmx_bool bRead,int file_version)
1651 {
1652 int i,idum,dum_nra,*dum_a;
1653 gmx_bool bDum=TRUE;
1654
1655 if (file_version < 44)
1656 for(i=0; i<MAXNODES; i++)
1657 gmx_fio_do_int(fio,idum);
1658 gmx_fio_do_int(fio,block->nr);
1659 if (file_version < 51)
1660 gmx_fio_do_int(fio,dum_nra);
1661 if (bRead) {
1662 block->nalloc_index = block->nr+1;
1663 snew(block->index,block->nalloc_index);
1664 }
1665 bDum=gmx_fio_ndo_int(fio,block->index,block->nr+1);
1666
1667 if (file_version < 51 && dum_nra > 0) {
1668 snew(dum_a,dum_nra);
1669 bDum=gmx_fio_ndo_int(fio,dum_a,dum_nra);
1670 sfree(dum_a);
1671 }
1672 }
1673
1674 static void do_blocka(t_fileio *fio, t_blocka *block,gmx_bool bRead,
1675 int file_version)
1676 {
1677 int i,idum;
1678 gmx_bool bDum=TRUE;
1679
1680 if (file_version < 44)
1681 for(i=0; i<MAXNODES; i++)
1682 gmx_fio_do_int(fio,idum);
1683 gmx_fio_do_int(fio,block->nr);
1684 gmx_fio_do_int(fio,block->nra);
1685 if (bRead) {
1686 block->nalloc_index = block->nr+1;
1687 snew(block->index,block->nalloc_index);
1688 block->nalloc_a = block->nra;
1689 snew(block->a,block->nalloc_a);
1690 }
1691 bDum=gmx_fio_ndo_int(fio,block->index,block->nr+1);
1692 bDum=gmx_fio_ndo_int(fio,block->a,block->nra);
1693 }
1694
1695 static void do_atom(t_fileio *fio, t_atom *atom,int ngrp,gmx_bool bRead,
1696 int file_version, gmx_groups_t *groups,int atnr)
1697 {
1698 int i,myngrp;
1699
1700 gmx_fio_do_real(fio,atom->m);
1701 gmx_fio_do_real(fio,atom->q);
1702 gmx_fio_do_real(fio,atom->mB);
1703 gmx_fio_do_real(fio,atom->qB);
1704 gmx_fio_do_ushort(fio, atom->type);
1705 gmx_fio_do_ushort(fio, atom->typeB);
1706 gmx_fio_do_int(fio,atom->ptype);
1707 gmx_fio_do_int(fio,atom->resind);
1708 if (file_version >= 52)
1709 gmx_fio_do_int(fio,atom->atomnumber);
1710 else if (bRead)
1711 atom->atomnumber = NOTSET;
1712 if (file_version < 23)
1713 myngrp = 8;
1714 else if (file_version < 39)
1715 myngrp = 9;
1716 else
1717 myngrp = ngrp;
1718
1719 if (file_version < 57) {
1720 unsigned char uchar[egcNR];
1721 gmx_fio_ndo_uchar(fio,uchar,myngrp);
1722 for(i=myngrp; (i<ngrp); i++) {
1723 uchar[i] = 0;
1724 }
1725 /* Copy the old data format to the groups struct */
1726 for(i=0; i<ngrp; i++) {
1727 groups->grpnr[i][atnr] = uchar[i];
1728 }
1729 }
1730 }
1731
1732 static void do_grps(t_fileio *fio, int ngrp,t_grps grps[],gmx_bool bRead,
1733 int file_version)
1734 {
1735 int i,j,myngrp;
1736 gmx_bool bDum=TRUE;
1737
1738 if (file_version < 23)
1739 myngrp = 8;
1740 else if (file_version < 39)
1741 myngrp = 9;
1742 else
1743 myngrp = ngrp;
1744
1745 for(j=0; (j<ngrp); j++) {
1746 if (j<myngrp) {
1747 gmx_fio_do_int(fio,grps[j].nr);
1748 if (bRead)
1749 snew(grps[j].nm_ind,grps[j].nr);
1750 bDum=gmx_fio_ndo_int(fio,grps[j].nm_ind,grps[j].nr);
1751 }
1752 else {
1753 grps[j].nr = 1;
1754 snew(grps[j].nm_ind,grps[j].nr);
1755 }
1756 }
1757 }
1758
1759 static void do_symstr(t_fileio *fio, char ***nm,gmx_bool bRead,t_symtab *symtab)
1760 {
1761 int ls;
1762
1763 if (bRead) {
1764 gmx_fio_do_int(fio,ls);
1765 *nm = get_symtab_handle(symtab,ls);
1766 }
1767 else {
1768 ls = lookup_symtab(symtab,*nm);
1769 gmx_fio_do_int(fio,ls);
1770 }
1771 }
1772
1773 static void do_strstr(t_fileio *fio, int nstr,char ***nm,gmx_bool bRead,
1774 t_symtab *symtab)
1775 {
1776 int j;
1777
1778 for (j=0; (j<nstr); j++)
1779 do_symstr(fio, &(nm[j]),bRead,symtab);
1780 }
1781
1782 static void do_resinfo(t_fileio *fio, int n,t_resinfo *ri,gmx_bool bRead,
1783 t_symtab *symtab, int file_version)
1784 {
1785 int j;
1786
1787 for (j=0; (j<n); j++) {
1788 do_symstr(fio, &(ri[j].name),bRead,symtab);
1789 if (file_version >= 63) {
1790 gmx_fio_do_int(fio,ri[j].nr);
1791 gmx_fio_do_uchar(fio, ri[j].ic);
1792 } else {
1793 ri[j].nr = j + 1;
1794 ri[j].ic = ' ';
1795 }
1796 }
1797 }
1798
1799 static void do_atoms(t_fileio *fio, t_atoms *atoms,gmx_bool bRead,t_symtab *symtab,
1800 int file_version,
1801 gmx_groups_t *groups)
1802 {
1803 int i;
1804
1805 gmx_fio_do_int(fio,atoms->nr);
1806 gmx_fio_do_int(fio,atoms->nres);
1807 if (file_version < 57) {
1808 gmx_fio_do_int(fio,groups->ngrpname);
1809 for(i=0; i<egcNR; i++) {
1810 groups->ngrpnr[i] = atoms->nr;
1811 snew(groups->grpnr[i],groups->ngrpnr[i]);
1812 }
1813 }
1814 if (bRead) {
1815 snew(atoms->atom,atoms->nr);
1816 snew(atoms->atomname,atoms->nr);
1817 snew(atoms->atomtype,atoms->nr);
1818 snew(atoms->atomtypeB,atoms->nr);
1819 snew(atoms->resinfo,atoms->nres);
1820 if (file_version < 57) {
1821 snew(groups->grpname,groups->ngrpname);
1822 }
1823 atoms->pdbinfo = NULL;
1824 }
1825 for(i=0; (i<atoms->nr); i++) {
1826 do_atom(fio, &atoms->atom[i],egcNR,bRead, file_version,groups,i);
1827 }
1828 do_strstr(fio, atoms->nr,atoms->atomname,bRead,symtab);
1829 if (bRead && (file_version <= 20)) {
1830 for(i=0; i<atoms->nr; i++) {
1831 atoms->atomtype[i] = put_symtab(symtab,"?");
1832 atoms->atomtypeB[i] = put_symtab(symtab,"?");
1833 }
1834 } else {
1835 do_strstr(fio, atoms->nr,atoms->atomtype,bRead,symtab);
1836 do_strstr(fio, atoms->nr,atoms->atomtypeB,bRead,symtab);
1837 }
1838 do_resinfo(fio, atoms->nres,atoms->resinfo,bRead,symtab,file_version);
1839
1840 if (file_version < 57) {
1841 do_strstr(fio, groups->ngrpname,groups->grpname,bRead,symtab);
1842
1843 do_grps(fio, egcNR,groups->grps,bRead,file_version);
1844 }
1845 }
1846
1847 static void do_groups(t_fileio *fio, gmx_groups_t *groups,
1848 gmx_bool bRead,t_symtab *symtab,
1849 int file_version)
1850 {
1851 int g,n,i;
1852 gmx_bool bDum=TRUE;
1853
1854 do_grps(fio, egcNR,groups->grps,bRead,file_version);
1855 gmx_fio_do_int(fio,groups->ngrpname);
1856 if (bRead) {
1857 snew(groups->grpname,groups->ngrpname);
1858 }
1859 do_strstr(fio, groups->ngrpname,groups->grpname,bRead,symtab);
1860 for(g=0; g<egcNR; g++) {
1861 gmx_fio_do_int(fio,groups->ngrpnr[g]);
1862 if (groups->ngrpnr[g] == 0) {
1863 if (bRead) {
1864 groups->grpnr[g] = NULL;
1865 }
1866 } else {
1867 if (bRead) {
1868 snew(groups->grpnr[g],groups->ngrpnr[g]);
1869 }
1870 bDum=gmx_fio_ndo_uchar(fio, groups->grpnr[g],groups->ngrpnr[g]);
1871 }
1872 }
1873 }
1874
1875 static void do_atomtypes(t_fileio *fio, t_atomtypes *atomtypes,gmx_bool bRead,
1876 t_symtab *symtab,int file_version)
1877 {
1878 int i,j;
1879 gmx_bool bDum = TRUE;
1880
1881 if (file_version > 25) {
1882 gmx_fio_do_int(fio,atomtypes->nr);
1883 j=atomtypes->nr;
1884 if (bRead) {
1885 snew(atomtypes->radius,j);
1886 snew(atomtypes->vol,j);
1887 snew(atomtypes->surftens,j);
1888 snew(atomtypes->atomnumber,j);
1889 snew(atomtypes->gb_radius,j);
1890 snew(atomtypes->S_hct,j);
1891 }
1892 bDum=gmx_fio_ndo_real(fio,atomtypes->radius,j);
1893 bDum=gmx_fio_ndo_real(fio,atomtypes->vol,j);
1894 bDum=gmx_fio_ndo_real(fio,atomtypes->surftens,j);
1895 if(file_version >= 40)
1896 {
1897 bDum=gmx_fio_ndo_int(fio,atomtypes->atomnumber,j);
1898 }
1899 if(file_version >= 60)
1900 {
1901 bDum=gmx_fio_ndo_real(fio,atomtypes->gb_radius,j);
1902 bDum=gmx_fio_ndo_real(fio,atomtypes->S_hct,j);
1903 }
1904 } else {
1905 /* File versions prior to 26 cannot do GBSA,
1906 * so they dont use this structure
1907 */
1908 atomtypes->nr = 0;
1909 atomtypes->radius = NULL;
1910 atomtypes->vol = NULL;
1911 atomtypes->surftens = NULL;
1912 atomtypes->atomnumber = NULL;
1913 atomtypes->gb_radius = NULL;
1914 atomtypes->S_hct = NULL;
1915 }
1916 }
1917
1918 static void do_symtab(t_fileio *fio, t_symtab *symtab,gmx_bool bRead)
1919 {
1920 int i,nr;
1921 t_symbuf *symbuf;
1922 char buf[STRLEN];
1923
1924 gmx_fio_do_int(fio,symtab->nr);
1925 nr = symtab->nr;
1926 if (bRead) {
1927 snew(symtab->symbuf,1);
1928 symbuf = symtab->symbuf;
1929 symbuf->bufsize = nr;
1930 snew(symbuf->buf,nr);
1931 for (i=0; (i<nr); i++) {
1932 gmx_fio_do_string(fio,buf);
1933 symbuf->buf[i]=strdup(buf);
1934 }
1935 }
1936 else {
1937 symbuf = symtab->symbuf;
1938 while (symbuf!=NULL) {
1939 for (i=0; (i<symbuf->bufsize) && (i<nr); i++)
1940 gmx_fio_do_string(fio,symbuf->buf[i]);
1941 nr-=i;
1942 symbuf=symbuf->next;
1943 }
1944 if (nr != 0)
1945 gmx_fatal(FARGS,"nr of symtab strings left: %d",nr);
1946 }
1947 }
1948
1949 static void do_cmap(t_fileio *fio, gmx_cmap_t *cmap_grid, gmx_bool bRead)
1950 {
1951 int i,j,ngrid,gs,nelem;
1952
1953 gmx_fio_do_int(fio,cmap_grid->ngrid);
1954 gmx_fio_do_int(fio,cmap_grid->grid_spacing);
1955
1956 ngrid = cmap_grid->ngrid;
1957 gs = cmap_grid->grid_spacing;
1958 nelem = gs * gs;
1959
1960 if(bRead)
1961 {
1962 snew(cmap_grid->cmapdata,ngrid);
1963
1964 for(i=0;i<cmap_grid->ngrid;i++)
1965 {
1966 snew(cmap_grid->cmapdata[i].cmap,4*nelem);
1967 }
1968 }
1969
1970 for(i=0;i<cmap_grid->ngrid;i++)
1971 {
1972 for(j=0;j<nelem;j++)
1973 {
1974 gmx_fio_do_real(fio,cmap_grid->cmapdata[i].cmap[j*4]);
1975 gmx_fio_do_real(fio,cmap_grid->cmapdata[i].cmap[j*4+1]);
1976 gmx_fio_do_real(fio,cmap_grid->cmapdata[i].cmap[j*4+2]);
1977 gmx_fio_do_real(fio,cmap_grid->cmapdata[i].cmap[j*4+3]);
1978 }
1979 }
1980 }
1981
1982
1983 void tpx_make_chain_identifiers(t_atoms *atoms,t_block *mols)
1984 {
1985 int m,a,a0,a1,r;
1986 char c,chainid;
1987 int chainnum;
1988
1989 /* We always assign a new chain number, but save the chain id characters
1990 * for larger molecules.
1991 */
1992 #define CHAIN_MIN_ATOMS 15
1993
1994 chainnum=0;
1995 chainid='A';
1996 for(m=0; m<mols->nr; m++)
1997 {
1998 a0=mols->index[m];
1999 a1=mols->index[m+1];
2000 if ((a1-a0 >= CHAIN_MIN_ATOMS) && (chainid <= 'Z'))
2001 {
2002 c=chainid;
2003 chainid++;
2004 }
2005 else
2006 {
2007 c=' ';
2008 }
2009 for(a=a0; a<a1; a++)
2010 {
2011 atoms->resinfo[atoms->atom[a].resind].chainnum = chainnum;
2012 atoms->resinfo[atoms->atom[a].resind].chainid = c;
2013 }
2014 chainnum++;
2015 }
2016
2017 /* Blank out the chain id if there was only one chain */
2018 if (chainid == 'B')
2019 {
2020 for(r=0; r<atoms->nres; r++)
2021 {
2022 atoms->resinfo[r].chainid = ' ';
2023 }
2024 }
2025 }
2026
2027 static void do_moltype(t_fileio *fio, gmx_moltype_t *molt,gmx_bool bRead,
2028 t_symtab *symtab, int file_version,
2029 gmx_groups_t *groups)
2030 {
2031 int i;
2032
2033 if (file_version >= 57) {
2034 do_symstr(fio, &(molt->name),bRead,symtab);
2035 }
2036
2037 do_atoms(fio, &molt->atoms, bRead, symtab, file_version, groups);
2038
2039 if (bRead && gmx_debug_at) {
2040 pr_atoms(debug,0,"atoms",&molt->atoms,TRUE);
2041 }
2042
2043 if (file_version >= 57) {
2044 do_ilists(fio, molt->ilist,bRead,file_version);
2045
2046 do_block(fio, &molt->cgs,bRead,file_version);
2047 if (bRead && gmx_debug_at) {
2048 pr_block(debug,0,"cgs",&molt->cgs,TRUE);
2049 }
2050 }
2051
2052 /* This used to be in the atoms struct */
2053 do_blocka(fio, &molt->excls, bRead, file_version);
2054 }
2055
2056 static void do_molblock(t_fileio *fio, gmx_molblock_t *molb,gmx_bool bRead,
2057 int file_version)
2058 {
2059 int i;
2060
2061 gmx_fio_do_int(fio,molb->type);
2062 gmx_fio_do_int(fio,molb->nmol);
2063 gmx_fio_do_int(fio,molb->natoms_mol);
2064 /* Position restraint coordinates */
2065 gmx_fio_do_int(fio,molb->nposres_xA);
2066 if (molb->nposres_xA > 0) {
2067 if (bRead) {
2068 snew(molb->posres_xA,molb->nposres_xA);
2069 }
2070 gmx_fio_ndo_rvec(fio,molb->posres_xA,molb->nposres_xA);
2071 }
2072 gmx_fio_do_int(fio,molb->nposres_xB);
2073 if (molb->nposres_xB > 0) {
2074 if (bRead) {
2075 snew(molb->posres_xB,molb->nposres_xB);
2076 }
2077 gmx_fio_ndo_rvec(fio,molb->posres_xB,molb->nposres_xB);
2078 }
2079
2080 }
2081
2082 static t_block mtop_mols(gmx_mtop_t *mtop)
2083 {
2084 int mb,m,a,mol;
2085 t_block mols;
2086
2087 mols.nr = 0;
2088 for(mb=0; mb<mtop->nmolblock; mb++) {
2089 mols.nr += mtop->molblock[mb].nmol;
2090 }
2091 mols.nalloc_index = mols.nr + 1;
2092 snew(mols.index,mols.nalloc_index);
2093
2094 a = 0;
2095 m = 0;
2096 mols.index[m] = a;
2097 for(mb=0; mb<mtop->nmolblock; mb++) {
2098 for(mol=0; mol<mtop->molblock[mb].nmol; mol++) {
2099 a += mtop->molblock[mb].natoms_mol;
2100 m++;
2101 mols.index[m] = a;
2102 }
2103 }
2104
2105 return mols;
2106 }
2107
2108 static void add_posres_molblock(gmx_mtop_t *mtop)
2109 {
2110 t_ilist *il;
2111 int am,i,mol,a;
2112 gmx_bool bFE;
2113 gmx_molblock_t *molb;
2114 t_iparams *ip;
2115
2116 il = &mtop->moltype[0].ilist[F_POSRES];
2117 if (il->nr == 0) {
2118 return;
2119 }
2120 am = 0;
2121 bFE = FALSE;
2122 for(i=0; i<il->nr; i+=2) {
2123 ip = &mtop->ffparams.iparams[il->iatoms[i]];
2124 am = max(am,il->iatoms[i+1]);
2125 if (ip->posres.pos0B[XX] != ip->posres.pos0A[XX] ||
2126 ip->posres.pos0B[YY] != ip->posres.pos0A[YY] ||
2127 ip->posres.pos0B[ZZ] != ip->posres.pos0A[ZZ]) {
2128 bFE = TRUE;
2129 }
2130 }
2131 /* Make the posres coordinate block end at a molecule end */
2132 mol = 0;
2133 while(am >= mtop->mols.index[mol+1]) {
2134 mol++;
2135 }
2136 molb = &mtop->molblock[0];
2137 molb->nposres_xA = mtop->mols.index[mol+1];
2138 snew(molb->posres_xA,molb->nposres_xA);
2139 if (bFE) {
2140 molb->nposres_xB = molb->nposres_xA;
2141 snew(molb->posres_xB,molb->nposres_xB);
2142 } else {
2143 molb->nposres_xB = 0;
2144 }
2145 for(i=0; i<il->nr; i+=2) {
2146 ip = &mtop->ffparams.iparams[il->iatoms[i]];
2147 a = il->iatoms[i+1];
2148 molb->posres_xA[a][XX] = ip->posres.pos0A[XX];
2149 molb->posres_xA[a][YY] = ip->posres.pos0A[YY];
2150 molb->posres_xA[a][ZZ] = ip->posres.pos0A[ZZ];
2151 if (bFE) {
2152 molb->posres_xB[a][XX] = ip->posres.pos0B[XX];
2153 molb->posres_xB[a][YY] = ip->posres.pos0B[YY];
2154 molb->posres_xB[a][ZZ] = ip->posres.pos0B[ZZ];
2155 }
2156 }
2157 }
2158
2159 static void set_disres_npair(gmx_mtop_t *mtop)
2160 {
2161 int mt,i,npair;
2162 t_iparams *ip;
2163 t_ilist *il;
2164 t_iatom *a;
2165
2166 ip = mtop->ffparams.iparams;
2167
2168 for(mt=0; mt<mtop->nmoltype; mt++) {
2169 il = &mtop->moltype[mt].ilist[F_DISRES];
2170 if (il->nr > 0) {
2171 a = il->iatoms;
2172 npair = 0;
2173 for(i=0; i<il->nr; i+=3) {
2174 npair++;
2175 if (i+3 == il->nr || ip[a[i]].disres.label != ip[a[i+3]].disres.label) {
2176 ip[a[i]].disres.npair = npair;
2177 npair = 0;
2178 }
2179 }
2180 }
2181 }
2182 }
2183
2184 static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop,gmx_bool bRead,
2185 int file_version)
2186 {
2187 int mt,mb,i;
2188 t_blocka dumb;
2189
2190 if (bRead)
2191 init_mtop(mtop);
2192 do_symtab(fio, &(mtop->symtab),bRead);
2193 if (bRead && debug)
2194 pr_symtab(debug,0,"symtab",&mtop->symtab);
2195
2196 do_symstr(fio, &(mtop->name),bRead,&(mtop->symtab));
2197
2198 if (file_version >= 57) {
2199 do_ffparams(fio, &mtop->ffparams,bRead,file_version);
2200
2201 gmx_fio_do_int(fio,mtop->nmoltype);
2202 } else {
2203 mtop->nmoltype = 1;
2204 }
2205 if (bRead) {
2206 snew(mtop->moltype,mtop->nmoltype);
2207 if (file_version < 57) {
2208 mtop->moltype[0].name = mtop->name;
2209 }
2210 }
2211 for(mt=0; mt<mtop->nmoltype; mt++) {
2212 do_moltype(fio, &mtop->moltype[mt],bRead,&mtop->symtab,file_version,
2213 &mtop->groups);
2214 }
2215
2216 if (file_version >= 57) {
2217 gmx_fio_do_int(fio,mtop->nmolblock);
2218 } else {
2219 mtop->nmolblock = 1;
2220 }
2221 if (bRead) {
2222 snew(mtop->molblock,mtop->nmolblock);
2223 }
2224 if (file_version >= 57) {
2225 for(mb=0; mb<mtop->nmolblock; mb++) {
2226 do_molblock(fio, &mtop->molblock[mb],bRead,file_version);
2227 }
2228 gmx_fio_do_int(fio,mtop->natoms);
2229 } else {
2230 mtop->molblock[0].type = 0;
2231 mtop->molblock[0].nmol = 1;
2232 mtop->molblock[0].natoms_mol = mtop->moltype[0].atoms.nr;
2233 mtop->molblock[0].nposres_xA = 0;
2234 mtop->molblock[0].nposres_xB = 0;
2235 }
2236
2237 do_atomtypes (fio, &(mtop->atomtypes),bRead,&(mtop->symtab), file_version);
2238 if (bRead && debug)
2239 pr_atomtypes(debug,0,"atomtypes",&mtop->atomtypes,TRUE);
2240
2241 if (file_version < 57) {
2242 /* Debug statements are inside do_idef */
2243 do_idef (fio, &mtop->ffparams,&mtop->moltype[0],bRead,file_version);
2244 mtop->natoms = mtop->moltype[0].atoms.nr;
2245 }
2246
2247 if(file_version >= 65)
2248 {
2249 do_cmap(fio, &mtop->ffparams.cmap_grid,bRead);
2250 }
2251 else
2252 {
2253 mtop->ffparams.cmap_grid.ngrid = 0;
2254 mtop->ffparams.cmap_grid.grid_spacing = 0;
2255 mtop->ffparams.cmap_grid.cmapdata = NULL;
2256 }
2257
2258 if (file_version >= 57) {
2259 do_groups(fio, &mtop->groups,bRead,&(mtop->symtab),file_version);
2260 }
2261
2262 if (file_version < 57) {
2263 do_block(fio, &mtop->moltype[0].cgs,bRead,file_version);
2264 if (bRead && gmx_debug_at) {
2265 pr_block(debug,0,"cgs",&mtop->moltype[0].cgs,TRUE);
2266 }
2267 do_block(fio, &mtop->mols,bRead,file_version);
2268 /* Add the posres coordinates to the molblock */
2269 add_posres_molblock(mtop);
2270 }
2271 if (bRead) {
2272 if (file_version >= 57) {
2273 mtop->mols = mtop_mols(mtop);
2274 }
2275 if (gmx_debug_at) {
2276 pr_block(debug,0,"mols",&mtop->mols,TRUE);
2277 }
2278 }
2279
2280 if (file_version < 51) {
2281 /* Here used to be the shake blocks */
2282 do_blocka(fio, &dumb,bRead,file_version);
2283 if (dumb.nr > 0)
2284 sfree(dumb.index);
2285 if (dumb.nra > 0)
2286 sfree(dumb.a);
2287 }
2288
2289 if (bRead) {
2290 close_symtab(&(mtop->symtab));
2291 }
2292 }
2293
2294 /* If TopOnlyOK is TRUE then we can read even future versions
2295 * of tpx files, provided the file_generation hasn't changed.
2296 * If it is FALSE, we need the inputrecord too, and bail out
2297 * if the file is newer than the program.
2298 *
2299 * The version and generation if the topology (see top of this file)
2300 * are returned in the two last arguments.
2301 *
2302 * If possible, we will read the inputrec even when TopOnlyOK is TRUE.
2303 */
2304 static void do_tpxheader(t_fileio *fio,gmx_bool bRead,t_tpxheader *tpx,
2305 gmx_bool TopOnlyOK, int *file_version,
2306 int *file_generation)
2307 {
2308 char buf[STRLEN];
2309 char file_tag[STRLEN];
2310 gmx_bool bDouble;
2311 int precision;
2312 int fver,fgen;
2313 int idum=0;
2314 real rdum=0;
2315
2316 gmx_fio_checktype(fio);
2317 gmx_fio_setdebug(fio,bDebugMode());
2318
2319 /* NEW! XDR tpb file */
2320 precision = sizeof(real);
2321 if (bRead) {
2322 gmx_fio_do_string(fio,buf);
2323 if (strncmp(buf,"VERSION",7))
2324 gmx_fatal(FARGS,"Can not read file %s,\n"
2325 " this file is from a Gromacs version which is older than 2.0\n"
2326 " Make a new one with grompp or use a gro or pdb file, if possible",
2327 gmx_fio_getname(fio));
2328 gmx_fio_do_int(fio,precision);
2329 bDouble = (precision == sizeof(double));
2330 if ((precision != sizeof(float)) && !bDouble)
2331 gmx_fatal(FARGS,"Unknown precision in file %s: real is %d bytes "
2332 "instead of %d or %d",
2333 gmx_fio_getname(fio),precision,sizeof(float),sizeof(double));
2334 gmx_fio_setprecision(fio,bDouble);
2335 fprintf(stderr,"Reading file %s, %s (%s precision)\n",
2336 gmx_fio_getname(fio),buf,bDouble ? "double" : "single");
2337 }
2338 else {
2339 gmx_fio_write_string(fio,GromacsVersion());
2340 bDouble = (precision == sizeof(double));
2341 gmx_fio_setprecision(fio,bDouble);
2342 gmx_fio_do_int(fio,precision);
2343 fver = tpx_version;
2344 sprintf(file_tag,"%s",tpx_tag);
2345 fgen = tpx_generation;
2346 }
2347
2348 /* Check versions! */
2349 gmx_fio_do_int(fio,fver);
2350
2351 if (fver >= 77)
2352 {
2353 gmx_fio_do_string(fio,file_tag);
2354 }
2355 if (bRead)
2356 {
2357 if (fver < 77)
2358 {
2359 /* Versions before 77 don't have the tag, set it to release */
2360 sprintf(file_tag,"%s",TPX_TAG_RELEASE);
2361 }
2362
2363 if (strcmp(file_tag,tpx_tag) != 0)
2364 {
2365 fprintf(stderr,"Note: file tpx tag '%s', software tpx tag '%s'\n",
2366 file_tag,tpx_tag);
2367
2368 /* We only support reading tpx files with the same tag as the code
2369 * or tpx files with the release tag and with lower version number.
2370 */
2371 if (!(strcmp(file_tag,TPX_TAG_RELEASE) == 0 && fver < tpx_version))
2372 {
2373 gmx_fatal(FARGS,"tpx tag/version mismatch: reading tpx file (%s) version %d, tag '%s' with program for tpx version %d, tag '%s'",
2374 gmx_fio_getname(fio),fver,file_tag,
2375 tpx_version,tpx_tag);
2376 }
2377 }
2378 }
2379
2380 if (fver >= 26)
2381 {
2382 gmx_fio_do_int(fio,fgen);
2383 }
2384 else
2385 {
2386 fgen=0;
2387 }
2388
2389 if (file_version != NULL)
2390 {
2391 *file_version = fver;
2392 }
2393 if (file_generation != NULL)
2394 {
2395 *file_generation = fgen;
2396 }
2397
2398
2399 if ((fver <= tpx_incompatible_version) ||
2400 ((fver > tpx_version) && !TopOnlyOK) ||
2401 (fgen > tpx_generation))
2402 gmx_fatal(FARGS,"reading tpx file (%s) version %d with version %d program",
2403 gmx_fio_getname(fio),fver,tpx_version);
2404
2405 do_section(fio,eitemHEADER,bRead);
2406 gmx_fio_do_int(fio,tpx->natoms);
2407 if (fver >= 28)
2408 gmx_fio_do_int(fio,tpx->ngtc);
2409 else
2410 tpx->ngtc = 0;
2411 if (fver < 62) {
2412 gmx_fio_do_int(fio,idum);
2413 gmx_fio_do_real(fio,rdum);
2414 }
2415 /*a better decision will eventually (5.0 or later) need to be made
2416 on how to treat the alchemical state of the system, which can now
2417 vary through a simulation, and cannot be completely described
2418 though a single lambda variable, or even a single state
2419 index. Eventually, should probably be a vector. MRS*/
2420 if (fver >= 79)
2421 {
2422 gmx_fio_do_int(fio,tpx->fep_state);
2423 }
2424 gmx_fio_do_real(fio,tpx->lambda);
2425 gmx_fio_do_int(fio,tpx->bIr);
2426 gmx_fio_do_int(fio,tpx->bTop);
2427 gmx_fio_do_int(fio,tpx->bX);
2428 gmx_fio_do_int(fio,tpx->bV);
2429 gmx_fio_do_int(fio,tpx->bF);
2430 gmx_fio_do_int(fio,tpx->bBox);
2431
2432 if((fgen > tpx_generation)) {
2433 /* This can only happen if TopOnlyOK=TRUE */
2434 tpx->bIr=FALSE;
2435 }
2436 }
2437
2438 static int do_tpx(t_fileio *fio, gmx_bool bRead,
2439 t_inputrec *ir,t_state *state,rvec *f,gmx_mtop_t *mtop,
2440 gmx_bool bXVallocated)
2441 {
2442 t_tpxheader tpx;
2443 t_inputrec dum_ir;
2444 gmx_mtop_t dum_top;
2445 gmx_bool TopOnlyOK,bDum=TRUE;
2446 int file_version,file_generation;
2447 int i;
2448 rvec *xptr,*vptr;
2449 int ePBC;
2450 gmx_bool bPeriodicMols;
2451
2452 if (!bRead) {
2453 tpx.natoms = state->natoms;
2454 tpx.ngtc = state->ngtc; /* need to add nnhpres here? */
2455 tpx.fep_state = state->fep_state;
2456 tpx.lambda = state->lambda[efptFEP];
2457 tpx.bIr = (ir != NULL);
2458 tpx.bTop = (mtop != NULL);
2459 tpx.bX = (state->x != NULL);
2460 tpx.bV = (state->v != NULL);
2461 tpx.bF = (f != NULL);
2462 tpx.bBox = TRUE;
2463 }
2464
2465 TopOnlyOK = (ir==NULL);
2466
2467 do_tpxheader(fio,bRead,&tpx,TopOnlyOK,&file_version,&file_generation);
2468
2469 if (bRead) {
2470 state->flags = 0;
2471 /* state->lambda = tpx.lambda;*/ /*remove this eventually? */
2472 /* The init_state calls initialize the Nose-Hoover xi integrals to zero */
2473 if (bXVallocated) {
2474 xptr = state->x;
2475 vptr = state->v;
2476 init_state(state,0,tpx.ngtc,0,0,0); /* nose-hoover chains */ /* eventually, need to add nnhpres here? */
2477 state->natoms = tpx.natoms;
2478 state->nalloc = tpx.natoms;
2479 state->x = xptr;
2480 state->v = vptr;
2481 } else {
2482 init_state(state,tpx.natoms,tpx.ngtc,0,0,0); /* nose-hoover chains */
2483 }
2484 }
2485
2486 #define do_test(fio,b,p) if (bRead && (p!=NULL) && !b) gmx_fatal(FARGS,"No %s in %s",#p,gmx_fio_getname(fio))
2487
2488 do_test(fio,tpx.bBox,state->box);
2489 do_section(fio,eitemBOX,bRead);
2490 if (tpx.bBox) {
2491 gmx_fio_ndo_rvec(fio,state->box,DIM);
2492 if (file_version >= 51) {
2493 gmx_fio_ndo_rvec(fio,state->box_rel,DIM);
2494 } else {
2495 /* We initialize box_rel after reading the inputrec */
2496 clear_mat(state->box_rel);
2497 }
2498 if (file_version >= 28) {
2499 gmx_fio_ndo_rvec(fio,state->boxv,DIM);
2500 if (file_version < 56) {
2501 matrix mdum;
2502 gmx_fio_ndo_rvec(fio,mdum,DIM);
2503 }
2504 }
2505 }
2506
2507 if (state->ngtc > 0 && file_version >= 28) {
2508 real *dumv;
2509 /*ndo_double(state->nosehoover_xi,state->ngtc,bDum);*/
2510 /*ndo_double(state->nosehoover_vxi,state->ngtc,bDum);*/
2511 /*ndo_double(state->therm_integral,state->ngtc,bDum);*/
2512 snew(dumv,state->ngtc);
2513 if (file_version < 69) {
2514 bDum=gmx_fio_ndo_real(fio,dumv,state->ngtc);
2515 }
2516 /* These used to be the Berendsen tcoupl_lambda's */
2517 bDum=gmx_fio_ndo_real(fio,dumv,state->ngtc);
2518 sfree(dumv);
2519 }
2520
2521 /* Prior to tpx version 26, the inputrec was here.
2522 * I moved it to enable partial forward-compatibility
2523 * for analysis/viewer programs.
2524 */
2525 if(file_version<26) {
2526 do_test(fio,tpx.bIr,ir);
2527 do_section(fio,eitemIR,bRead);
2528 if (tpx.bIr) {
2529 if (ir) {
2530 do_inputrec(fio, ir,bRead,file_version,
2531 mtop ? &mtop->ffparams.fudgeQQ : NULL);
2532 if (bRead && debug)
2533 pr_inputrec(debug,0,"inputrec",ir,FALSE);
2534 }
2535 else {
2536 do_inputrec(fio, &dum_ir,bRead,file_version,
2537 mtop ? &mtop->ffparams.fudgeQQ :NULL);
2538 if (bRead && debug)
2539 pr_inputrec(debug,0,"inputrec",&dum_ir,FALSE);
2540 done_inputrec(&dum_ir);
2541 }
2542
2543 }
2544 }
2545
2546 do_test(fio,tpx.bTop,mtop);
2547 do_section(fio,eitemTOP,bRead);
2548 if (tpx.bTop) {
2549 if (mtop) {
2550 do_mtop(fio,mtop,bRead, file_version);
2551 } else {
2552 do_mtop(fio,&dum_top,bRead,file_version);
2553 done_mtop(&dum_top,TRUE);
2554 }
2555 }
2556 do_test(fio,tpx.bX,state->x);
2557 do_section(fio,eitemX,bRead);
2558 if (tpx.bX) {
2559 if (bRead) {
2560 state->flags |= (1<<estX);
2561 }
2562 gmx_fio_ndo_rvec(fio,state->x,state->natoms);
2563 }
2564
2565 do_test(fio,tpx.bV,state->v);
2566 do_section(fio,eitemV,bRead);
2567 if (tpx.bV) {
2568 if (bRead) {
2569 state->flags |= (1<<estV);
2570 }
2571 gmx_fio_ndo_rvec(fio,state->v,state->natoms);
2572 }
2573
2574 do_test(fio,tpx.bF,f);
2575 do_section(fio,eitemF,bRead);
2576 if (tpx.bF) gmx_fio_ndo_rvec(fio,f,state->natoms);
2577
2578 /* Starting with tpx version 26, we have the inputrec
2579 * at the end of the file, so we can ignore it
2580 * if the file is never than the software (but still the
2581 * same generation - see comments at the top of this file.
2582 *
2583 *
2584 */
2585 ePBC = -1;
2586 bPeriodicMols = FALSE;
2587 if (file_version >= 26) {
2588 do_test(fio,tpx.bIr,ir);
2589 do_section(fio,eitemIR,bRead);
2590 if (tpx.bIr) {
2591 if (file_version >= 53) {
2592 /* Removed the pbc info from do_inputrec, since we always want it */
2593 if (!bRead) {
2594 ePBC = ir->ePBC;
2595 bPeriodicMols = ir->bPeriodicMols;
2596 }
2597 gmx_fio_do_int(fio,ePBC);
2598 gmx_fio_do_gmx_bool(fio,bPeriodicMols);
2599 }
2600 if (file_generation <= tpx_generation && ir) {
2601 do_inputrec(fio, ir,bRead,file_version,mtop ? &mtop->ffparams.fudgeQQ : NULL);
2602 if (bRead && debug)
2603 pr_inputrec(debug,0,"inputrec",ir,FALSE);
2604 if (file_version < 51)
2605 set_box_rel(ir,state);
2606 if (file_version < 53) {
2607 ePBC = ir->ePBC;
2608 bPeriodicMols = ir->bPeriodicMols;
2609 }
2610 }
2611 if (bRead && ir && file_version >= 53) {
2612 /* We need to do this after do_inputrec, since that initializes ir */
2613 ir->ePBC = ePBC;
2614 ir->bPeriodicMols = bPeriodicMols;
2615 }
2616 }
2617 }
2618
2619 if (bRead)
2620 {
2621 if (tpx.bIr && ir)
2622 {
2623 if (state->ngtc == 0)
2624 {
2625 /* Reading old version without tcoupl state data: set it */
2626 init_gtc_state(state,ir->opts.ngtc,0,ir->opts.nhchainlength);
2627 }
2628 if (tpx.bTop && mtop)
2629 {
2630 if (file_version < 57)
2631 {
2632 if (mtop->moltype[0].ilist[F_DISRES].nr > 0)
2633 {
2634 ir->eDisre = edrSimple;
2635 }
2636 else
2637 {
2638 ir->eDisre = edrNone;
2639 }
2640 }
2641 set_disres_npair(mtop);
2642 }
2643 }
2644
2645 if (tpx.bTop && mtop)
2646 {
2647 gmx_mtop_finalize(mtop);
2648 }
2649
2650 if (file_version >= 57)
2651 {
2652 char *env;
2653 int ienv;
2654 env = getenv("GMX_NOCHARGEGROUPS");
2655 if (env != NULL)
2656 {
2657 sscanf(env,"%d",&ienv);
2658 fprintf(stderr,"\nFound env.var. GMX_NOCHARGEGROUPS = %d\n",
2659 ienv);
2660 if (ienv > 0)
2661 {
2662 fprintf(stderr,
2663 "Will make single atomic charge groups in non-solvent%s\n",
2664 ienv > 1 ? " and solvent" : "");
2665 gmx_mtop_make_atomic_charge_groups(mtop,ienv==1);
2666 }
2667 fprintf(stderr,"\n");
2668 }
2669 }
2670 }
2671
2672 return ePBC;
2673 }
2674
2675 /************************************************************
2676 *
2677 * The following routines are the exported ones
2678 *
2679 ************************************************************/
2680
2681 t_fileio *open_tpx(const char *fn,const char *mode)
2682 {
2683 return gmx_fio_open(fn,mode);
2684 }
2685
2686 void close_tpx(t_fileio *fio)
2687 {
2688 gmx_fio_close(fio);
2689 }
2690
2691 void read_tpxheader(const char *fn, t_tpxheader *tpx, gmx_bool TopOnlyOK,
2692 int *file_version, int *file_generation)
2693 {
2694 t_fileio *fio;
2695
2696 fio = open_tpx(fn,"r");
2697 do_tpxheader(fio,TRUE,tpx,TopOnlyOK,file_version,file_generation);
2698 close_tpx(fio);
2699 }
2700
2701 void write_tpx_state(const char *fn,
2702 t_inputrec *ir,t_state *state,gmx_mtop_t *mtop)
2703 {
2704 t_fileio *fio;
2705
2706 fio = open_tpx(fn,"w");
2707 do_tpx(fio,FALSE,ir,state,NULL,mtop,FALSE);
2708 close_tpx(fio);
2709 }
2710
2711 void read_tpx_state(const char *fn,
2712 t_inputrec *ir,t_state *state,rvec *f,gmx_mtop_t *mtop)
2713 {
2714 t_fileio *fio;
2715
2716 fio = open_tpx(fn,"r");
2717 do_tpx(fio,TRUE,ir,state,f,mtop,FALSE);
2718 close_tpx(fio);
2719 }
2720
2721 int read_tpx(const char *fn,
2722 t_inputrec *ir, matrix box,int *natoms,
2723 rvec *x,rvec *v,rvec *f,gmx_mtop_t *mtop)
2724 {
2725 t_fileio *fio;
2726 t_state state;
2727 int ePBC;
2728
2729 state.x = x;
2730 state.v = v;
2731 fio = open_tpx(fn,"r");
2732 ePBC = do_tpx(fio,TRUE,ir,&state,f,mtop,TRUE);
2733 close_tpx(fio);
2734 *natoms = state.natoms;
2735 if (box)
2736 copy_mat(state.box,box);
2737 state.x = NULL;
2738 state.v = NULL;
2739 done_state(&state);
2740
2741 return ePBC;
2742 }
2743
2744 int read_tpx_top(const char *fn,
2745 t_inputrec *ir, matrix box,int *natoms,
2746 rvec *x,rvec *v,rvec *f,t_topology *top)
2747 {
2748 gmx_mtop_t mtop;
2749 t_topology *ltop;
2750 int ePBC;
2751
2752 ePBC = read_tpx(fn,ir,box,natoms,x,v,f,&mtop);
2753
2754 *top = gmx_mtop_t_to_t_topology(&mtop);
2755
2756 return ePBC;
2757 }
2758
2759 gmx_bool fn2bTPX(const char *file)
2760 {
2761 switch (fn2ftp(file)) {
2762 case efTPR:
2763 case efTPB:
2764 case efTPA:
2765 return TRUE;
2766 default:
2767 return FALSE;
2768 }
2769 }
2770
2771 gmx_bool read_tps_conf(const char *infile,char *title,t_topology *top,int *ePBC,
2772 rvec **x,rvec **v,matrix box,gmx_bool bMass)
2773 {
2774 t_tpxheader header;
2775 int natoms,i,version,generation;
2776 gmx_bool bTop,bXNULL=FALSE;
2777 gmx_mtop_t *mtop;
2778 t_topology *topconv;
2779 gmx_atomprop_t aps;
2780
2781 bTop = fn2bTPX(infile);
2782 *ePBC = -1;
2783 if (bTop) {
2784 read_tpxheader(infile,&header,TRUE,&version,&generation);
2785 if (x)
2786 snew(*x,header.natoms);
2787 if (v)
2788 snew(*v,header.natoms);
2789 snew(mtop,1);
2790 *ePBC = read_tpx(infile,NULL,box,&natoms,
2791 (x==NULL) ? NULL : *x,(v==NULL) ? NULL : *v,NULL,mtop);
2792 *top = gmx_mtop_t_to_t_topology(mtop);
2793 sfree(mtop);
2794 strcpy(title,*top->name);
2795 tpx_make_chain_identifiers(&top->atoms,&top->mols);
2796 }
2797 else {
2798 get_stx_coordnum(infile,&natoms);
2799 init_t_atoms(&top->atoms,natoms,(fn2ftp(infile) == efPDB));
2800 if (x == NULL)
2801 {
2802 snew(x,1);
2803 bXNULL = TRUE;
2804 }
2805 snew(*x,natoms);
2806 if (v)
2807 snew(*v,natoms);
2808 read_stx_conf(infile,title,&top->atoms,*x,(v==NULL) ? NULL : *v,ePBC,box);
2809 if (bXNULL)
2810 {
2811 sfree(*x);
2812 sfree(x);
2813 }
2814 if (bMass) {
2815 aps = gmx_atomprop_init();
2816 for(i=0; (i<natoms); i++)
2817 if (!gmx_atomprop_query(aps,epropMass,
2818 *top->atoms.resinfo[top->atoms.atom[i].resind].name,
2819 *top->atoms.atomname[i],
2820 &(top->atoms.atom[i].m))) {
2821 if (debug)
2822 fprintf(debug,"Can not find mass for atom %s %d %s, setting to 1\n",
2823 *top->atoms.resinfo[top->atoms.atom[i].resind].name,
2824 top->atoms.resinfo[top->atoms.atom[i].resind].nr,
2825 *top->atoms.atomname[i]);
2826 }
2827 gmx_atomprop_destroy(aps);
2828 }
2829 top->idef.ntypes=-1;
2830 }
2831
2832 return bTop;
2833 }
The ultimate molecular dynamics simulation package