Reimplement constant acceleration groups
[alexxy/gromacs.git] / src / gromacs / fileio / tpxio.cpp
index 448bb9dd58aef89cf2ee7c09cbf4fc0087bf7821..d190d663ed65439cde1a942b82d33e4199ba6572 100644 (file)
@@ -49,6 +49,7 @@
 #include <memory>
 #include <vector>
 
+#include "gromacs/applied_forces/awh/read_params.h"
 #include "gromacs/fileio/filetypes.h"
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/fileio/gmxfio_xdr.h"
@@ -75,6 +76,7 @@
 #include "gromacs/utility/futil.h"
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/inmemoryserializer.h"
+#include "gromacs/utility/iserializer.h"
 #include "gromacs/utility/keyvaluetreebuilder.h"
 #include "gromacs/utility/keyvaluetreeserializer.h"
 #include "gromacs/utility/smalloc.h"
@@ -93,7 +95,7 @@
  * merging with mainstream GROMACS, set this tag string back to
  * TPX_TAG_RELEASE, and instead add an element to tpxv.
  */
-static const char* tpx_tag = TPX_TAG_RELEASE;
+static const std::string tpx_tag = TPX_TAG_RELEASE;
 
 /*! \brief Enum of values that describe the contents of a tpr file
  * whose format matches a version number
@@ -135,6 +137,9 @@ enum tpxv
     tpxv_VSite1,                                  /**< Added 1 type virtual site */
     tpxv_MTS,                                     /**< Added multiple time stepping */
     tpxv_RemovedConstantAcceleration, /**< Removed support for constant acceleration NEMD. */
+    tpxv_TransformationPullCoord,     /**< Support for transformation pull coordinates */
+    tpxv_SoftcoreGapsys,              /**< Added gapsys softcore function */
+    tpxv_ReaddedConstantAcceleration, /**< Re-added support for constant acceleration NEMD. */
     tpxv_Count                        /**< the total number of tpxv versions */
 };
 
@@ -245,9 +250,6 @@ static const t_ftupd ftupd[] = {
 };
 #define NFTUPD asize(ftupd)
 
-/* Needed for backward compatibility */
-#define MAXNODES 256
-
 /**************************************************************
  *
  * Now the higer level routines that do io of the structures and arrays
@@ -302,17 +304,7 @@ static void do_pull_coord(gmx::ISerializer* serializer,
         {
             if (pcrd->eType == PullingAlgorithm::External)
             {
-                std::string buf;
-                if (serializer->reading())
-                {
-                    serializer->doString(&buf);
-                    pcrd->externalPotentialProvider = gmx_strdup(buf.c_str());
-                }
-                else
-                {
-                    buf = pcrd->externalPotentialProvider;
-                    serializer->doString(&buf);
-                }
+                serializer->doString(&pcrd->externalPotentialProvider);
             }
             else
             {
@@ -351,6 +343,17 @@ static void do_pull_coord(gmx::ISerializer* serializer,
             pcrd->ngroup = 0;
         }
         serializer->doIvec(&pcrd->dim.as_vec());
+        if (file_version >= tpxv_TransformationPullCoord)
+        {
+            serializer->doString(&pcrd->expression);
+        }
+        else
+        {
+            if (serializer->reading())
+            {
+                pcrd->expression.clear();
+            }
+        }
     }
     else
     {
@@ -445,9 +448,9 @@ static void do_simtempvals(gmx::ISerializer* serializer, t_simtemp* simtemp, int
         {
             if (serializer->reading())
             {
-                snew(simtemp->temperatures, n_lambda);
+                simtemp->temperatures.resize(n_lambda);
             }
-            serializer->doRealArray(simtemp->temperatures, n_lambda);
+            serializer->doRealArray(simtemp->temperatures.data(), n_lambda);
         }
     }
 }
@@ -496,9 +499,7 @@ static void do_fepvals(gmx::ISerializer* serializer, t_lambda* fepvals, int file
             {
                 fepvals->all_lambda[g].resize(fepvals->n_lambda);
                 serializer->doDoubleArray(fepvals->all_lambda[g].data(), fepvals->n_lambda);
-                serializer->doBoolArray(
-                        fepvals->separate_dvdl.begin(),
-                        gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real>::size());
+                serializer->doBoolArray(fepvals->separate_dvdl.begin(), fepvals->separate_dvdl.size());
             }
             else if (fepvals->init_lambda >= 0)
             {
@@ -623,6 +624,20 @@ static void do_fepvals(gmx::ISerializer* serializer, t_lambda* fepvals, int file
     {
         fepvals->edHdLPrintEnergy = FreeEnergyPrintEnergy::No;
     }
+    if (file_version >= tpxv_SoftcoreGapsys)
+    {
+        serializer->doInt(reinterpret_cast<int*>(&fepvals->softcoreFunction));
+        serializer->doReal(&fepvals->scGapsysScaleLinpointLJ);
+        serializer->doReal(&fepvals->scGapsysScaleLinpointQ);
+        serializer->doReal(&fepvals->scGapsysSigmaLJ);
+    }
+    else
+    {
+        fepvals->softcoreFunction        = SoftcoreType::Beutler;
+        fepvals->scGapsysScaleLinpointLJ = 0.85;
+        fepvals->scGapsysScaleLinpointQ  = 0.3;
+        fepvals->scGapsysSigmaLJ         = 0.3;
+    }
 
     /* handle lambda_neighbors */
     if ((file_version >= 83 && file_version < 90) || file_version >= 92)
