Allow for "computational electrophysiology" simulations (CompEl)
[alexxy/gromacs.git] / src / gromacs / gmxlib / checkpoint.c
index 65984af7025b2eed4b009f3893e9353631a01b62..901bb1234a3d49175bcc5c1006b10ad1c24c7781 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2008,2009,2010,2011,2012,2013, by the GROMACS development team, led by
+ * Copyright (c) 2008,2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -99,7 +99,7 @@
  * But old code can not read a new entry that is present in the file
  * (but can read a new format when new entries are not present).
  */
-static const int cpt_version = 15;
+static const int cpt_version = 16;
 
 
 const char *est_names[estNR] =
@@ -810,7 +810,7 @@ static void do_cpt_header(XDR *xd, gmx_bool bRead, int *file_version,
                           int *natoms, int *ngtc, int *nnhpres, int *nhchainlength,
                           int *nlambda, int *flags_state,
                           int *flags_eks, int *flags_enh, int *flags_dfh,
-                          int *nED,
+                          int *nED, int *eSwapCoords,
                           FILE *list)
 {
     bool_t res = 0;
@@ -953,6 +953,10 @@ static void do_cpt_header(XDR *xd, gmx_bool bRead, int *file_version,
     {
         *nED = 0;
     }
+    if (*file_version >= 16)
+    {
+        do_cpt_int_err(xd, "swap", eSwapCoords, list);
+    }
 }
 
 static int do_cpt_footer(XDR *xd, int file_version)
@@ -1090,6 +1094,123 @@ static int do_cpt_ekinstate(XDR *xd, int fflags, ekinstate_t *ekins,
 }
 
 
+static int do_cpt_swapstate(XDR *xd, gmx_bool bRead, swapstate_t *swapstate, FILE *list)
+{
+    int ii, ic, j;
+    int ret              = 0;
+    int swap_cpt_version = 1;
+
+
+    if (eswapNO == swapstate->eSwapCoords)
+    {
+        return ret;
+    }
+
+    swapstate->bFromCpt = bRead;
+
+    do_cpt_int_err(xd, "swap checkpoint version", &swap_cpt_version, list);
+    do_cpt_int_err(xd, "swap coupling steps", &swapstate->nAverage, list);
+
+    /* When reading, init_swapcoords has not been called yet,
+     * so we have to allocate memory first. */
+
+    for (ic = 0; ic < eCompNR; ic++)
+    {
+        for (ii = 0; ii < eIonNR; ii++)
+        {
+            if (bRead)
+            {
+                do_cpt_int_err(xd, "swap requested atoms", &swapstate->nat_req[ic][ii], list);
+            }
+            else
+            {
+                do_cpt_int_err(xd, "swap requested atoms p", swapstate->nat_req_p[ic][ii], list);
+            }
+
+            if (bRead)
+            {
+                do_cpt_int_err(xd, "swap influx netto", &swapstate->inflow_netto[ic][ii], list);
+            }
+            else
+            {
+                do_cpt_int_err(xd, "swap influx netto p", swapstate->inflow_netto_p[ic][ii], list);
+            }
+
+            if (bRead && (NULL == swapstate->nat_past[ic][ii]) )
+            {
+                snew(swapstate->nat_past[ic][ii], swapstate->nAverage);
+            }
+
+            for (j = 0; j < swapstate->nAverage; j++)
+            {
+                if (bRead)
+                {
+                    do_cpt_int_err(xd, "swap past atom counts", &swapstate->nat_past[ic][ii][j], list);
+                }
+                else
+                {
+                    do_cpt_int_err(xd, "swap past atom counts p", &swapstate->nat_past_p[ic][ii][j], list);
+                }
+            }
+        }
+    }
+
+    /* Ion flux per channel */
+    for (ic = 0; ic < eChanNR; ic++)
+    {
+        for (ii = 0; ii < eIonNR; ii++)
+        {
+            if (bRead)
+            {
+                do_cpt_int_err(xd, "channel flux", &swapstate->fluxfromAtoB[ic][ii], list);
+            }
+            else
+            {
+                do_cpt_int_err(xd, "channel flux p", swapstate->fluxfromAtoB_p[ic][ii], list);
+            }
+        }
+    }
+
+    /* Ion flux leakage */
+    if (bRead)
+    {
+        snew(swapstate->fluxleak, 1);
+    }
+    do_cpt_int_err(xd, "flux leakage", swapstate->fluxleak, list);
+
+    /* Ion history */
+    do_cpt_int_err(xd, "number of ions", &swapstate->nions, list);
+
+    if (bRead)
+    {
+        snew(swapstate->channel_label, swapstate->nions);
+        snew(swapstate->comp_from, swapstate->nions);
+    }
+
+    do_cpt_u_chars(xd, "channel history", swapstate->nions, swapstate->channel_label, list);
+    do_cpt_u_chars(xd, "domain history", swapstate->nions, swapstate->comp_from, list);
+
+    /* Save the last known whole positions to checkpoint
+     * file to be able to also make multimeric channels whole in PBC */
+    do_cpt_int_err(xd, "Ch0 atoms", &swapstate->nat[eChan0], list);
+    do_cpt_int_err(xd, "Ch1 atoms", &swapstate->nat[eChan1], list);
+    if (bRead)
+    {
+        snew(swapstate->xc_old_whole[eChan0], swapstate->nat[eChan0]);
+        snew(swapstate->xc_old_whole[eChan1], swapstate->nat[eChan1]);
+        do_cpt_n_rvecs_err(xd, "Ch0 whole x", swapstate->nat[eChan0], swapstate->xc_old_whole[eChan0], list);
+        do_cpt_n_rvecs_err(xd, "Ch1 whole x", swapstate->nat[eChan1], swapstate->xc_old_whole[eChan1], list);
+    }
+    else
+    {
+        do_cpt_n_rvecs_err(xd, "Ch0 whole x", swapstate->nat[eChan0], *swapstate->xc_old_whole_p[eChan0], list);
+        do_cpt_n_rvecs_err(xd, "Ch1 whole x", swapstate->nat[eChan1], *swapstate->xc_old_whole_p[eChan1], list);
+    }
+
+    return ret;
+}
+
+
 static int do_cpt_enerhist(XDR *xd, gmx_bool bRead,
                            int fflags, energyhistory_t *enerhist,
                            FILE *list)
@@ -1532,7 +1653,7 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
                   DOMAINDECOMP(cr) ? cr->dd->nc : NULL, &npmenodes,
                   &state->natoms, &state->ngtc, &state->nnhpres,
                   &state->nhchainlength, &(state->dfhist.nlambda), &state->flags, &flags_eks, &flags_enh, &flags_dfh,
-                  &state->edsamstate.nED,
+                  &state->edsamstate.nED, &state->swapstate.eSwapCoords,
                   NULL);
 
     sfree(version);
