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