927f82fa384661a8eb3456882a3575435cbc4790
[alexxy/gromacs.git] / src / gmxlib / checkpoint.c
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2  *
3  * 
4  * This file is part of Gromacs        Copyright (c) 1991-2008
5  * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
6  *
7  * This program is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License
9  * as published by the Free Software Foundation; either version 2
10  * of the License, or (at your option) any later version.
11  *
12  * To help us fund GROMACS development, we humbly ask that you cite
13  * the research papers on the package. Check out http://www.gromacs.org
14  * 
15  * And Hey:
16  * Gnomes, ROck Monsters And Chili Sauce
17  */
18
19 /* The source code in this file should be thread-safe. 
20  Please keep it that way. */
21
22
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 #include "gmx_header_config.h"
27
28 #include <string.h>
29 #include <time.h>
30
31 #ifdef HAVE_SYS_TIME_H
32 #include <sys/time.h>
33 #endif
34
35
36 #ifdef GMX_NATIVE_WINDOWS
37 /* _chsize_s */
38 #include <io.h>
39 #include <sys/locking.h>
40 #endif
41
42
43 #include "filenm.h"
44 #include "names.h"
45 #include "typedefs.h"
46 #include "smalloc.h"
47 #include "gmxfio.h"
48 #include "xdrf.h"
49 #include "statutil.h"
50 #include "txtdump.h"
51 #include "vec.h"
52 #include "network.h"
53 #include "gmx_random.h"
54 #include "checkpoint.h"
55 #include "futil.h"
56 #include "string2.h"
57 #include <fcntl.h>
58
59
60 #ifdef GMX_FAHCORE
61 #include "corewrap.h"
62 #endif
63
64
65 /* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
66 char *
67 gmx_ctime_r(const time_t *clock,char *buf, int n);
68
69
70 #define CPT_MAGIC1 171817
71 #define CPT_MAGIC2 171819
72 #define CPTSTRLEN 1024
73
74 #ifdef GMX_DOUBLE
75 #define GMX_CPT_BUILD_DP 1
76 #else
77 #define GMX_CPT_BUILD_DP 0
78 #endif
79
80 /* cpt_version should normally only be changed
81  * when the header of footer format changes.
82  * The state data format itself is backward and forward compatible.
83  * But old code can not read a new entry that is present in the file
84  * (but can read a new format when new entries are not present).
85  */
86 static const int cpt_version = 14;
87
88
89 const char *est_names[estNR]=
90 {
91     "FE-lambda",
92     "box", "box-rel", "box-v", "pres_prev",
93     "nosehoover-xi", "thermostat-integral",
94     "x", "v", "SDx", "CGp", "LD-rng", "LD-rng-i",
95     "disre_initf", "disre_rm3tav",
96     "orire_initf", "orire_Dtav",
97     "svir_prev", "nosehoover-vxi", "v_eta", "vol0", "nhpres_xi", "nhpres_vxi", "fvir_prev","fep_state", "MC-rng", "MC-rng-i"
98 };
99
100 enum { eeksEKIN_N, eeksEKINH, eeksDEKINDL, eeksMVCOS, eeksEKINF, eeksEKINO, eeksEKINSCALEF, eeksEKINSCALEH, eeksVSCALE, eeksEKINTOTAL, eeksNR };
101
102 const char *eeks_names[eeksNR]=
103 {
104     "Ekin_n", "Ekinh", "dEkindlambda", "mv_cos",
105     "Ekinf", "Ekinh_old", "EkinScaleF_NHC", "EkinScaleH_NHC","Vscale_NHC","Ekin_Total"
106 };
107
108 enum { eenhENERGY_N, eenhENERGY_AVER, eenhENERGY_SUM, eenhENERGY_NSUM,
109        eenhENERGY_SUM_SIM, eenhENERGY_NSUM_SIM,
110        eenhENERGY_NSTEPS, eenhENERGY_NSTEPS_SIM, 
111        eenhENERGY_DELTA_H_NN,
112        eenhENERGY_DELTA_H_LIST, 
113        eenhENERGY_DELTA_H_STARTTIME, 
114        eenhENERGY_DELTA_H_STARTLAMBDA, 
115        eenhNR };
116
117 const char *eenh_names[eenhNR]=
118 {
119     "energy_n", "energy_aver", "energy_sum", "energy_nsum",
120     "energy_sum_sim", "energy_nsum_sim",
121     "energy_nsteps", "energy_nsteps_sim", 
122     "energy_delta_h_nn",
123     "energy_delta_h_list", 
124     "energy_delta_h_start_time", 
125     "energy_delta_h_start_lambda"
126 };
127
128 /* free energy history variables -- need to be preserved over checkpoint */
129 enum { edfhBEQUIL,edfhNATLAMBDA,edfhWLHISTO,edfhWLDELTA,edfhSUMWEIGHTS,edfhSUMDG,edfhSUMMINVAR,edfhSUMVAR,
130        edfhACCUMP,edfhACCUMM,edfhACCUMP2,edfhACCUMM2,edfhTIJ,edfhTIJEMP,edfhNR };
131 /* free energy history variable names  */
132 const char *edfh_names[edfhNR]=
133 {
134     "bEquilibrated","N_at_state", "Wang-Landau_Histogram", "Wang-Landau-delta", "Weights", "Free Energies", "minvar","variance",
135     "accumulated_plus", "accumulated_minus", "accumulated_plus_2",  "accumulated_minus_2", "Tij", "Tij_empirical"
136 };
137
138 #ifdef GMX_NATIVE_WINDOWS
139 static int
140 gmx_wintruncate(const char *filename, __int64 size)
141 {
142 #ifdef GMX_FAHCORE
143     /*we do this elsewhere*/
144     return 0;
145 #else
146     FILE *fp;
147     int   rc;
148     
149     fp=fopen(filename,"rb+");
150     
151     if(fp==NULL)
152     {
153         return -1;
154     }
155     
156     return _chsize_s( fileno(fp), size);
157 #endif
158 }
159 #endif
160
161
162 enum { ecprREAL, ecprRVEC, ecprMATRIX };
163
164 enum { cptpEST, cptpEEKS, cptpEENH, cptpEDFH };
165 /* enums for the different components of checkpoint variables, replacing the hard coded ones.
166    cptpEST - state variables.
167    cptpEEKS - Kinetic energy state variables.
168    cptpEENH - Energy history state variables.
169    cptpEDFH - free energy history variables.
170 */
171
172
173 static const char *st_names(int cptp,int ecpt)
174 {
175     switch (cptp)
176     {
177     case cptpEST: return est_names [ecpt]; break;
178     case cptpEEKS: return eeks_names[ecpt]; break;
179     case cptpEENH: return eenh_names[ecpt]; break;
180     case cptpEDFH: return edfh_names[ecpt]; break;
181     }
182
183     return NULL;
184 }
185
186 static void cp_warning(FILE *fp)
187 {
188     fprintf(fp,"\nWARNING: Checkpoint file is corrupted or truncated\n\n");
189 }
190
191 static void cp_error()
192 {
193     gmx_fatal(FARGS,"Checkpoint file corrupted/truncated, or maybe you are out of disk space?");
194 }
195
196 static void do_cpt_string_err(XDR *xd,gmx_bool bRead,const char *desc,char **s,FILE *list)
197 {
198     bool_t res=0;
199     
200     if (bRead)
201     {
202         snew(*s,CPTSTRLEN);
203     }
204     res = xdr_string(xd,s,CPTSTRLEN);
205     if (res == 0)
206     {
207         cp_error();
208     }
209     if (list)
210     {
211         fprintf(list,"%s = %s\n",desc,*s);
212         sfree(*s);
213     }
214 }
215
216 static int do_cpt_int(XDR *xd,const char *desc,int *i,FILE *list)
217 {
218     bool_t res=0;
219     
220     res = xdr_int(xd,i);
221     if (res == 0)
222     {
223         return -1;
224     }
225     if (list)
226     {
227         fprintf(list,"%s = %d\n",desc,*i);
228     }
229     return 0;
230 }
231
232 static int do_cpt_u_chars(XDR *xd,const char *desc,int n,unsigned char *i,FILE *list)
233 {
234     bool_t res=1;
235     int j;
236     if (list)
237     {
238         fprintf(list,"%s = ",desc);
239     }
240     for (j=0; j<n && res; j++)
241     {
242         res &= xdr_u_char(xd,&i[j]);
243         if (list)
244         {
245             fprintf(list,"%02x",i[j]);
246         }
247     }
248     if (list)
249     {
250         fprintf(list,"\n");
251     }
252     if (res == 0)
253     {
254         return -1;
255     }
256
257     return 0;
258 }
259
260 static void do_cpt_int_err(XDR *xd,const char *desc,int *i,FILE *list)
261 {
262     if (do_cpt_int(xd,desc,i,list) < 0)
263     {
264         cp_error();
265     }
266 }
267
268 static void do_cpt_step_err(XDR *xd,const char *desc,gmx_large_int_t *i,FILE *list)
269 {
270     bool_t res=0;
271     char   buf[STEPSTRSIZE];
272
273     res = xdr_gmx_large_int(xd,i,"reading checkpoint file");
274     if (res == 0)
275     {
276         cp_error();
277     }
278     if (list)
279     {
280         fprintf(list,"%s = %s\n",desc,gmx_step_str(*i,buf));
281     }
282 }
283
284 static void do_cpt_double_err(XDR *xd,const char *desc,double *f,FILE *list)
285 {
286     bool_t res=0;
287     
288     res = xdr_double(xd,f);
289     if (res == 0)
290     {
291         cp_error();
292     }
293     if (list)
294     {
295         fprintf(list,"%s = %f\n",desc,*f);
296     }
297 }
298
299 /* If nval >= 0, nval is used; on read this should match the passed value.
300  * If nval n<0, *nptr is used; on read the value is stored in nptr
301  */
302 static int do_cpte_reals_low(XDR *xd,int cptp,int ecpt,int sflags,
303                              int nval,int *nptr,real **v,
304                              FILE *list,int erealtype)
305 {
306     bool_t res=0;
307 #ifndef GMX_DOUBLE
308     int  dtc=xdr_datatype_float; 
309 #else
310     int  dtc=xdr_datatype_double;
311 #endif
312     real *vp,*va=NULL;
313     float  *vf;
314     double *vd;
315     int  nf,dt,i;
316     
317     if (list == NULL)
318     {
319         if (nval >= 0)
320         {
321             nf = nval;
322         }
323         else
324         {
325         if (nptr == NULL)
326         {
327             gmx_incons("*ntpr=NULL in do_cpte_reals_low");
328         }
329         nf = *nptr;
330         }
331     }
332     res = xdr_int(xd,&nf);
333     if (res == 0)
334     {
335         return -1;
336     }
337     if (list == NULL)
338     {
339         if (nval >= 0)
340         {
341             if (nf != nval)
342             {
343                 gmx_fatal(FARGS,"Count mismatch for state entry %s, code count is %d, file count is %d\n",st_names(cptp,ecpt),nval,nf);
344             }
345         }
346         else
347         {
348             *nptr = nf;
349         }
350     }
351     dt = dtc;
352     res = xdr_int(xd,&dt);
353     if (res == 0)
354     {
355         return -1;
356     }
357     if (dt != dtc)
358     {
359         fprintf(stderr,"Precision mismatch for state entry %s, code precision is %s, file precision is %s\n",
360                 st_names(cptp,ecpt),xdr_datatype_names[dtc],
361                 xdr_datatype_names[dt]);
362     }
363     if (list || !(sflags & (1<<ecpt)))
364     {
365         snew(va,nf);
366         vp = va;
367     }
368     else
369     {
370         if (*v == NULL)
371         {
372             snew(*v,nf);
373         }
374         vp = *v;
375     }
376     if (dt == xdr_datatype_float)
377     {
378         if (dtc == xdr_datatype_float)
379         {
380             vf = (float *)vp;
381         }
382         else
383         {
384             snew(vf,nf);
385         }
386         res = xdr_vector(xd,(char *)vf,nf,
387                          (unsigned int)sizeof(float),(xdrproc_t)xdr_float);
388         if (res == 0)
389         {
390             return -1;
391         }
392         if (dtc != xdr_datatype_float)
393         {
394             for(i=0; i<nf; i++)
395             {
396                 vp[i] = vf[i];
397             }
398             sfree(vf);
399         }
400     }
401     else
402     {
403         if (dtc == xdr_datatype_double)
404         {
405             vd = (double *)vp;
406         }
407         else
408         {
409             snew(vd,nf);
410         }
411         res = xdr_vector(xd,(char *)vd,nf,
412                          (unsigned int)sizeof(double),(xdrproc_t)xdr_double);
413         if (res == 0)
414         {
415             return -1;
416         }
417         if (dtc != xdr_datatype_double)
418         {
419             for(i=0; i<nf; i++)
420             {
421                 vp[i] = vd[i];
422             }
423             sfree(vd);
424         }
425     }
426     
427     if (list)
428     {
429         switch (erealtype)
430         {
431         case ecprREAL:
432             pr_reals(list,0,st_names(cptp,ecpt),vp,nf);
433             break;
434         case ecprRVEC:
435             pr_rvecs(list,0,st_names(cptp,ecpt),(rvec *)vp,nf/3);
436             break;
437         default:
438             gmx_incons("Unknown checkpoint real type");
439         }
440     }
441     if (va)
442     {
443         sfree(va);
444     }
445
446     return 0;
447 }
448
449
450 /* This function stores n along with the reals for reading,
451  * but on reading it assumes that n matches the value in the checkpoint file,
452  * a fatal error is generated when this is not the case.
453  */
454 static int do_cpte_reals(XDR *xd,int cptp,int ecpt,int sflags,
455                          int n,real **v,FILE *list)
456 {
457     return do_cpte_reals_low(xd,cptp,ecpt,sflags,n,NULL,v,list,ecprREAL);
458 }
459
460 /* This function does the same as do_cpte_reals,
461  * except that on reading it ignores the passed value of *n
462  * and stored the value read from the checkpoint file in *n.
463  */
464 static int do_cpte_n_reals(XDR *xd,int cptp,int ecpt,int sflags,
465                            int *n,real **v,FILE *list)
466 {
467     return do_cpte_reals_low(xd,cptp,ecpt,sflags,-1,n,v,list,ecprREAL);
468 }
469
470 static int do_cpte_real(XDR *xd,int cptp,int ecpt,int sflags,
471                         real *r,FILE *list)
472 {
473     int n;
474
475     return do_cpte_reals_low(xd,cptp,ecpt,sflags,1,NULL,&r,list,ecprREAL);
476 }
477
478 static int do_cpte_ints(XDR *xd,int cptp,int ecpt,int sflags,
479                         int n,int **v,FILE *list)
480 {
481     bool_t res=0;
482     int  dtc=xdr_datatype_int;
483     int *vp,*va=NULL;
484     int  nf,dt,i;
485     
486     nf = n;
487     res = xdr_int(xd,&nf);
488     if (res == 0)
489     {
490         return -1;
491     }
492     if (list == NULL && v != NULL && nf != n)
493     {
494         gmx_fatal(FARGS,"Count mismatch for state entry %s, code count is %d, file count is %d\n",st_names(cptp,ecpt),n,nf);
495     }
496     dt = dtc;
497     res = xdr_int(xd,&dt);
498     if (res == 0)
499     {
500         return -1;
501     }
502     if (dt != dtc)
503     {
504         gmx_fatal(FARGS,"Type mismatch for state entry %s, code type is %s, file type is %s\n",
505                   st_names(cptp,ecpt),xdr_datatype_names[dtc],
506                   xdr_datatype_names[dt]);
507     }
508     if (list || !(sflags & (1<<ecpt)) || v == NULL)
509     {
510         snew(va,nf);
511         vp = va;
512     }
513     else
514     {
515         if (*v == NULL)
516         {
517             snew(*v,nf);
518         }
519         vp = *v;
520     }
521     res = xdr_vector(xd,(char *)vp,nf,
522                      (unsigned int)sizeof(int),(xdrproc_t)xdr_int);
523     if (res == 0)
524     {
525         return -1;
526     }
527     if (list)
528     {
529         pr_ivec(list,0,st_names(cptp,ecpt),vp,nf,TRUE);
530     }
531     if (va)
532     {
533         sfree(va);
534     }
535
536     return 0;
537 }
538
539 static int do_cpte_int(XDR *xd,int cptp,int ecpt,int sflags,
540                        int *i,FILE *list)
541 {
542     return do_cpte_ints(xd,cptp,ecpt,sflags,1,&i,list);
543 }
544
545 static int do_cpte_doubles(XDR *xd,int cptp,int ecpt,int sflags,
546                            int n,double **v,FILE *list)
547 {
548     bool_t res=0;
549     int  dtc=xdr_datatype_double;
550     double *vp,*va=NULL;
551     int  nf,dt,i;
552     
553     nf = n;
554     res = xdr_int(xd,&nf);
555     if (res == 0)
556     {
557         return -1;
558     }
559     if (list == NULL && nf != n)
560     {
561         gmx_fatal(FARGS,"Count mismatch for state entry %s, code count is %d, file count is %d\n",st_names(cptp,ecpt),n,nf);
562     }
563     dt = dtc;
564     res = xdr_int(xd,&dt);
565     if (res == 0)
566     {
567         return -1;
568     }
569     if (dt != dtc)
570     {
571         gmx_fatal(FARGS,"Precision mismatch for state entry %s, code precision is %s, file precision is %s\n",
572                   st_names(cptp,ecpt),xdr_datatype_names[dtc],
573                   xdr_datatype_names[dt]);
574     }
575     if (list || !(sflags & (1<<ecpt)))
576     {
577         snew(va,nf);
578         vp = va;
579     }
580     else
581     {
582         if (*v == NULL)
583         {
584             snew(*v,nf);
585         }
586         vp = *v;
587     }
588     res = xdr_vector(xd,(char *)vp,nf,
589                      (unsigned int)sizeof(double),(xdrproc_t)xdr_double);
590     if (res == 0)
591     {
592         return -1;
593     }
594     if (list)
595     {
596         pr_doubles(list,0,st_names(cptp,ecpt),vp,nf);
597     }
598     if (va)
599     {
600         sfree(va);
601     }
602
603     return 0;
604 }
605
606 static int do_cpte_double(XDR *xd,int cptp,int ecpt,int sflags,
607                           double *r,FILE *list)
608 {
609     return do_cpte_doubles(xd,cptp,ecpt,sflags,1,&r,list);
610 }
611
612
613 static int do_cpte_rvecs(XDR *xd,int cptp,int ecpt,int sflags,
614                          int n,rvec **v,FILE *list)
615 {
616     int n3;
617
618     return do_cpte_reals_low(xd,cptp,ecpt,sflags,
619                              n*DIM,NULL,(real **)v,list,ecprRVEC);
620 }
621
622 static int do_cpte_matrix(XDR *xd,int cptp,int ecpt,int sflags,
623                           matrix v,FILE *list)
624 {
625     real *vr;
626     real ret;
627
628     vr = (real *)&(v[0][0]);
629     ret = do_cpte_reals_low(xd,cptp,ecpt,sflags,
630                             DIM*DIM,NULL,&vr,NULL,ecprMATRIX);
631     
632     if (list && ret == 0)
633     {
634         pr_rvecs(list,0,st_names(cptp,ecpt),v,DIM);
635     }
636     
637     return ret;
638 }
639
640
641 static int do_cpte_nmatrix(XDR *xd,int cptp,int ecpt,int sflags,
642                            int n, real **v,FILE *list)
643 {
644     int i;
645     real *vr;
646     real ret,reti;
647     char name[CPTSTRLEN];
648
649     ret = 0;
650     if (v==NULL)
651     {
652         snew(v,n);
653     }
654     for (i=0;i<n;i++)
655     {
656         reti = 0;
657         vr = v[i];
658         reti = do_cpte_reals_low(xd,cptp,ecpt,sflags,n,NULL,&(v[i]),NULL,ecprREAL);
659         if (list && reti == 0)
660         {
661             sprintf(name,"%s[%d]",st_names(cptp,ecpt),i);
662             pr_reals(list,0,name,v[i],n);
663         }
664         if (reti == 0)
665         {
666             ret = 0;
667         }
668     }
669     return ret;
670 }
671
672 static int do_cpte_matrices(XDR *xd,int cptp,int ecpt,int sflags,
673                             int n,matrix **v,FILE *list)
674 {
675     bool_t res=0;
676     matrix *vp,*va=NULL;
677     real *vr;
678     int  nf,i,j,k;
679     int  ret;
680
681     nf = n;
682     res = xdr_int(xd,&nf);
683     if (res == 0)
684     {
685         return -1;
686     }
687     if (list == NULL && nf != n)
688     {
689         gmx_fatal(FARGS,"Count mismatch for state entry %s, code count is %d, file count is %d\n",st_names(cptp,ecpt),n,nf);
690     }
691     if (list || !(sflags & (1<<ecpt)))
692     {
693         snew(va,nf);
694         vp = va;
695     }
696     else
697     {
698         if (*v == NULL)
699         {
700             snew(*v,nf);
701         }
702         vp = *v;
703     }
704     snew(vr,nf*DIM*DIM);
705     for(i=0; i<nf; i++)
706     {
707         for(j=0; j<DIM; j++)
708         {
709             for(k=0; k<DIM; k++)
710             {
711                 vr[(i*DIM+j)*DIM+k] = vp[i][j][k];
712             }
713         }
714     }
715     ret = do_cpte_reals_low(xd,cptp,ecpt,sflags,
716                             nf*DIM*DIM,NULL,&vr,NULL,ecprMATRIX);
717     for(i=0; i<nf; i++)
718     {
719         for(j=0; j<DIM; j++)
720         {
721             for(k=0; k<DIM; k++)
722             {
723                 vp[i][j][k] = vr[(i*DIM+j)*DIM+k];
724             }
725         }
726     }
727     sfree(vr);
728     
729     if (list && ret == 0)
730     {
731         for(i=0; i<nf; i++)
732         {
733             pr_rvecs(list,0,st_names(cptp,ecpt),vp[i],DIM);
734         }
735     }
736     if (va)
737     {
738         sfree(va);
739     }
740     
741     return ret;
742 }
743
744 static void do_cpt_header(XDR *xd,gmx_bool bRead,int *file_version,
745                           char **version,char **btime,char **buser,char **bhost,
746                           int *double_prec,
747                           char **fprog,char **ftime,
748                           int *eIntegrator,int *simulation_part,
749                           gmx_large_int_t *step,double *t,
750                           int *nnodes,int *dd_nc,int *npme,
751                           int *natoms,int *ngtc, int *nnhpres, int *nhchainlength,
752                           int *nlambda, int *flags_state,
753                           int *flags_eks,int *flags_enh, int *flags_dfh,
754                           FILE *list)
755 {
756     bool_t res=0;
757     int  magic;
758     int  idum=0;
759     int  i;
760     char *fhost;
761
762     if (bRead)
763     {
764         magic = -1;
765     }
766     else
767     {
768         magic = CPT_MAGIC1;
769     }
770     res = xdr_int(xd,&magic);
771     if (res == 0)
772     {
773         gmx_fatal(FARGS,"The checkpoint file is empty/corrupted, or maybe you are out of disk space?");
774     }
775     if (magic != CPT_MAGIC1)
776     {
777         gmx_fatal(FARGS,"Start of file magic number mismatch, checkpoint file has %d, should be %d\n"
778                   "The checkpoint file is corrupted or not a checkpoint file",
779                   magic,CPT_MAGIC1);
780     }
781     if (!bRead)
782     {
783         snew(fhost,255);
784 #ifdef HAVE_UNISTD_H
785         if (gethostname(fhost,255) != 0)
786         {
787             sprintf(fhost,"unknown");
788         }
789 #else
790         sprintf(fhost,"unknown");
791 #endif  
792     }
793     do_cpt_string_err(xd,bRead,"GROMACS version"           ,version,list);
794     do_cpt_string_err(xd,bRead,"GROMACS build time"        ,btime,list);
795     do_cpt_string_err(xd,bRead,"GROMACS build user"        ,buser,list);
796     do_cpt_string_err(xd,bRead,"GROMACS build host"        ,bhost,list);
797     do_cpt_string_err(xd,bRead,"generating program"        ,fprog,list);
798     do_cpt_string_err(xd,bRead,"generation time"           ,ftime,list);
799     *file_version = cpt_version;
800     do_cpt_int_err(xd,"checkpoint file version",file_version,list);
801     if (*file_version > cpt_version)
802     {
803         gmx_fatal(FARGS,"Attempting to read a checkpoint file of version %d with code of version %d\n",*file_version,cpt_version);
804     }
805     if (*file_version >= 13)
806     {
807         do_cpt_int_err(xd,"GROMACS double precision",double_prec,list);
808     }
809     else
810     {
811         *double_prec = -1;
812     }
813     if (*file_version >= 12)
814     {
815         do_cpt_string_err(xd,bRead,"generating host"           ,&fhost,list);
816         if (list == NULL)
817         {
818             sfree(fhost);
819         }
820     }
821     do_cpt_int_err(xd,"#atoms"            ,natoms     ,list);
822     do_cpt_int_err(xd,"#T-coupling groups",ngtc       ,list);
823     if (*file_version >= 10) 
824     {
825         do_cpt_int_err(xd,"#Nose-Hoover T-chains",nhchainlength,list);
826     }
827     else
828     {
829         *nhchainlength = 1;
830     }
831     if (*file_version >= 11)
832     {
833         do_cpt_int_err(xd,"#Nose-Hoover T-chains for barostat ",nnhpres,list);
834     }
835     else
836     {
837         *nnhpres = 0;
838     }
839     if (*file_version >= 14)
840     {
841         do_cpt_int_err(xd,"# of total lambda states ",nlambda,list);
842     }
843     else
844     {
845         *nlambda = 0;
846     }
847     do_cpt_int_err(xd,"integrator"        ,eIntegrator,list);
848         if (*file_version >= 3)
849         {
850                 do_cpt_int_err(xd,"simulation part #", simulation_part,list);
851         }
852         else
853         {
854                 *simulation_part = 1;
855         }
856     if (*file_version >= 5)
857     {
858         do_cpt_step_err(xd,"step"         ,step       ,list);
859     }
860     else
861     {
862         do_cpt_int_err(xd,"step"          ,&idum      ,list);
863         *step = idum;
864     }
865     do_cpt_double_err(xd,"t"              ,t          ,list);
866     do_cpt_int_err(xd,"#PP-nodes"         ,nnodes     ,list);
867     idum = 1;
868     do_cpt_int_err(xd,"dd_nc[x]",dd_nc ? &(dd_nc[0]) : &idum,list);
869     do_cpt_int_err(xd,"dd_nc[y]",dd_nc ? &(dd_nc[1]) : &idum,list);
870     do_cpt_int_err(xd,"dd_nc[z]",dd_nc ? &(dd_nc[2]) : &idum,list);
871     do_cpt_int_err(xd,"#PME-only nodes",npme,list);
872     do_cpt_int_err(xd,"state flags",flags_state,list);
873         if (*file_version >= 4)
874     {
875         do_cpt_int_err(xd,"ekin data flags",flags_eks,list);
876         do_cpt_int_err(xd,"energy history flags",flags_enh,list);
877     }
878     else
879     {
880         *flags_eks  = 0;
881         *flags_enh   = (*flags_state >> (estORIRE_DTAV+1));
882         *flags_state = (*flags_state & ~((1<<(estORIRE_DTAV+1)) |
883                                          (1<<(estORIRE_DTAV+2)) |
884                                          (1<<(estORIRE_DTAV+3))));
885     }
886         if (*file_version >= 14)
887     {
888         do_cpt_int_err(xd,"df history flags",flags_dfh,list);
889     } else {
890         *flags_dfh = 0;
891     }
892 }
893
894 static int do_cpt_footer(XDR *xd,gmx_bool bRead,int file_version)
895 {
896     bool_t res=0;
897     int  magic;
898     
899     if (file_version >= 2)
900     {
901         magic = CPT_MAGIC2;
902         res = xdr_int(xd,&magic);
903         if (res == 0)
904         {
905             cp_error();
906         }
907         if (magic != CPT_MAGIC2)
908         {
909             return -1;
910         }
911     }
912
913     return 0;
914 }
915
916 static int do_cpt_state(XDR *xd,gmx_bool bRead,
917                         int fflags,t_state *state,
918                         gmx_bool bReadRNG,FILE *list)
919 {
920     int  sflags;
921     int  **rng_p,**rngi_p;
922     int  i;
923     int  ret;
924     int  nnht,nnhtp;
925
926     ret = 0;
927     
928     nnht = state->nhchainlength*state->ngtc;
929     nnhtp = state->nhchainlength*state->nnhpres;
930
931     if (bReadRNG)
932     {
933         rng_p  = (int **)&state->ld_rng;
934         rngi_p = &state->ld_rngi;
935     }
936     else
937     {
938         /* Do not read the RNG data */
939         rng_p  = NULL;
940         rngi_p = NULL;
941     }
942     /* We want the MC_RNG the same across all the notes for now -- lambda MC is global */
943
944     sflags = state->flags;
945     for(i=0; (i<estNR && ret == 0); i++)
946     {
947         if (fflags & (1<<i))
948         {
949             switch (i)
950             {
951             case estLAMBDA:  ret = do_cpte_reals(xd,cptpEST,i,sflags,efptNR,&(state->lambda),list); break;
952             case estFEPSTATE: ret = do_cpte_int (xd,cptpEST,i,sflags,&state->fep_state,list); break;
953             case estBOX:     ret = do_cpte_matrix(xd,cptpEST,i,sflags,state->box,list); break;
954             case estBOX_REL: ret = do_cpte_matrix(xd,cptpEST,i,sflags,state->box_rel,list); break;
955             case estBOXV:    ret = do_cpte_matrix(xd,cptpEST,i,sflags,state->boxv,list); break;
956             case estPRES_PREV: ret = do_cpte_matrix(xd,cptpEST,i,sflags,state->pres_prev,list); break;
957             case estSVIR_PREV:  ret = do_cpte_matrix(xd,cptpEST,i,sflags,state->svir_prev,list); break;
958             case estFVIR_PREV:  ret = do_cpte_matrix(xd,cptpEST,i,sflags,state->fvir_prev,list); break;
959             case estNH_XI:   ret = do_cpte_doubles(xd,cptpEST,i,sflags,nnht,&state->nosehoover_xi,list); break;
960             case estNH_VXI:  ret = do_cpte_doubles(xd,cptpEST,i,sflags,nnht,&state->nosehoover_vxi,list); break;
961             case estNHPRES_XI:   ret = do_cpte_doubles(xd,cptpEST,i,sflags,nnhtp,&state->nhpres_xi,list); break;
962             case estNHPRES_VXI:  ret = do_cpte_doubles(xd,cptpEST,i,sflags,nnhtp,&state->nhpres_vxi,list); break;
963             case estTC_INT:  ret = do_cpte_doubles(xd,cptpEST,i,sflags,state->ngtc,&state->therm_integral,list); break;
964             case estVETA:    ret = do_cpte_real(xd,cptpEST,i,sflags,&state->veta,list); break;
965             case estVOL0:    ret = do_cpte_real(xd,cptpEST,i,sflags,&state->vol0,list); break;
966             case estX:       ret = do_cpte_rvecs(xd,cptpEST,i,sflags,state->natoms,&state->x,list); break;
967             case estV:       ret = do_cpte_rvecs(xd,cptpEST,i,sflags,state->natoms,&state->v,list); break;
968             case estSDX:     ret = do_cpte_rvecs(xd,cptpEST,i,sflags,state->natoms,&state->sd_X,list); break;
969             case estLD_RNG:  ret = do_cpte_ints(xd,cptpEST,i,sflags,state->nrng,rng_p,list); break;
970             case estLD_RNGI: ret = do_cpte_ints(xd,cptpEST,i,sflags,state->nrngi,rngi_p,list); break;
971             case estMC_RNG:  ret = do_cpte_ints(xd,cptpEST,i,sflags,state->nmcrng,(int **)&state->mc_rng,list); break;
972             case estMC_RNGI: ret = do_cpte_ints(xd,cptpEST,i,sflags,1,&state->mc_rngi,list); break;
973             case estDISRE_INITF:  ret = do_cpte_real (xd,cptpEST,i,sflags,&state->hist.disre_initf,list); break;
974             case estDISRE_RM3TAV: ret = do_cpte_reals(xd,cptpEST,i,sflags,state->hist.ndisrepairs,&state->hist.disre_rm3tav,list); break;
975             case estORIRE_INITF:  ret = do_cpte_real (xd,cptpEST,i,sflags,&state->hist.orire_initf,list); break;
976             case estORIRE_DTAV:   ret = do_cpte_reals(xd,cptpEST,i,sflags,state->hist.norire_Dtav,&state->hist.orire_Dtav,list); break;
977             default:
978                 gmx_fatal(FARGS,"Unknown state entry %d\n"
979                           "You are probably reading a new checkpoint file with old code",i);
980             }
981         }
982     }
983     
984     return ret;
985 }
986
987 static int do_cpt_ekinstate(XDR *xd,gmx_bool bRead,
988                             int fflags,ekinstate_t *ekins,
989                             FILE *list)
990 {
991     int  i;
992     int  ret;
993
994     ret = 0;
995
996     for(i=0; (i<eeksNR && ret == 0); i++)
997     {
998         if (fflags & (1<<i))
999         {
1000             switch (i)
1001             {
1002                 
1003                         case eeksEKIN_N:     ret = do_cpte_int(xd,cptpEEKS,i,fflags,&ekins->ekin_n,list); break;
1004                         case eeksEKINH :     ret = do_cpte_matrices(xd,cptpEEKS,i,fflags,ekins->ekin_n,&ekins->ekinh,list); break;
1005                         case eeksEKINF:      ret = do_cpte_matrices(xd,cptpEEKS,i,fflags,ekins->ekin_n,&ekins->ekinf,list); break;
1006                         case eeksEKINO:      ret = do_cpte_matrices(xd,cptpEEKS,i,fflags,ekins->ekin_n,&ekins->ekinh_old,list); break;
1007             case eeksEKINTOTAL:  ret = do_cpte_matrix(xd,cptpEEKS,i,fflags,ekins->ekin_total,list); break;
1008             case eeksEKINSCALEF: ret = do_cpte_doubles(xd,cptpEEKS,i,fflags,ekins->ekin_n,&ekins->ekinscalef_nhc,list); break;
1009             case eeksVSCALE:     ret = do_cpte_doubles(xd,1,cptpEEKS,fflags,ekins->ekin_n,&ekins->vscale_nhc,list); break;
1010             case eeksEKINSCALEH: ret = do_cpte_doubles(xd,1,cptpEEKS,fflags,ekins->ekin_n,&ekins->ekinscaleh_nhc,list); break;
1011                         case eeksDEKINDL :   ret = do_cpte_real(xd,1,cptpEEKS,fflags,&ekins->dekindl,list); break;
1012             case eeksMVCOS:      ret = do_cpte_real(xd,1,cptpEEKS,fflags,&ekins->mvcos,list); break;
1013             default:
1014                 gmx_fatal(FARGS,"Unknown ekin data state entry %d\n"
1015                           "You are probably reading a new checkpoint file with old code",i);
1016             }
1017         }
1018     }
1019     
1020     return ret;
1021 }
1022
1023
1024 static int do_cpt_enerhist(XDR *xd,gmx_bool bRead,
1025                            int fflags,energyhistory_t *enerhist,
1026                            FILE *list)
1027 {
1028     int  i;
1029     int  j;
1030     int  ret;
1031
1032     ret = 0;
1033
1034     if (bRead)
1035     {
1036         enerhist->nsteps     = 0;
1037         enerhist->nsum       = 0;
1038         enerhist->nsteps_sim = 0;
1039         enerhist->nsum_sim   = 0;
1040         enerhist->dht        = NULL;
1041
1042         if (fflags & (1<< eenhENERGY_DELTA_H_NN) )
1043         {
1044             snew(enerhist->dht,1);
1045             enerhist->dht->ndh = NULL;
1046             enerhist->dht->dh = NULL;
1047             enerhist->dht->start_lambda_set=FALSE;
1048         }
1049     }
1050
1051     for(i=0; (i<eenhNR && ret == 0); i++)
1052     {
1053         if (fflags & (1<<i))
1054         {
1055             switch (i)
1056             {
1057                         case eenhENERGY_N:     ret = do_cpte_int(xd,cptpEENH,i,fflags,&enerhist->nener,list); break;
1058                         case eenhENERGY_AVER:  ret = do_cpte_doubles(xd,cptpEENH,i,fflags,enerhist->nener,&enerhist->ener_ave,list); break;
1059                         case eenhENERGY_SUM:   ret = do_cpte_doubles(xd,cptpEENH,i,fflags,enerhist->nener,&enerhist->ener_sum,list); break;
1060             case eenhENERGY_NSUM:  do_cpt_step_err(xd,eenh_names[i],&enerhist->nsum,list); break;
1061             case eenhENERGY_SUM_SIM: ret = do_cpte_doubles(xd,cptpEENH,i,fflags,enerhist->nener,&enerhist->ener_sum_sim,list); break;
1062             case eenhENERGY_NSUM_SIM:   do_cpt_step_err(xd,eenh_names[i],&enerhist->nsum_sim,list); break;
1063             case eenhENERGY_NSTEPS:     do_cpt_step_err(xd,eenh_names[i],&enerhist->nsteps,list); break;
1064             case eenhENERGY_NSTEPS_SIM: do_cpt_step_err(xd,eenh_names[i],&enerhist->nsteps_sim,list); break;
1065             case eenhENERGY_DELTA_H_NN: do_cpt_int_err(xd,eenh_names[i], &(enerhist->dht->nndh), list);
1066                 if (bRead) /* now allocate memory for it */
1067                 {
1068                     snew(enerhist->dht->dh, enerhist->dht->nndh);
1069                     snew(enerhist->dht->ndh, enerhist->dht->nndh);
1070                     for(j=0;j<enerhist->dht->nndh;j++)
1071                     {
1072                         enerhist->dht->ndh[j] = 0;
1073                         enerhist->dht->dh[j] = NULL;
1074                     }
1075                 }
1076                 break;
1077             case eenhENERGY_DELTA_H_LIST:
1078                 for(j=0;j<enerhist->dht->nndh;j++)
1079                 {
1080                     ret=do_cpte_n_reals(xd, cptpEENH, i, fflags, &enerhist->dht->ndh[j], &(enerhist->dht->dh[j]), list);
1081                 }
1082                 break;
1083             case eenhENERGY_DELTA_H_STARTTIME:
1084                 ret=do_cpte_double(xd, cptpEENH, i, fflags, &(enerhist->dht->start_time), list); break;
1085             case eenhENERGY_DELTA_H_STARTLAMBDA:
1086                 ret=do_cpte_double(xd, cptpEENH, i, fflags, &(enerhist->dht->start_lambda), list); break;
1087             default:
1088                 gmx_fatal(FARGS,"Unknown energy history entry %d\n"
1089                           "You are probably reading a new checkpoint file with old code",i);
1090             }
1091         }
1092     }
1093
1094     if ((fflags & (1<<eenhENERGY_SUM)) && !(fflags & (1<<eenhENERGY_SUM_SIM)))
1095     {
1096         /* Assume we have an old file format and copy sum to sum_sim */
1097         srenew(enerhist->ener_sum_sim,enerhist->nener);
1098         for(i=0; i<enerhist->nener; i++)
1099         {
1100             enerhist->ener_sum_sim[i] = enerhist->ener_sum[i];
1101         }
1102         fflags |= (1<<eenhENERGY_SUM_SIM);
1103     }
1104     
1105     if ( (fflags & (1<<eenhENERGY_NSUM)) &&
1106         !(fflags & (1<<eenhENERGY_NSTEPS)))
1107     {
1108         /* Assume we have an old file format and copy nsum to nsteps */
1109         enerhist->nsteps = enerhist->nsum;
1110         fflags |= (1<<eenhENERGY_NSTEPS);
1111     }
1112     if ( (fflags & (1<<eenhENERGY_NSUM_SIM)) &&
1113         !(fflags & (1<<eenhENERGY_NSTEPS_SIM)))
1114     {
1115         /* Assume we have an old file format and copy nsum to nsteps */
1116         enerhist->nsteps_sim = enerhist->nsum_sim;
1117         fflags |= (1<<eenhENERGY_NSTEPS_SIM);
1118     }
1119
1120     return ret;
1121 }
1122
1123 static int do_cpt_df_hist(XDR *xd,gmx_bool bRead,int fflags,df_history_t *dfhist,FILE *list)
1124 {
1125     int  i,nlambda;
1126     int  ret;
1127
1128     nlambda = dfhist->nlambda;
1129     ret = 0;
1130
1131     for(i=0; (i<edfhNR && ret == 0); i++)
1132     {
1133         if (fflags & (1<<i))
1134         {
1135             switch (i)
1136             {
1137             case edfhBEQUIL:       ret = do_cpte_int(xd,cptpEDFH,i,fflags,&dfhist->bEquil,list); break;
1138             case edfhNATLAMBDA:    ret = do_cpte_ints(xd,cptpEDFH,i,fflags,nlambda,&dfhist->n_at_lam,list); break;
1139             case edfhWLHISTO:      ret = do_cpte_reals(xd,cptpEDFH,i,fflags,nlambda,&dfhist->wl_histo,list); break;
1140             case edfhWLDELTA:      ret = do_cpte_real(xd,cptpEDFH,i,fflags,&dfhist->wl_delta,list); break;
1141             case edfhSUMWEIGHTS:   ret = do_cpte_reals(xd,cptpEDFH,i,fflags,nlambda,&dfhist->sum_weights,list); break;
1142             case edfhSUMDG:        ret = do_cpte_reals(xd,cptpEDFH,i,fflags,nlambda,&dfhist->sum_dg,list); break;
1143             case edfhSUMMINVAR:    ret = do_cpte_reals(xd,cptpEDFH,i,fflags,nlambda,&dfhist->sum_minvar,list); break;
1144             case edfhSUMVAR:       ret = do_cpte_reals(xd,cptpEDFH,i,fflags,nlambda,&dfhist->sum_variance,list); break;
1145             case edfhACCUMP:       ret = do_cpte_nmatrix(xd,cptpEDFH,i,fflags,nlambda,dfhist->accum_p,list); break;
1146             case edfhACCUMM:       ret = do_cpte_nmatrix(xd,cptpEDFH,i,fflags,nlambda,dfhist->accum_m,list); break;
1147             case edfhACCUMP2:      ret = do_cpte_nmatrix(xd,cptpEDFH,i,fflags,nlambda,dfhist->accum_p2,list); break;
1148             case edfhACCUMM2:      ret = do_cpte_nmatrix(xd,cptpEDFH,i,fflags,nlambda,dfhist->accum_m2,list); break;
1149             case edfhTIJ:          ret = do_cpte_nmatrix(xd,cptpEDFH,i,fflags,nlambda,dfhist->Tij,list); break;
1150             case edfhTIJEMP:       ret = do_cpte_nmatrix(xd,cptpEDFH,i,fflags,nlambda,dfhist->Tij_empirical,list); break;
1151
1152             default:
1153                 gmx_fatal(FARGS,"Unknown df history entry %d\n"
1154                           "You are probably reading a new checkpoint file with old code",i);
1155             }
1156         }
1157     }
1158
1159     return ret;
1160 }
1161
1162 static int do_cpt_files(XDR *xd, gmx_bool bRead, 
1163                         gmx_file_position_t **p_outputfiles, int *nfiles, 
1164                         FILE *list, int file_version)
1165 {
1166     int    i,j;
1167     gmx_off_t  offset;
1168     gmx_off_t  mask = 0xFFFFFFFFL;
1169     int    offset_high,offset_low;
1170     char   *buf;
1171     gmx_file_position_t *outputfiles;
1172
1173     if (do_cpt_int(xd,"number of output files",nfiles,list) != 0)
1174     {
1175         return -1;
1176     }
1177
1178     if(bRead)
1179     {
1180         snew(*p_outputfiles,*nfiles);
1181     }
1182
1183     outputfiles = *p_outputfiles;
1184
1185     for(i=0;i<*nfiles;i++)
1186     {
1187         /* 64-bit XDR numbers are not portable, so it is stored as separate high/low fractions */
1188         if(bRead)
1189         {
1190             do_cpt_string_err(xd,bRead,"output filename",&buf,list);
1191             strncpy(outputfiles[i].filename,buf,CPTSTRLEN-1);
1192             if(list==NULL)
1193             {
1194                 sfree(buf);                     
1195             }
1196
1197             if (do_cpt_int(xd,"file_offset_high",&offset_high,list) != 0)
1198             {
1199                 return -1;
1200             }
1201             if (do_cpt_int(xd,"file_offset_low",&offset_low,list) != 0)
1202             {
1203                 return -1;
1204             }
1205 #if (SIZEOF_GMX_OFF_T > 4)
1206             outputfiles[i].offset = ( ((gmx_off_t) offset_high) << 32 ) | ( (gmx_off_t) offset_low & mask );
1207 #else
1208             outputfiles[i].offset = offset_low;
1209 #endif
1210         }
1211         else
1212         {
1213             buf = outputfiles[i].filename;
1214             do_cpt_string_err(xd,bRead,"output filename",&buf,list);
1215             /* writing */
1216             offset      = outputfiles[i].offset;
1217             if (offset == -1)
1218             {
1219                 offset_low  = -1;
1220                 offset_high = -1;
1221             }
1222             else
1223             {
1224 #if (SIZEOF_GMX_OFF_T > 4)
1225                 offset_low  = (int) (offset & mask);
1226                 offset_high = (int) ((offset >> 32) & mask);
1227 #else
1228                 offset_low  = offset;
1229                 offset_high = 0;
1230 #endif
1231             }
1232             if (do_cpt_int(xd,"file_offset_high",&offset_high,list) != 0)
1233             {
1234                 return -1;
1235             }
1236             if (do_cpt_int(xd,"file_offset_low",&offset_low,list) != 0)
1237             {
1238                 return -1;
1239             }
1240         }
1241         if (file_version >= 8)
1242         {
1243             if (do_cpt_int(xd,"file_checksum_size",&(outputfiles[i].chksum_size),
1244                            list) != 0)
1245             {
1246                 return -1;
1247             }
1248             if (do_cpt_u_chars(xd,"file_checksum",16,outputfiles[i].chksum,list) != 0)
1249             {
1250                 return -1;
1251             }
1252         } 
1253         else 
1254         {
1255             outputfiles[i].chksum_size = -1;
1256         }
1257     }
1258     return 0;
1259 }
1260
1261
1262 void write_checkpoint(const char *fn,gmx_bool bNumberAndKeep,
1263                       FILE *fplog,t_commrec *cr,
1264                       int eIntegrator,int simulation_part,
1265                       gmx_bool bExpanded, int elamstats,
1266                       gmx_large_int_t step,double t,t_state *state)
1267 {
1268     t_fileio *fp;
1269     int  file_version;
1270     char *version;
1271     char *btime;
1272     char *buser;
1273     char *bhost;
1274     int  double_prec;
1275     char *fprog;
1276     char *fntemp; /* the temporary checkpoint file name */
1277     time_t now;
1278     char timebuf[STRLEN];
1279     int  nppnodes,npmenodes,flag_64bit;
1280     char buf[1024],suffix[5+STEPSTRSIZE],sbuf[STEPSTRSIZE];
1281     gmx_file_position_t *outputfiles;
1282     int  noutputfiles;
1283     char *ftime;
1284     int  flags_eks,flags_enh,flags_dfh,i;
1285     t_fileio *ret;
1286                 
1287     if (PAR(cr))
1288     {
1289         if (DOMAINDECOMP(cr))
1290         {
1291             nppnodes  = cr->dd->nnodes;
1292             npmenodes = cr->npmenodes;
1293         }
1294         else
1295         {
1296             nppnodes  = cr->nnodes;
1297             npmenodes = 0;
1298         }
1299     }
1300     else
1301     {
1302         nppnodes  = 1;
1303         npmenodes = 0;
1304     }
1305
1306     /* make the new temporary filename */
1307     snew(fntemp, strlen(fn)+5+STEPSTRSIZE);
1308     strcpy(fntemp,fn);
1309     fntemp[strlen(fn) - strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
1310     sprintf(suffix,"_%s%s","step",gmx_step_str(step,sbuf));
1311     strcat(fntemp,suffix);
1312     strcat(fntemp,fn+strlen(fn) - strlen(ftp2ext(fn2ftp(fn))) - 1);
1313    
1314     time(&now);
1315     gmx_ctime_r(&now,timebuf,STRLEN);
1316
1317     if (fplog)
1318     { 
1319         fprintf(fplog,"Writing checkpoint, step %s at %s\n\n",
1320                 gmx_step_str(step,buf),timebuf);
1321     }
1322     
1323     /* Get offsets for open files */
1324     gmx_fio_get_output_file_positions(&outputfiles, &noutputfiles);
1325
1326     fp = gmx_fio_open(fntemp,"w");
1327         
1328     if (state->ekinstate.bUpToDate)
1329     {
1330         flags_eks =
1331             ((1<<eeksEKIN_N) | (1<<eeksEKINH) | (1<<eeksEKINF) | 
1332              (1<<eeksEKINO) | (1<<eeksEKINSCALEF) | (1<<eeksEKINSCALEH) | 
1333              (1<<eeksVSCALE) | (1<<eeksDEKINDL) | (1<<eeksMVCOS));
1334     }
1335     else
1336     {
1337         flags_eks = 0;
1338     }
1339
1340     flags_enh = 0;
1341     if (state->enerhist.nsum > 0 || state->enerhist.nsum_sim > 0)
1342     {
1343         flags_enh |= (1<<eenhENERGY_N);
1344         if (state->enerhist.nsum > 0)
1345         {
1346             flags_enh |= ((1<<eenhENERGY_AVER) | (1<<eenhENERGY_SUM) |
1347                           (1<<eenhENERGY_NSTEPS) | (1<<eenhENERGY_NSUM));
1348         }
1349         if (state->enerhist.nsum_sim > 0)
1350         {
1351             flags_enh |= ((1<<eenhENERGY_SUM_SIM) | (1<<eenhENERGY_NSTEPS_SIM) |
1352                           (1<<eenhENERGY_NSUM_SIM));
1353         }
1354         if (state->enerhist.dht)
1355         {
1356             flags_enh |= ( (1<< eenhENERGY_DELTA_H_NN) |
1357                            (1<< eenhENERGY_DELTA_H_LIST) | 
1358                            (1<< eenhENERGY_DELTA_H_STARTTIME) |
1359                            (1<< eenhENERGY_DELTA_H_STARTLAMBDA) );
1360         }
1361     }
1362
1363     if (bExpanded)
1364     {
1365         flags_dfh = ((1<<edfhBEQUIL) | (1<<edfhNATLAMBDA) | (1<<edfhSUMWEIGHTS) |  (1<<edfhSUMDG)  |
1366                      (1<<edfhTIJ) | (1<<edfhTIJEMP));
1367         if (EWL(elamstats))
1368         {
1369             flags_dfh |= ((1<<edfhWLDELTA) | (1<<edfhWLHISTO));
1370         }
1371         if ((elamstats == elamstatsMINVAR) || (elamstats == elamstatsBARKER) || (elamstats == elamstatsMETROPOLIS))
1372         {
1373             flags_dfh |= ((1<<edfhACCUMP) | (1<<edfhACCUMM) | (1<<edfhACCUMP2) | (1<<edfhACCUMM2)
1374                           | (1<<edfhSUMMINVAR) | (1<<edfhSUMVAR));
1375         }
1376     } else {
1377         flags_dfh = 0;
1378     }
1379     
1380     /* We can check many more things now (CPU, acceleration, etc), but
1381      * it is highly unlikely to have two separate builds with exactly
1382      * the same version, user, time, and build host!
1383      */
1384
1385     version = gmx_strdup(VERSION);
1386     btime   = gmx_strdup(BUILD_TIME);
1387     buser   = gmx_strdup(BUILD_USER);
1388     bhost   = gmx_strdup(BUILD_HOST);
1389
1390     double_prec = GMX_CPT_BUILD_DP;
1391     fprog   = gmx_strdup(Program());
1392
1393     ftime   = &(timebuf[0]);
1394     
1395     do_cpt_header(gmx_fio_getxdr(fp),FALSE,&file_version,
1396                   &version,&btime,&buser,&bhost,&double_prec,&fprog,&ftime,
1397                   &eIntegrator,&simulation_part,&step,&t,&nppnodes,
1398                   DOMAINDECOMP(cr) ? cr->dd->nc : NULL,&npmenodes,
1399                   &state->natoms,&state->ngtc,&state->nnhpres,
1400                   &state->nhchainlength,&(state->dfhist.nlambda),&state->flags,&flags_eks,&flags_enh,&flags_dfh,
1401                   NULL);
1402     
1403     sfree(version);
1404     sfree(btime);
1405     sfree(buser);
1406     sfree(bhost);
1407     sfree(fprog);
1408
1409     if((do_cpt_state(gmx_fio_getxdr(fp),FALSE,state->flags,state,TRUE,NULL) < 0)        ||
1410        (do_cpt_ekinstate(gmx_fio_getxdr(fp),FALSE,flags_eks,&state->ekinstate,NULL) < 0)||
1411        (do_cpt_enerhist(gmx_fio_getxdr(fp),FALSE,flags_enh,&state->enerhist,NULL) < 0)  ||
1412        (do_cpt_df_hist(gmx_fio_getxdr(fp),FALSE,flags_dfh,&state->dfhist,NULL) < 0)  ||
1413        (do_cpt_files(gmx_fio_getxdr(fp),FALSE,&outputfiles,&noutputfiles,NULL,
1414                      file_version) < 0))
1415     {
1416         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
1417     }
1418
1419     do_cpt_footer(gmx_fio_getxdr(fp),FALSE,file_version);
1420
1421     /* we really, REALLY, want to make sure to physically write the checkpoint, 
1422        and all the files it depends on, out to disk. Because we've
1423        opened the checkpoint with gmx_fio_open(), it's in our list
1424        of open files.  */
1425     ret=gmx_fio_all_output_fsync();
1426
1427     if (ret)
1428     {
1429         char buf[STRLEN];
1430         sprintf(buf,
1431                 "Cannot fsync '%s'; maybe you are out of disk space?",
1432                 gmx_fio_getname(ret));
1433
1434         if (getenv(GMX_IGNORE_FSYNC_FAILURE_ENV)==NULL)
1435         {
1436             gmx_file(buf);
1437         }
1438         else
1439         {
1440             gmx_warning(buf);
1441         }
1442     }
1443
1444     if( gmx_fio_close(fp) != 0)
1445     {
1446         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
1447     }
1448
1449     /* we don't move the checkpoint if the user specified they didn't want it,
1450        or if the fsyncs failed */
1451     if (!bNumberAndKeep && !ret)
1452     {
1453         if (gmx_fexist(fn))
1454         {
1455             /* Rename the previous checkpoint file */
1456             strcpy(buf,fn);
1457             buf[strlen(fn) - strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
1458             strcat(buf,"_prev");
1459             strcat(buf,fn+strlen(fn) - strlen(ftp2ext(fn2ftp(fn))) - 1);
1460 #ifndef GMX_FAHCORE
1461             /* we copy here so that if something goes wrong between now and
1462              * the rename below, there's always a state.cpt.
1463              * If renames are atomic (such as in POSIX systems),
1464              * this copying should be unneccesary.
1465              */
1466             gmx_file_copy(fn, buf, FALSE);
1467             /* We don't really care if this fails: 
1468              * there's already a new checkpoint.
1469              */
1470 #else
1471             gmx_file_rename(fn, buf);
1472 #endif
1473         }
1474         if (gmx_file_rename(fntemp, fn) != 0)
1475         {
1476             gmx_file("Cannot rename checkpoint file; maybe you are out of disk space?");
1477         }
1478     }
1479
1480     sfree(outputfiles);
1481     sfree(fntemp);
1482
1483 #ifdef GMX_FAHCORE
1484     /*code for alternate checkpointing scheme.  moved from top of loop over 
1485       steps */
1486     fcRequestCheckPoint();
1487     if ( fcCheckPointParallel( cr->nodeid, NULL,0) == 0 ) {
1488         gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", step );
1489     }
1490 #endif /* end GMX_FAHCORE block */
1491 }
1492
1493 static void print_flag_mismatch(FILE *fplog,int sflags,int fflags)
1494 {
1495     int i;
1496     
1497     fprintf(fplog,"\nState entry mismatch between the simulation and the checkpoint file\n");
1498     fprintf(fplog,"Entries which are not present in the checkpoint file will not be updated\n");
1499     fprintf(fplog,"  %24s    %11s    %11s\n","","simulation","checkpoint");
1500     for(i=0; i<estNR; i++)
1501     {
1502         if ((sflags & (1<<i)) || (fflags & (1<<i)))
1503         {
1504             fprintf(fplog,"  %24s    %11s    %11s\n",
1505                     est_names[i],
1506                     (sflags & (1<<i)) ? "  present  " : "not present",
1507                     (fflags & (1<<i)) ? "  present  " : "not present");
1508         }
1509     }
1510 }
1511
1512 static void check_int(FILE *fplog,const char *type,int p,int f,gmx_bool *mm)
1513 {
1514         FILE *fp = fplog ? fplog : stderr;
1515
1516     if (p != f)
1517     {
1518                 fprintf(fp,"  %s mismatch,\n",type);
1519                 fprintf(fp,"    current program: %d\n",p);
1520                 fprintf(fp,"    checkpoint file: %d\n",f);
1521                 fprintf(fp,"\n");
1522         *mm = TRUE;
1523     }
1524 }
1525
1526 static void check_string(FILE *fplog,const char *type,const char *p,
1527                          const char *f,gmx_bool *mm)
1528 {
1529         FILE *fp = fplog ? fplog : stderr;
1530         
1531     if (strcmp(p,f) != 0)
1532     {
1533                 fprintf(fp,"  %s mismatch,\n",type);
1534                 fprintf(fp,"    current program: %s\n",p);
1535                 fprintf(fp,"    checkpoint file: %s\n",f);
1536                 fprintf(fp,"\n");
1537         *mm = TRUE;
1538     }
1539 }
1540
1541 static void check_match(FILE *fplog,
1542                         char *version,
1543                         char *btime,char *buser,char *bhost,int double_prec,
1544                         char *fprog,
1545                         t_commrec *cr,gmx_bool bPartDecomp,int npp_f,int npme_f,
1546                         ivec dd_nc,ivec dd_nc_f)
1547 {
1548     int  npp;
1549     gmx_bool mm;
1550     
1551     mm = FALSE;
1552     
1553     check_string(fplog,"Version"      ,VERSION      ,version,&mm);
1554     check_string(fplog,"Build time"   ,BUILD_TIME   ,btime  ,&mm);
1555     check_string(fplog,"Build user"   ,BUILD_USER   ,buser  ,&mm);
1556     check_string(fplog,"Build host"   ,BUILD_HOST   ,bhost  ,&mm);
1557     check_int   (fplog,"Double prec." ,GMX_CPT_BUILD_DP,double_prec,&mm);
1558     check_string(fplog,"Program name" ,Program()    ,fprog  ,&mm);
1559     
1560     check_int   (fplog,"#nodes"       ,cr->nnodes   ,npp_f+npme_f ,&mm);
1561     if (bPartDecomp)
1562     {
1563         dd_nc[XX] = 1;
1564         dd_nc[YY] = 1;
1565         dd_nc[ZZ] = 1;
1566     }
1567     if (cr->nnodes > 1)
1568     {
1569         check_int (fplog,"#PME-nodes"  ,cr->npmenodes,npme_f     ,&mm);
1570
1571         npp = cr->nnodes;
1572         if (cr->npmenodes >= 0)
1573         {
1574             npp -= cr->npmenodes;
1575         }
1576         if (npp == npp_f)
1577         {
1578             check_int (fplog,"#DD-cells[x]",dd_nc[XX]    ,dd_nc_f[XX],&mm);
1579             check_int (fplog,"#DD-cells[y]",dd_nc[YY]    ,dd_nc_f[YY],&mm);
1580             check_int (fplog,"#DD-cells[z]",dd_nc[ZZ]    ,dd_nc_f[ZZ],&mm);
1581         }
1582     }
1583     
1584     if (mm)
1585     {
1586                 fprintf(stderr,
1587                                 "Gromacs binary or parallel settings not identical to previous run.\n"
1588                                 "Continuation is exact, but is not guaranteed to be binary identical%s.\n\n",
1589                                 fplog ? ",\n see the log file for details" : "");
1590                 
1591         if (fplog)
1592         {
1593                         fprintf(fplog,
1594                                         "Gromacs binary or parallel settings not identical to previous run.\n"
1595                                         "Continuation is exact, but is not guaranteed to be binary identical.\n\n");
1596                 }
1597     }
1598 }
1599
1600 static void read_checkpoint(const char *fn,FILE **pfplog,
1601                             t_commrec *cr,gmx_bool bPartDecomp,ivec dd_nc,
1602                             int eIntegrator, int *init_fep_state, gmx_large_int_t *step,double *t,
1603                             t_state *state,gmx_bool *bReadRNG,gmx_bool *bReadEkin,
1604                             int *simulation_part,
1605                             gmx_bool bAppendOutputFiles,gmx_bool bForceAppend)
1606 {
1607     t_fileio *fp;
1608     int  i,j,rc;
1609     int  file_version;
1610     char *version,*btime,*buser,*bhost,*fprog,*ftime;
1611     int  double_prec;
1612         char filename[STRLEN],buf[STEPSTRSIZE];
1613     int  nppnodes,eIntegrator_f,nppnodes_f,npmenodes_f;
1614     ivec dd_nc_f;
1615     int  natoms,ngtc,nnhpres,nhchainlength,nlambda,fflags,flags_eks,flags_enh,flags_dfh;
1616     int  d;
1617     int  ret;
1618     gmx_file_position_t *outputfiles;
1619     int  nfiles;
1620     t_fileio *chksum_file;
1621     FILE* fplog = *pfplog;
1622     unsigned char digest[16];
1623 #ifndef GMX_NATIVE_WINDOWS
1624     struct flock fl;  /* don't initialize here: the struct order is OS 
1625                          dependent! */
1626 #endif
1627
1628     const char *int_warn=
1629               "WARNING: The checkpoint file was generated with integrator %s,\n"
1630               "         while the simulation uses integrator %s\n\n";
1631     const char *sd_note=
1632         "NOTE: The checkpoint file was for %d nodes doing SD or BD,\n"
1633         "      while the simulation uses %d SD or BD nodes,\n"
1634         "      continuation will be exact, except for the random state\n\n";
1635     
1636 #ifndef GMX_NATIVE_WINDOWS
1637     fl.l_type=F_WRLCK;
1638     fl.l_whence=SEEK_SET;
1639     fl.l_start=0;
1640     fl.l_len=0;
1641     fl.l_pid=0;
1642 #endif
1643
1644     if (PARTDECOMP(cr))
1645     {
1646         gmx_fatal(FARGS,
1647                   "read_checkpoint not (yet) supported with particle decomposition");
1648     }
1649     
1650     fp = gmx_fio_open(fn,"r");
1651     do_cpt_header(gmx_fio_getxdr(fp),TRUE,&file_version,
1652                   &version,&btime,&buser,&bhost,&double_prec,&fprog,&ftime,
1653                   &eIntegrator_f,simulation_part,step,t,
1654                   &nppnodes_f,dd_nc_f,&npmenodes_f,
1655                   &natoms,&ngtc,&nnhpres,&nhchainlength,&nlambda,
1656                   &fflags,&flags_eks,&flags_enh,&flags_dfh,NULL);
1657
1658     if (bAppendOutputFiles &&
1659         file_version >= 13 && double_prec != GMX_CPT_BUILD_DP)
1660     {
1661         gmx_fatal(FARGS,"Output file appending requested, but the code and checkpoint file precision (single/double) don't match");
1662     }
1663     
1664     if (cr == NULL || MASTER(cr))
1665     {
1666         fprintf(stderr,"\nReading checkpoint file %s generated: %s\n\n",
1667                 fn,ftime);
1668     }
1669         
1670         /* This will not be written if we do appending, since fplog is still NULL then */
1671     if (fplog)
1672     {
1673         fprintf(fplog,"\n");
1674         fprintf(fplog,"Reading checkpoint file %s\n",fn);
1675         fprintf(fplog,"  file generated by:     %s\n",fprog);  
1676         fprintf(fplog,"  file generated at:     %s\n",ftime);  
1677         fprintf(fplog,"  GROMACS build time:    %s\n",btime);  
1678         fprintf(fplog,"  GROMACS build user:    %s\n",buser);  
1679         fprintf(fplog,"  GROMACS build host:    %s\n",bhost);
1680         fprintf(fplog,"  GROMACS double prec.:  %d\n",double_prec);
1681         fprintf(fplog,"  simulation part #:     %d\n",*simulation_part);
1682         fprintf(fplog,"  step:                  %s\n",gmx_step_str(*step,buf));
1683         fprintf(fplog,"  time:                  %f\n",*t);  
1684         fprintf(fplog,"\n");
1685     }
1686     
1687     if (natoms != state->natoms)
1688     {
1689         gmx_fatal(FARGS,"Checkpoint file is for a system of %d atoms, while the current system consists of %d atoms",natoms,state->natoms);
1690     }
1691     if (ngtc != state->ngtc)
1692     {
1693         gmx_fatal(FARGS,"Checkpoint file is for a system of %d T-coupling groups, while the current system consists of %d T-coupling groups",ngtc,state->ngtc);
1694     }
1695     if (nnhpres != state->nnhpres)
1696     {
1697         gmx_fatal(FARGS,"Checkpoint file is for a system of %d NH-pressure-coupling variables, while the current system consists of %d NH-pressure-coupling variables",nnhpres,state->nnhpres);
1698     }
1699
1700     if (nlambda != state->dfhist.nlambda)
1701     {
1702         gmx_fatal(FARGS,"Checkpoint file is for a system with %d lambda states, while the current system consists of %d lambda states",nlambda,state->dfhist.nlambda);
1703     }
1704
1705     init_gtc_state(state,state->ngtc,state->nnhpres,nhchainlength); /* need to keep this here to keep the tpr format working */
1706     /* write over whatever was read; we use the number of Nose-Hoover chains from the checkpoint */
1707     
1708     if (eIntegrator_f != eIntegrator)
1709     {
1710         if (MASTER(cr))
1711         {
1712             fprintf(stderr,int_warn,EI(eIntegrator_f),EI(eIntegrator));
1713         }
1714                 if(bAppendOutputFiles)
1715                 {
1716                         gmx_fatal(FARGS,
1717                                           "Output file appending requested, but input/checkpoint integrators do not match.\n"
1718                                           "Stopping the run to prevent you from ruining all your data...\n"
1719                                           "If you _really_ know what you are doing, try with the -noappend option.\n");
1720                 }
1721         if (fplog)
1722         {
1723             fprintf(fplog,int_warn,EI(eIntegrator_f),EI(eIntegrator));
1724         }
1725     }
1726
1727     if (!PAR(cr))
1728     {
1729         nppnodes = 1;
1730         cr->npmenodes = 0;
1731     }
1732     else if (bPartDecomp)
1733     {
1734         nppnodes = cr->nnodes;
1735         cr->npmenodes = 0;
1736     }
1737     else if (cr->nnodes == nppnodes_f + npmenodes_f)
1738     {
1739         if (cr->npmenodes < 0)
1740         {
1741             cr->npmenodes = npmenodes_f;
1742         }
1743         nppnodes = cr->nnodes - cr->npmenodes;
1744         if (nppnodes == nppnodes_f)
1745         {
1746             for(d=0; d<DIM; d++)
1747             {
1748                 if (dd_nc[d] == 0)
1749                 {
1750                     dd_nc[d] = dd_nc_f[d];
1751                 }
1752             }
1753         }
1754     }
1755     else
1756     {
1757         /* The number of PP nodes has not been set yet */
1758         nppnodes = -1;
1759     }
1760
1761     if ((EI_SD(eIntegrator) || eIntegrator == eiBD) && nppnodes > 0)
1762     {
1763         /* Correct the RNG state size for the number of PP nodes.
1764          * Such assignments should all be moved to one central function.
1765          */
1766         state->nrng  = nppnodes*gmx_rng_n();
1767         state->nrngi = nppnodes;
1768     }
1769     
1770     *bReadRNG = TRUE;
1771     if (fflags != state->flags)
1772     {
1773                 
1774         if (MASTER(cr))
1775         {
1776                         if(bAppendOutputFiles)
1777                         {
1778                                 gmx_fatal(FARGS,
1779                                                   "Output file appending requested, but input and checkpoint states are not identical.\n"
1780                                                   "Stopping the run to prevent you from ruining all your data...\n"
1781                                                   "You can try with the -noappend option, and get more info in the log file.\n");
1782                         }
1783                         
1784             if (getenv("GMX_ALLOW_CPT_MISMATCH") == NULL)
1785             {
1786                 gmx_fatal(FARGS,"You seem to have switched ensemble, integrator, T and/or P-coupling algorithm between the cpt and tpr file. The recommended way of doing this is passing the cpt file to grompp (with option -t) instead of to mdrun. If you know what you are doing, you can override this error by setting the env.var. GMX_ALLOW_CPT_MISMATCH");
1787             }
1788             else
1789             {
1790                 fprintf(stderr,
1791                         "WARNING: The checkpoint state entries do not match the simulation,\n"
1792                         "         see the log file for details\n\n");
1793             }
1794         }
1795                 
1796                 if(fplog)
1797                 {
1798                         print_flag_mismatch(fplog,state->flags,fflags);
1799                 }
1800     }
1801     else
1802     {
1803         if ((EI_SD(eIntegrator) || eIntegrator == eiBD) &&
1804             nppnodes != nppnodes_f)
1805         {
1806             *bReadRNG = FALSE;
1807             if (MASTER(cr))
1808             {
1809                 fprintf(stderr,sd_note,nppnodes_f,nppnodes);
1810             }
1811             if (fplog)
1812             {
1813                 fprintf(fplog ,sd_note,nppnodes_f,nppnodes);
1814             }
1815         }
1816         if (MASTER(cr))
1817         {
1818             check_match(fplog,version,btime,buser,bhost,double_prec,fprog,
1819                         cr,bPartDecomp,nppnodes_f,npmenodes_f,dd_nc,dd_nc_f);
1820         }
1821     }
1822     ret = do_cpt_state(gmx_fio_getxdr(fp),TRUE,fflags,state,*bReadRNG,NULL);
1823     *init_fep_state = state->fep_state;  /* there should be a better way to do this than setting it here.
1824                                             Investigate for 5.0. */
1825     if (ret)
1826     {
1827         cp_error();
1828     }
1829     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp),TRUE,
1830                            flags_eks,&state->ekinstate,NULL);
1831     if (ret)
1832     {
1833         cp_error();
1834     }
1835     *bReadEkin = ((flags_eks & (1<<eeksEKINH)) || (flags_eks & (1<<eeksEKINF)) || (flags_eks & (1<<eeksEKINO)) ||
1836                   ((flags_eks & (1<<eeksEKINSCALEF)) | (flags_eks & (1<<eeksEKINSCALEH)) | (flags_eks & (1<<eeksVSCALE))));
1837     
1838     ret = do_cpt_enerhist(gmx_fio_getxdr(fp),TRUE,
1839                           flags_enh,&state->enerhist,NULL);
1840     if (ret)
1841     {
1842         cp_error();
1843     }
1844
1845     if (file_version < 6)
1846     {
1847         const char *warn="Reading checkpoint file in old format, assuming that the run that generated this file started at step 0, if this is not the case the averages stored in the energy file will be incorrect.";
1848
1849         fprintf(stderr,"\nWARNING: %s\n\n",warn);
1850         if (fplog)
1851         {
1852             fprintf(fplog,"\nWARNING: %s\n\n",warn);
1853         }
1854         state->enerhist.nsum     = *step;
1855         state->enerhist.nsum_sim = *step;
1856     }
1857
1858     ret = do_cpt_df_hist(gmx_fio_getxdr(fp),TRUE,
1859                          flags_dfh,&state->dfhist,NULL);
1860     if (ret)
1861     {
1862         cp_error();
1863     }
1864
1865         ret = do_cpt_files(gmx_fio_getxdr(fp),TRUE,&outputfiles,&nfiles,NULL,file_version);
1866         if (ret)
1867         {
1868                 cp_error();
1869         }
1870                                            
1871     ret = do_cpt_footer(gmx_fio_getxdr(fp),TRUE,file_version);
1872     if (ret)
1873     {
1874         cp_error();
1875     }
1876     if( gmx_fio_close(fp) != 0)
1877         {
1878         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
1879         }
1880     
1881     sfree(fprog);
1882     sfree(ftime);
1883     sfree(btime);
1884     sfree(buser);
1885     sfree(bhost);
1886         
1887         /* If the user wants to append to output files,
1888      * we use the file pointer positions of the output files stored
1889      * in the checkpoint file and truncate the files such that any frames
1890      * written after the checkpoint time are removed.
1891      * All files are md5sum checked such that we can be sure that
1892      * we do not truncate other (maybe imprortant) files.
1893          */
1894     if (bAppendOutputFiles)
1895     {
1896         if (fn2ftp(outputfiles[0].filename)!=efLOG)
1897         {
1898             /* make sure first file is log file so that it is OK to use it for 
1899              * locking
1900              */
1901             gmx_fatal(FARGS,"The first output file should always be the log "
1902                       "file but instead is: %s. Cannot do appending because of this condition.", outputfiles[0].filename);
1903         }
1904         for(i=0;i<nfiles;i++)
1905         {
1906             if (outputfiles[i].offset < 0)
1907             {
1908                 gmx_fatal(FARGS,"The original run wrote a file called '%s' which "
1909                     "is larger than 2 GB, but mdrun did not support large file"
1910                     " offsets. Can not append. Run mdrun with -noappend",
1911                     outputfiles[i].filename);
1912             }
1913 #ifdef GMX_FAHCORE
1914             chksum_file=gmx_fio_open(outputfiles[i].filename,"a");
1915
1916 #else
1917             chksum_file=gmx_fio_open(outputfiles[i].filename,"r+");
1918
1919             /* lock log file */                
1920             if (i==0)
1921             {
1922                 /* Note that there are systems where the lock operation
1923                  * will succeed, but a second process can also lock the file.
1924                  * We should probably try to detect this.
1925                  */
1926 #ifndef GMX_NATIVE_WINDOWS
1927                 if (fcntl(fileno(gmx_fio_getfp(chksum_file)), F_SETLK, &fl)
1928                     ==-1)
1929 #else
1930                 if (_locking(fileno(gmx_fio_getfp(chksum_file)), _LK_NBLCK, LONG_MAX)==-1)
1931 #endif
1932                 {
1933                     if (errno == ENOSYS)
1934                     {
1935                         if (!bForceAppend)
1936                         {
1937                             gmx_fatal(FARGS,"File locking is not supported on this system. Use -noappend or specify -append explicitly to append anyhow.");
1938                         }
1939                         else
1940                         {
1941                             fprintf(stderr,"\nNOTE: File locking is not supported on this system, will not lock %s\n\n",outputfiles[i].filename);
1942                             if (fplog)
1943                             {
1944                                 fprintf(fplog,"\nNOTE: File locking not supported on this system, will not lock %s\n\n",outputfiles[i].filename);
1945                             }
1946                         }
1947                     }
1948                     else if (errno == EACCES || errno == EAGAIN)
1949                     {
1950                         gmx_fatal(FARGS,"Failed to lock: %s. Already running "
1951                                   "simulation?", outputfiles[i].filename);
1952                     }
1953                     else
1954                     {
1955                         gmx_fatal(FARGS,"Failed to lock: %s. %s.",
1956                                   outputfiles[i].filename, strerror(errno));
1957                     }
1958                 }
1959             }
1960             
1961             /* compute md5 chksum */ 
1962             if (outputfiles[i].chksum_size != -1)
1963             {
1964                 if (gmx_fio_get_file_md5(chksum_file,outputfiles[i].offset,
1965                                      digest) != outputfiles[i].chksum_size)  /*at the end of the call the file position is at the end of the file*/
1966                 {
1967                     gmx_fatal(FARGS,"Can't read %d bytes of '%s' to compute checksum. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
1968                               outputfiles[i].chksum_size, 
1969                               outputfiles[i].filename);
1970                 }
1971             } 
1972             if (i==0)  /*log file needs to be seeked in case we need to truncate (other files are truncated below)*/
1973             {
1974                 if (gmx_fio_seek(chksum_file,outputfiles[i].offset))
1975                 {
1976                         gmx_fatal(FARGS,"Seek error! Failed to truncate log-file: %s.", strerror(errno));
1977                 }
1978             }
1979 #endif
1980             
1981             if (i==0) /*open log file here - so that lock is never lifted 
1982                         after chksum is calculated */
1983             {
1984                 *pfplog = gmx_fio_getfp(chksum_file);
1985             }
1986             else
1987             {
1988                 gmx_fio_close(chksum_file);
1989             }
1990 #ifndef GMX_FAHCORE            
1991             /* compare md5 chksum */
1992             if (outputfiles[i].chksum_size != -1 &&
1993                 memcmp(digest,outputfiles[i].chksum,16)!=0) 
1994             {
1995                 if (debug)
1996                 {
1997                     fprintf(debug,"chksum for %s: ",outputfiles[i].filename);
1998                     for (j=0; j<16; j++)
1999                     {
2000                         fprintf(debug,"%02x",digest[j]);
2001                     }
2002                     fprintf(debug,"\n");
2003                 }
2004                 gmx_fatal(FARGS,"Checksum wrong for '%s'. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
2005                           outputfiles[i].filename);
2006             }
2007 #endif        
2008
2009               
2010             if (i!=0) /*log file is already seeked to correct position */
2011             {
2012 #ifdef GMX_NATIVE_WINDOWS
2013                 rc = gmx_wintruncate(outputfiles[i].filename,outputfiles[i].offset);
2014 #else            
2015                 rc = truncate(outputfiles[i].filename,outputfiles[i].offset);
2016 #endif
2017                 if(rc!=0)
2018                 {
2019                     gmx_fatal(FARGS,"Truncation of file %s failed. Cannot do appending because of this failure.",outputfiles[i].filename);
2020                 }
2021             }
2022         }
2023     }
2024
2025     sfree(outputfiles);
2026 }
2027
2028
2029 void load_checkpoint(const char *fn,FILE **fplog,
2030                      t_commrec *cr,gmx_bool bPartDecomp,ivec dd_nc,
2031                      t_inputrec *ir,t_state *state,
2032                      gmx_bool *bReadRNG,gmx_bool *bReadEkin,
2033                      gmx_bool bAppend,gmx_bool bForceAppend)
2034 {
2035     gmx_large_int_t step;
2036     double t;
2037
2038     if (SIMMASTER(cr)) {
2039       /* Read the state from the checkpoint file */
2040       read_checkpoint(fn,fplog,
2041                       cr,bPartDecomp,dd_nc,
2042                       ir->eI,&(ir->fepvals->init_fep_state),&step,&t,state,bReadRNG,bReadEkin,
2043                       &ir->simulation_part,bAppend,bForceAppend);
2044     }
2045     if (PAR(cr)) {
2046       gmx_bcast(sizeof(cr->npmenodes),&cr->npmenodes,cr);
2047       gmx_bcast(DIM*sizeof(dd_nc[0]),dd_nc,cr);
2048       gmx_bcast(sizeof(step),&step,cr);
2049       gmx_bcast(sizeof(*bReadRNG),bReadRNG,cr);
2050       gmx_bcast(sizeof(*bReadEkin),bReadEkin,cr);
2051     }
2052     ir->bContinuation    = TRUE;
2053     if (ir->nsteps >= 0)
2054     {
2055         ir->nsteps          += ir->init_step - step;
2056     }
2057     ir->init_step        = step;
2058         ir->simulation_part += 1;
2059 }
2060
2061 static void read_checkpoint_data(t_fileio *fp,int *simulation_part,
2062                                  gmx_large_int_t *step,double *t,t_state *state,
2063                                  gmx_bool bReadRNG,
2064                                  int *nfiles,gmx_file_position_t **outputfiles)
2065 {
2066     int  file_version;
2067     char *version,*btime,*buser,*bhost,*fprog,*ftime;
2068     int  double_prec;
2069     int  eIntegrator;
2070     int  nppnodes,npme;
2071     ivec dd_nc;
2072     int  flags_eks,flags_enh,flags_dfh;
2073     int  nfiles_loc;
2074     gmx_file_position_t *files_loc=NULL;
2075     int  ret;
2076         
2077     do_cpt_header(gmx_fio_getxdr(fp),TRUE,&file_version,
2078                   &version,&btime,&buser,&bhost,&double_prec,&fprog,&ftime,
2079                   &eIntegrator,simulation_part,step,t,&nppnodes,dd_nc,&npme,
2080                   &state->natoms,&state->ngtc,&state->nnhpres,&state->nhchainlength,
2081                   &(state->dfhist.nlambda),&state->flags,&flags_eks,&flags_enh,&flags_dfh,NULL);
2082     ret =
2083         do_cpt_state(gmx_fio_getxdr(fp),TRUE,state->flags,state,bReadRNG,NULL);
2084     if (ret)
2085     {
2086         cp_error();
2087     }
2088     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp),TRUE,
2089                            flags_eks,&state->ekinstate,NULL);
2090     if (ret)
2091     {
2092         cp_error();
2093     }
2094     ret = do_cpt_enerhist(gmx_fio_getxdr(fp),TRUE,
2095                           flags_enh,&state->enerhist,NULL);
2096     if (ret)
2097     {
2098         cp_error();
2099     }
2100     ret = do_cpt_df_hist(gmx_fio_getxdr(fp),TRUE,
2101                           flags_dfh,&state->dfhist,NULL);
2102     if (ret)
2103     {
2104         cp_error();
2105     }
2106
2107     ret = do_cpt_files(gmx_fio_getxdr(fp),TRUE,
2108                        outputfiles != NULL ? outputfiles : &files_loc,
2109                        outputfiles != NULL ? nfiles : &nfiles_loc,
2110                        NULL,file_version);
2111     if (files_loc != NULL)
2112     {
2113         sfree(files_loc);
2114     }
2115         
2116     if (ret)
2117     {
2118         cp_error();
2119     }
2120         
2121     ret = do_cpt_footer(gmx_fio_getxdr(fp),TRUE,file_version);
2122     if (ret)
2123     {
2124         cp_error();
2125     }
2126
2127     sfree(fprog);
2128     sfree(ftime);
2129     sfree(btime);
2130     sfree(buser);
2131     sfree(bhost);
2132 }
2133
2134 void 
2135 read_checkpoint_state(const char *fn,int *simulation_part,
2136                       gmx_large_int_t *step,double *t,t_state *state)
2137 {
2138     t_fileio *fp;
2139     
2140     fp = gmx_fio_open(fn,"r");
2141     read_checkpoint_data(fp,simulation_part,step,t,state,FALSE,NULL,NULL);
2142     if( gmx_fio_close(fp) != 0)
2143         {
2144         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2145         }
2146 }
2147
2148 void read_checkpoint_trxframe(t_fileio *fp,t_trxframe *fr)
2149 {
2150     t_state state;
2151     int simulation_part;
2152     gmx_large_int_t step;
2153     double t;
2154     
2155     init_state(&state,0,0,0,0,0);
2156     
2157     read_checkpoint_data(fp,&simulation_part,&step,&t,&state,FALSE,NULL,NULL);
2158     
2159     fr->natoms  = state.natoms;
2160     fr->bTitle  = FALSE;
2161     fr->bStep   = TRUE;
2162     fr->step    = gmx_large_int_to_int(step,
2163                                     "conversion of checkpoint to trajectory");
2164     fr->bTime   = TRUE;
2165     fr->time    = t;
2166     fr->bLambda = TRUE;
2167     fr->lambda  = state.lambda[efptFEP];
2168     fr->fep_state  = state.fep_state;
2169     fr->bAtoms  = FALSE;
2170     fr->bX      = (state.flags & (1<<estX));
2171     if (fr->bX)
2172     {
2173         fr->x     = state.x;
2174         state.x   = NULL;
2175     }
2176     fr->bV      = (state.flags & (1<<estV));
2177     if (fr->bV)
2178     {
2179         fr->v     = state.v;
2180         state.v   = NULL;
2181     }
2182     fr->bF      = FALSE;
2183     fr->bBox    = (state.flags & (1<<estBOX));
2184     if (fr->bBox)
2185     {
2186         copy_mat(state.box,fr->box);
2187     }
2188     done_state(&state);
2189 }
2190
2191 void list_checkpoint(const char *fn,FILE *out)
2192 {
2193     t_fileio *fp;
2194     int  file_version;
2195     char *version,*btime,*buser,*bhost,*fprog,*ftime;
2196     int  double_prec;
2197     int  eIntegrator,simulation_part,nppnodes,npme;
2198     gmx_large_int_t step;
2199     double t;
2200     ivec dd_nc;
2201     t_state state;
2202     int  flags_eks,flags_enh,flags_dfh;
2203     int  indent;
2204     int  i,j;
2205     int  ret;
2206     gmx_file_position_t *outputfiles;
2207         int  nfiles;
2208         
2209     init_state(&state,-1,-1,-1,-1,0);
2210
2211     fp = gmx_fio_open(fn,"r");
2212     do_cpt_header(gmx_fio_getxdr(fp),TRUE,&file_version,
2213                   &version,&btime,&buser,&bhost,&double_prec,&fprog,&ftime,
2214                   &eIntegrator,&simulation_part,&step,&t,&nppnodes,dd_nc,&npme,
2215                   &state.natoms,&state.ngtc,&state.nnhpres,&state.nhchainlength,
2216                   &(state.dfhist.nlambda),&state.flags,
2217                   &flags_eks,&flags_enh,&flags_dfh,out);
2218     ret = do_cpt_state(gmx_fio_getxdr(fp),TRUE,state.flags,&state,TRUE,out);
2219     if (ret)
2220     {
2221         cp_error();
2222     }
2223     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp),TRUE,
2224                            flags_eks,&state.ekinstate,out);
2225     if (ret)
2226     {
2227         cp_error();
2228     }
2229     ret = do_cpt_enerhist(gmx_fio_getxdr(fp),TRUE,
2230                           flags_enh,&state.enerhist,out);
2231
2232     if (ret == 0)
2233     {
2234         init_df_history(&state.dfhist,state.dfhist.nlambda,0); /* reinitialize state with correct sizes */
2235         ret = do_cpt_df_hist(gmx_fio_getxdr(fp),TRUE,
2236                              flags_dfh,&state.dfhist,out);
2237     }
2238     if (ret == 0)
2239     {
2240                 do_cpt_files(gmx_fio_getxdr(fp),TRUE,&outputfiles,&nfiles,out,file_version);
2241         }
2242         
2243     if (ret == 0)
2244     {
2245         ret = do_cpt_footer(gmx_fio_getxdr(fp),TRUE,file_version);
2246     }
2247         
2248     if (ret)
2249     {
2250         cp_warning(out);
2251     }
2252     if( gmx_fio_close(fp) != 0)
2253         {
2254         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2255         }
2256     
2257     done_state(&state);
2258 }
2259
2260
2261 static gmx_bool exist_output_file(const char *fnm_cp,int nfile,const t_filenm fnm[])
2262 {
2263     int i;
2264
2265     /* Check if the output file name stored in the checkpoint file
2266      * is one of the output file names of mdrun.
2267      */
2268     i = 0;
2269     while (i < nfile &&
2270            !(is_output(&fnm[i]) && strcmp(fnm_cp,fnm[i].fns[0]) == 0))
2271     {
2272         i++;
2273     }
2274     
2275     return (i < nfile && gmx_fexist(fnm_cp));
2276 }
2277
2278 /* This routine cannot print tons of data, since it is called before the log file is opened. */
2279 gmx_bool read_checkpoint_simulation_part(const char *filename, int *simulation_part,
2280                                      gmx_large_int_t *cpt_step,t_commrec *cr,
2281                                      gmx_bool bAppendReq,
2282                                      int nfile,const t_filenm fnm[],
2283                                      const char *part_suffix,gmx_bool *bAddPart)
2284 {
2285     t_fileio *fp;
2286     gmx_large_int_t step=0;
2287         double t;
2288     t_state state;
2289     int  nfiles;
2290     gmx_file_position_t *outputfiles;
2291     int  nexist,f;
2292     gmx_bool bAppend;
2293     char *fn,suf_up[STRLEN];
2294
2295     bAppend = FALSE;
2296
2297     if (SIMMASTER(cr)) {
2298         if(!gmx_fexist(filename) || (!(fp = gmx_fio_open(filename,"r")) ))
2299         {
2300             *simulation_part = 0;
2301         }
2302         else 
2303         {
2304             init_state(&state,0,0,0,0,0);
2305
2306             read_checkpoint_data(fp,simulation_part,&step,&t,&state,FALSE,
2307                                  &nfiles,&outputfiles);
2308             if( gmx_fio_close(fp) != 0)
2309             {
2310                 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2311             }
2312             done_state(&state);
2313
2314             if (bAppendReq)
2315             {
2316                 nexist = 0;
2317                 for(f=0; f<nfiles; f++)
2318                 {
2319                     if (exist_output_file(outputfiles[f].filename,nfile,fnm))
2320                     {
2321                         nexist++;
2322                     }
2323                 }
2324                 if (nexist == nfiles)
2325                 {
2326                     bAppend = bAppendReq;
2327                 }
2328                 else if (nexist > 0)
2329                 {
2330                     fprintf(stderr,
2331                             "Output file appending has been requested,\n"
2332                             "but some output files listed in the checkpoint file %s\n"
2333                             "are not present or are named differently by the current program:\n",
2334                             filename);
2335                     fprintf(stderr,"output files present:");
2336                     for(f=0; f<nfiles; f++)
2337                     {
2338                         if (exist_output_file(outputfiles[f].filename,
2339                                               nfile,fnm))
2340                         {
2341                             fprintf(stderr," %s",outputfiles[f].filename);
2342                         }
2343                     }
2344                     fprintf(stderr,"\n");
2345                     fprintf(stderr,"output files not present or named differently:");
2346                     for(f=0; f<nfiles; f++)
2347                     {
2348                         if (!exist_output_file(outputfiles[f].filename,
2349                                                nfile,fnm))
2350                         {
2351                             fprintf(stderr," %s",outputfiles[f].filename);
2352                         }
2353                     }
2354                     fprintf(stderr,"\n");
2355                     
2356                     gmx_fatal(FARGS,"File appending requested, but only %d of the %d output files are present",nexist,nfiles);
2357                 }
2358             }
2359             
2360             if (bAppend)
2361             {
2362                 if (nfiles == 0)
2363                 {
2364                     gmx_fatal(FARGS,"File appending requested, but no output file information is stored in the checkpoint file");
2365                 }
2366                 fn = outputfiles[0].filename;
2367                 if (strlen(fn) < 4 ||
2368                     gmx_strcasecmp(fn+strlen(fn)-4,ftp2ext(efLOG)) == 0)
2369                 {
2370                     gmx_fatal(FARGS,"File appending requested, but the log file is not the first file listed in the checkpoint file");
2371                 }
2372                 /* Set bAddPart to whether the suffix string '.part' is present
2373                  * in the log file name.
2374                  */
2375                 strcpy(suf_up,part_suffix);
2376                 upstring(suf_up);
2377                 *bAddPart = (strstr(fn,part_suffix) != NULL ||
2378                              strstr(fn,suf_up) != NULL);
2379             }
2380
2381             sfree(outputfiles);
2382         }
2383     }
2384     if (PAR(cr))
2385     {
2386         gmx_bcast(sizeof(*simulation_part),simulation_part,cr);
2387
2388         if (*simulation_part > 0 && bAppendReq)
2389         {
2390             gmx_bcast(sizeof(bAppend),&bAppend,cr);
2391             gmx_bcast(sizeof(*bAddPart),bAddPart,cr);
2392         }
2393     }
2394     if (NULL != cpt_step)
2395     {
2396         *cpt_step = step;
2397     }
2398
2399     return bAppend;
2400 }