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