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