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