@@ -655,73 +670,6 @@ static void do_fepvals(gmx::ISerializer* serializer, t_lambda* fepvals, int file
     }
 }
 
-static void do_awhBias(gmx::ISerializer* serializer, gmx::AwhBiasParams* awhBiasParams, int gmx_unused file_version)
-{
-    serializer->doInt(&awhBiasParams->eTarget);
-    serializer->doDouble(&awhBiasParams->targetBetaScaling);
-    serializer->doDouble(&awhBiasParams->targetCutoff);
-    serializer->doInt(&awhBiasParams->eGrowth);
-    if (serializer->reading())
-    {
-        int temp = 0;
-        serializer->doInt(&temp);
-        awhBiasParams->bUserData = static_cast<bool>(temp);
-    }
-    else
-    {
-        int temp = static_cast<int>(awhBiasParams->bUserData);
-        serializer->doInt(&temp);
-    }
-    serializer->doDouble(&awhBiasParams->errorInitial);
-    serializer->doInt(&awhBiasParams->ndim);
-    serializer->doInt(&awhBiasParams->shareGroup);
-    serializer->doBool(&awhBiasParams->equilibrateHistogram);
-
-    if (serializer->reading())
-    {
-        snew(awhBiasParams->dimParams, awhBiasParams->ndim);
-    }
-
-    for (int d = 0; d < awhBiasParams->ndim; d++)
-    {
-        gmx::AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
-
-        serializer->doInt(&dimParams->eCoordProvider);
-        serializer->doInt(&dimParams->coordIndex);
-        serializer->doDouble(&dimParams->origin);
-        serializer->doDouble(&dimParams->end);
-        serializer->doDouble(&dimParams->period);
-        serializer->doDouble(&dimParams->forceConstant);
-        serializer->doDouble(&dimParams->diffusion);
-        serializer->doDouble(&dimParams->coordValueInit);
-        serializer->doDouble(&dimParams->coverDiameter);
-    }
-}
-
-static void do_awh(gmx::ISerializer* serializer, gmx::AwhParams* awhParams, int gmx_unused file_version)
-{
-    serializer->doInt(&awhParams->numBias);
-    serializer->doInt(&awhParams->nstOut);
-    serializer->doInt64(&awhParams->seed);
-    serializer->doInt(&awhParams->nstSampleCoord);
-    serializer->doInt(&awhParams->numSamplesUpdateFreeEnergy);
-    serializer->doInt(&awhParams->ePotential);
-    serializer->doBool(&awhParams->shareBiasMultisim);
-
-    if (awhParams->numBias > 0)
-    {
-        if (serializer->reading())
-        {
-            snew(awhParams->awhBiasParams, awhParams->numBias);
-        }
-
-        for (int k = 0; k < awhParams->numBias; k++)
-        {
-            do_awhBias(serializer, &awhParams->awhBiasParams[k], file_version);
-        }
-    }
-}
-
 static void do_pull(gmx::ISerializer* serializer, pull_params_t* pull, int file_version, PullingAlgorithm ePullOld)
 {
     PullGroupGeometry eGeomOld = PullGroupGeometry::Count;
@@ -828,6 +776,10 @@ static void do_pull(gmx::ISerializer* serializer, pull_params_t* pull, int file_
         for (g = 0; g < pull->ncoord; g++)
         {
             do_pull_coord(serializer, &pull->coord[g], file_version, ePullOld, eGeomOld, dimOld);
+            if (serializer->reading())
+            {
+                pull->coord[g].coordIndex = g;
+            }
         }
     }
     if (file_version >= tpxv_PullAverage)
@@ -926,7 +878,7 @@ static void do_swapgroup(gmx::ISerializer* serializer, t_swapGroup* g)
     serializer->doIntArray(g->ind, g->nat);
 
     /* Requested counts for compartments A and B */
-    serializer->doIntArray(g->nmolReq, eCompNR);
+    serializer->doIntArray(g->nmolReq.data(), static_cast<int>(Compartment::Count));
 }
 
 static void do_swapcoords_tpx(gmx::ISerializer* serializer, t_swapcoords* swap, int file_version)