@@ -1546,6 +1667,7 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
         (do_cpt_enerhist(gmx_fio_getxdr(fp), FALSE, flags_enh, &state->enerhist, NULL) < 0)  ||
         (do_cpt_df_hist(gmx_fio_getxdr(fp), flags_dfh, &state->dfhist, NULL) < 0)  ||
         (do_cpt_EDstate(gmx_fio_getxdr(fp), FALSE, &state->edsamstate, NULL) < 0)      ||
+        (do_cpt_swapstate(gmx_fio_getxdr(fp), FALSE, &state->swapstate, NULL) < 0) ||
         (do_cpt_files(gmx_fio_getxdr(fp), FALSE, &outputfiles, &noutputfiles, NULL,
                       file_version) < 0))
     {
@@ -1793,7 +1915,7 @@ static void read_checkpoint(const char *fn, FILE **pfplog,
                   &nppnodes_f, dd_nc_f, &npmenodes_f,
                   &natoms, &ngtc, &nnhpres, &nhchainlength, &nlambda,
                   &fflags, &flags_eks, &flags_enh, &flags_dfh,
-                  &state->edsamstate.nED, NULL);
+                  &state->edsamstate.nED, &state->swapstate.eSwapCoords, NULL);
 
     if (bAppendOutputFiles &&
         file_version >= 13 && double_prec != GMX_CPT_BUILD_DP)
@@ -2000,6 +2122,12 @@ static void read_checkpoint(const char *fn, FILE **pfplog,
         cp_error();
     }
 
+    ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, &state->swapstate, NULL);
+    if (ret)
+    {
+        cp_error();
+    }
+
     ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, &state->edsamstate, NULL);
     if (ret)
     {
@@ -2228,7 +2356,7 @@ static void read_checkpoint_data(t_fileio *fp, int *simulation_part,
                   &eIntegrator, simulation_part, step, t, &nppnodes, dd_nc, &npme,
                   &state->natoms, &state->ngtc, &state->nnhpres, &state->nhchainlength,
                   &(state->dfhist.nlambda), &state->flags, &flags_eks, &flags_enh, &flags_dfh,
-                  &state->edsamstate.nED, NULL);
+                  &state->edsamstate.nED, &state->swapstate.eSwapCoords, NULL);
     ret =
         do_cpt_state(gmx_fio_getxdr(fp), TRUE, state->flags, state, bReadRNG, NULL);
     if (ret)
@@ -2258,6 +2386,12 @@ static void read_checkpoint_data(t_fileio *fp, int *simulation_part,
         cp_error();
     }
 
+    ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, &state->swapstate, NULL);
+    if (ret)
+    {
+        cp_error();
+    }
+
     ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE,
                        outputfiles != NULL ? outputfiles : &files_loc,
                        outputfiles != NULL ? nfiles : &nfiles_loc,
@@ -2372,7 +2506,8 @@ void list_checkpoint(const char *fn, FILE *out)
                   &eIntegrator, &simulation_part, &step, &t, &nppnodes, dd_nc, &npme,
                   &state.natoms, &state.ngtc, &state.nnhpres, &state.nhchainlength,
                   &(state.dfhist.nlambda), &state.flags,
-                  &flags_eks, &flags_enh, &flags_dfh, &state.edsamstate.nED, out);
+                  &flags_eks, &flags_enh, &flags_dfh, &state.edsamstate.nED,
+                  &state.swapstate.eSwapCoords, out);
     ret = do_cpt_state(gmx_fio_getxdr(fp), TRUE, state.flags, &state, TRUE, out);
     if (ret)
     {
@@ -2397,6 +2532,11 @@ void list_checkpoint(const char *fn, FILE *out)
         ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, &state.edsamstate, out);
     }
 
+    if (ret == 0)
+    {
+        ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, &state.swapstate, out);
+    }
+
     if (ret == 0)
     {
         do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, &nfiles, out, file_version);