@@ -1568,9 +1520,12 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v
         {
             if (serializer->reading())
             {
-                snew(ir->awhParams, 1);
+                ir->awhParams = std::make_unique<gmx::AwhParams>(serializer);
+            }
+            else
+            {
+                ir->awhParams->serialize(serializer);
             }
-            do_awh(serializer, ir->awhParams, file_version);
         }
     }
     else
@@ -1625,10 +1580,14 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v
     {
         ir->opts.nhchainlength = 1;
     }
-    int removedOptsNgacc = 0;
-    if (serializer->reading() && file_version < tpxv_RemovedConstantAcceleration)
+    if (serializer->reading() && file_version >= tpxv_RemovedConstantAcceleration
+        && file_version < tpxv_ReaddedConstantAcceleration)
     {
-        serializer->doInt(&removedOptsNgacc);
+        ir->opts.ngacc = 0;
+    }
+    else
+    {
+        serializer->doInt(&ir->opts.ngacc);
     }
     serializer->doInt(&ir->opts.ngfrz);
     serializer->doInt(&ir->opts.ngener);
@@ -1643,6 +1602,7 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v
         snew(ir->opts.anneal_temp, ir->opts.ngtc);
         snew(ir->opts.tau_t, ir->opts.ngtc);
         snew(ir->opts.nFreeze, ir->opts.ngfrz);
+        snew(ir->opts.acceleration, ir->opts.ngacc);
         snew(ir->opts.egp_flags, ir->opts.ngener * ir->opts.ngener);
     }
     if (ir->opts.ngtc > 0)
@@ -1655,18 +1615,20 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v
     {
         serializer->doIvecArray(ir->opts.nFreeze, ir->opts.ngfrz);
     }
-    if (serializer->reading() && file_version < tpxv_RemovedConstantAcceleration && removedOptsNgacc > 0)
+    if (ir->opts.ngacc > 0)
     {
-        std::vector<gmx::RVec> dummy;
-        dummy.resize(removedOptsNgacc);
-        serializer->doRvecArray(reinterpret_cast<rvec*>(dummy.data()), removedOptsNgacc);
-        ir->useConstantAcceleration = std::any_of(dummy.begin(), dummy.end(), [](const gmx::RVec& vec) {
-            return vec[XX] != 0.0 || vec[YY] != 0.0 || vec[ZZ] != 0.0;
-        });
+        serializer->doRvecArray(ir->opts.acceleration, ir->opts.ngacc);
     }
-    else
+    if (serializer->reading())
     {
         ir->useConstantAcceleration = false;
+        for (int g = 0; g < ir->opts.ngacc; g++)
+        {
+            if (norm2(ir->opts.acceleration[g]) != 0)
+            {
+                ir->useConstantAcceleration = true;
+            }
+        }
     }
     serializer->doIntArray(ir->opts.egp_flags, ir->opts.ngener * ir->opts.ngener);
 
@@ -2593,15 +2555,11 @@ static void do_molblock(gmx::ISerializer* serializer, gmx_molblock_t* molb, int
 
 static void set_disres_npair(gmx_mtop_t* mtop)
 {
-    gmx_mtop_ilistloop_t iloop;
-    int                  nmol;
-
     gmx::ArrayRef<t_iparams> ip = mtop->ffparams.iparams;
 
-    iloop = gmx_mtop_ilistloop_init(mtop);
-    while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
+    for (const auto ilist : IListRange(*mtop))
     {
-        const InteractionList& il = (*ilist)[F_DISRES];
+        const InteractionList& il = ilist.list()[F_DISRES];
 
         if (!il.empty())
         {
@@ -2682,11 +2640,6 @@ static void do_mtop(gmx::ISerializer* serializer, gmx_mtop_t* mtop, int file_ver
     }
 
     do_groups(serializer, &mtop->groups, &(mtop->symtab));
-    if (file_version < tpxv_RemovedConstantAcceleration)
-    {
-        mtop->groups.groups[SimulationAtomGroupType::AccelerationUnused].clear();
-        mtop->groups.groupNumbers[SimulationAtomGroupType::AccelerationUnused].clear();
-    }
 
     mtop->haveMoleculeIndices = true;
 
@@ -2783,7 +2736,7 @@ static void do_tpxheader(gmx::FileIOXdrSerializer* serializer,
         serializer->doString(&buf);
         gmx_fio_setprecision(fio, tpx->isDouble);
         serializer->doInt(&precision);
-        fileTag = gmx::formatString("%s", tpx_tag);
+        fileTag = tpx_tag;
     }
 
     /* Check versions! */
@@ -2815,7 +2768,7 @@ static void do_tpxheader(gmx::FileIOXdrSerializer* serializer,
 
         if (fileTag != tpx_tag)
         {
-            fprintf(stderr, "Note: file tpx tag '%s', software tpx tag '%s'\n", fileTag.c_str(), tpx_tag);
+            fprintf(stderr, "Note: file tpx tag '%s', software tpx tag '%s'\n", fileTag.c_str(), tpx_tag.c_str());
 
             /* We only support reading tpx files with the same tag as the code
              * or tpx files with the release tag and with lower version number.
@@ -2829,7 +2782,7 @@ static void do_tpxheader(gmx::FileIOXdrSerializer* serializer,
                           tpx->fileVersion,
                           fileTag.c_str(),
                           tpx_version,
-                          tpx_tag);
+                          tpx_tag.c_str());
             }
         }
     }
@@ -3355,7 +3308,7 @@ TpxFileHeader readTpxHeader(const char* fileName, bool canReadTopologyOnly)
     return tpx;
 }
 
-void write_tpx_state(const char* fn, const t_inputrec* ir, const t_state* state, const gmx_mtop_t* mtop)
+void write_tpx_state(const char* fn, const t_inputrec* ir, const t_state* state, const gmx_mtop_t& mtop)
 {
     /* To write a state, we first need to write the state information to a buffer before
      * we append the raw bytes to the file. For this, the header information needs to be
@@ -3365,7 +3318,7 @@ void write_tpx_state(const char* fn, const t_inputrec* ir, const t_state* state,
 
     t_fileio* fio;
 
-    TpxFileHeader tpx = populateTpxHeader(*state, ir, mtop);
+    TpxFileHeader tpx = populateTpxHeader(*state, ir, &mtop);
     // Long-term we should move to use little endian in files to avoid extra byte swapping,
     // but since we just used the default XDR format (which is big endian) for the TPR
     // header it would cause third-party libraries reading our raw data to tear their hair
@@ -3379,7 +3332,7 @@ void write_tpx_state(const char* fn, const t_inputrec* ir, const t_state* state,
                 const_cast<t_state*>(state),
                 nullptr,
                 nullptr,
-                const_cast<gmx_mtop_t*>(mtop));
+                const_cast<gmx_mtop_t*>(&mtop));
 
     std::vector<char> tprBody = tprBodySerializer.finishAndGetBuffer();
     tpx.sizeOfTprBody         = tprBody.size();