src/gromacs/mdtypes/state.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #include "gmxpre.h"
  40
  41 #include "state.h"
  42
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/math/paddedvector.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/math/veccompare.h"
  50 #include "gromacs/mdtypes/awh_history.h"
  51 #include "gromacs/mdtypes/df_history.h"
  52 #include "gromacs/mdtypes/inputrec.h"
  53 #include "gromacs/mdtypes/md_enums.h"
  54 #include "gromacs/mdtypes/pull_params.h"
  55 #include "gromacs/mdtypes/swaphistory.h"
  56 #include "gromacs/pbcutil/boxutilities.h"
  57 #include "gromacs/pbcutil/pbc.h"
  58 #include "gromacs/utility/compare.h"
  59 #include "gromacs/utility/gmxassert.h"
  60 #include "gromacs/utility/smalloc.h"
  61 #include "gromacs/utility/stringutil.h"
  62
  63 #include "checkpointdata.h"
  64
  65 /* The source code in this file should be thread-safe.
  66       Please keep it that way. */
  67
  68 history_t::history_t() : disre_initf(0), orire_initf(0) {}
  69
  70 ekinstate_t::ekinstate_t() :
  71     bUpToDate(FALSE),
  72     ekin_n(0),
  73     ekinh(nullptr),
  74     ekinf(nullptr),
  75     ekinh_old(nullptr),
  76     ekin_total(),
  77
  78     dekindl(0),
  79     mvcos(0)
  80 {
  81     clear_mat(ekin_total);
  82 }
  83
  84 namespace
  85 {
  86 /*!
  87  * \brief Enum describing the contents ekinstate_t writes to modular checkpoint
  88  *
  89  * When changing the checkpoint content, add a new element just above Count, and adjust the
  90  * checkpoint functionality.
  91  */
  92 enum class CheckpointVersion
  93 {
  94     Base, //!< First version of modular checkpointing
  95     Count //!< Number of entries. Add new versions right above this!
  96 };
  97 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
  98 } // namespace
  99
 100 template<gmx::CheckpointDataOperation operation>
 101 void ekinstate_t::doCheckpoint(gmx::CheckpointData<operation> checkpointData)
 102 {
 103     gmx::checkpointVersion(&checkpointData, "ekinstate_t version", c_currentVersion);
 104
 105     checkpointData.scalar("bUpToDate", &bUpToDate);
 106     if (!bUpToDate)
 107     {
 108         return;
 109     }
 110     auto numOfTensors = ekin_n;
 111     checkpointData.scalar("ekin_n", &numOfTensors);
 112     if (operation == gmx::CheckpointDataOperation::Read)
 113     {
 114         // If this isn't matching, we haven't allocated the right amount of data
 115         GMX_RELEASE_ASSERT(numOfTensors == ekin_n,
 116                            "ekinstate_t checkpoint reading: Tensor size mismatch.");
 117     }
 118     for (int idx = 0; idx < numOfTensors; ++idx)
 119     {
 120         checkpointData.tensor(gmx::formatString("ekinh %d", idx), ekinh[idx]);
 121         checkpointData.tensor(gmx::formatString("ekinf %d", idx), ekinf[idx]);
 122         checkpointData.tensor(gmx::formatString("ekinh_old %d", idx), ekinh_old[idx]);
 123     }
 124     checkpointData.arrayRef("ekinscalef_nhc", gmx::makeCheckpointArrayRef<operation>(ekinscalef_nhc));
 125     checkpointData.arrayRef("ekinscaleh_nhc", gmx::makeCheckpointArrayRef<operation>(ekinscaleh_nhc));
 126     checkpointData.arrayRef("vscale_nhc", gmx::makeCheckpointArrayRef<operation>(vscale_nhc));
 127     checkpointData.scalar("dekindl", &dekindl);
 128     checkpointData.scalar("mvcos", &mvcos);
 129
 130     if (operation == gmx::CheckpointDataOperation::Read)
 131     {
 132         hasReadEkinState = true;
 133     }
 134 }
 135
 136 // Explicit template instantiation
 137 template void ekinstate_t::doCheckpoint(gmx::CheckpointData<gmx::CheckpointDataOperation::Read> checkpointData);
 138 template void ekinstate_t::doCheckpoint(gmx::CheckpointData<gmx::CheckpointDataOperation::Write> checkpointData);
 139
 140 void init_gtc_state(t_state* state, int ngtc, int nnhpres, int nhchainlength)
 141 {
 142     state->ngtc          = ngtc;
 143     state->nnhpres       = nnhpres;
 144     state->nhchainlength = nhchainlength;
 145     state->nosehoover_xi.resize(state->nhchainlength * state->ngtc, 0);
 146     state->nosehoover_vxi.resize(state->nhchainlength * state->ngtc, 0);
 147     state->therm_integral.resize(state->ngtc, 0);
 148     state->baros_integral = 0.0;
 149     state->nhpres_xi.resize(state->nhchainlength * nnhpres, 0);
 150     state->nhpres_vxi.resize(state->nhchainlength * nnhpres, 0);
 151 }
 152
 153
 154 /* Checkpoint code relies on this function having no effect if
 155    state->natoms is > 0 and passed as natoms. */
 156 void state_change_natoms(t_state* state, int natoms)
 157 {
 158     state->natoms = natoms;
 159
 160     /* We need padding, since we might use SIMD access, but the
 161      * containers here all ensure that. */
 162     if (state->flags & enumValueToBitMask(StateEntry::X))
 163     {
 164         state->x.resizeWithPadding(natoms);
 165     }
 166     if (state->flags & enumValueToBitMask(StateEntry::V))
 167     {
 168         state->v.resizeWithPadding(natoms);
 169     }
 170     if (state->flags & enumValueToBitMask(StateEntry::Cgp))
 171     {
 172         state->cg_p.resizeWithPadding(natoms);
 173     }
 174 }
 175
 176 void init_dfhist_state(t_state* state, int dfhistNumLambda)
 177 {
 178     if (dfhistNumLambda > 0)
 179     {
 180         snew(state->dfhist, 1);
 181         init_df_history(state->dfhist, dfhistNumLambda);
 182     }
 183     else
 184     {
 185         state->dfhist = nullptr;
 186     }
 187 }
 188
 189 void comp_state(const t_state* st1, const t_state* st2, gmx_bool bRMSD, real ftol, real abstol)
 190 {
 191     int i, j, nc;
 192
 193     fprintf(stdout, "comparing flags\n");
 194     cmp_int(stdout, "flags", -1, st1->flags, st2->flags);
 195     fprintf(stdout, "comparing box\n");
 196     cmp_rvecs(stdout, "box", DIM, st1->box, st2->box, FALSE, ftol, abstol);
 197     fprintf(stdout, "comparing box_rel\n");
 198     cmp_rvecs(stdout, "box_rel", DIM, st1->box_rel, st2->box_rel, FALSE, ftol, abstol);
 199     fprintf(stdout, "comparing boxv\n");
 200     cmp_rvecs(stdout, "boxv", DIM, st1->boxv, st2->boxv, FALSE, ftol, abstol);
 201     if (st1->flags & enumValueToBitMask(StateEntry::SVirPrev))
 202     {
 203         fprintf(stdout, "comparing shake vir_prev\n");
 204         cmp_rvecs(stdout, "svir_prev", DIM, st1->svir_prev, st2->svir_prev, FALSE, ftol, abstol);
 205     }
 206     if (st1->flags & enumValueToBitMask(StateEntry::FVirPrev))
 207     {
 208         fprintf(stdout, "comparing force vir_prev\n");
 209         cmp_rvecs(stdout, "fvir_prev", DIM, st1->fvir_prev, st2->fvir_prev, FALSE, ftol, abstol);
 210     }
 211     if (st1->flags & enumValueToBitMask(StateEntry::PressurePrevious))
 212     {
 213         fprintf(stdout, "comparing prev_pres\n");
 214         cmp_rvecs(stdout, "pres_prev", DIM, st1->pres_prev, st2->pres_prev, FALSE, ftol, abstol);
 215     }
 216     cmp_int(stdout, "ngtc", -1, st1->ngtc, st2->ngtc);
 217     cmp_int(stdout, "nhchainlength", -1, st1->nhchainlength, st2->nhchainlength);
 218     if (st1->ngtc == st2->ngtc && st1->nhchainlength == st2->nhchainlength)
 219     {
 220         for (i = 0; i < st1->ngtc; i++)
 221         {
 222             nc = i * st1->nhchainlength;
 223             for (j = 0; j < nc; j++)
 224             {
 225                 cmp_real(stdout, "nosehoover_xi", i, st1->nosehoover_xi[nc + j], st2->nosehoover_xi[nc + j], ftol, abstol);
 226             }
 227         }
 228     }
 229     cmp_int(stdout, "nnhpres", -1, st1->nnhpres, st2->nnhpres);
 230     if (st1->nnhpres == st2->nnhpres && st1->nhchainlength == st2->nhchainlength)
 231     {
 232         for (i = 0; i < st1->nnhpres; i++)
 233         {
 234             nc = i * st1->nhchainlength;
 235             for (j = 0; j < nc; j++)
 236             {
 237                 cmp_real(stdout, "nosehoover_xi", i, st1->nhpres_xi[nc + j], st2->nhpres_xi[nc + j], ftol, abstol);
 238             }
 239         }
 240     }
 241
 242     cmp_int(stdout, "natoms", -1, st1->natoms, st2->natoms);
 243     if (st1->natoms == st2->natoms)
 244     {
 245         if ((st1->flags & enumValueToBitMask(StateEntry::X))
 246             && (st2->flags & enumValueToBitMask(StateEntry::X)))
 247         {
 248             fprintf(stdout, "comparing x\n");
 249             cmp_rvecs(stdout, "x", st1->natoms, st1->x.rvec_array(), st2->x.rvec_array(), bRMSD, ftol, abstol);
 250         }
 251         if ((st1->flags & enumValueToBitMask(StateEntry::V))
 252             && (st2->flags & enumValueToBitMask(StateEntry::V)))
 253         {
 254             fprintf(stdout, "comparing v\n");
 255             cmp_rvecs(stdout, "v", st1->natoms, st1->v.rvec_array(), st2->v.rvec_array(), bRMSD, ftol, abstol);
 256         }
 257     }
 258 }
 259
 260 rvec* makeRvecArray(gmx::ArrayRef<const gmx::RVec> v, gmx::index n)
 261 {
 262     GMX_ASSERT(v.ssize() >= n, "We can't copy more elements than the vector size");
 263
 264     rvec* dest;
 265
 266     snew(dest, n);
 267
 268     const rvec* vPtr = as_rvec_array(v.data());
 269     for (gmx::index i = 0; i < n; i++)
 270     {
 271         copy_rvec(vPtr[i], dest[i]);
 272     }
 273
 274     return dest;
 275 }
 276
 277 t_state::t_state() :
 278     natoms(0),
 279     ngtc(0),
 280     nnhpres(0),
 281     nhchainlength(0),
 282     flags(0),
 283     fep_state(0),
 284     lambda(),
 285
 286     baros_integral(0),
 287     veta(0),
 288     vol0(0),
 289
 290     ekinstate(),
 291     hist(),
 292     dfhist(nullptr),
 293     awhHistory(nullptr),
 294     ddp_count(0),
 295     ddp_count_cg_gl(0)
 296
 297 {
 298     // It would be nicer to initialize these with {} or {{0}} in the
 299     // above initialization list, but uncrustify doesn't understand
 300     // that.
 301     // TODO Fix this if we switch to clang-format some time.
 302     lambda = { { 0 } };
 303     clear_mat(box);
 304     clear_mat(box_rel);
 305     clear_mat(boxv);
 306     clear_mat(pres_prev);
 307     clear_mat(svir_prev);
 308     clear_mat(fvir_prev);
 309 }
 310
 311 void set_box_rel(const t_inputrec* ir, t_state* state)
 312 {
 313     /* Make sure the box obeys the restrictions before we fix the ratios */
 314     correct_box(nullptr, 0, state->box);
 315
 316     clear_mat(state->box_rel);
 317
 318     if (inputrecPreserveShape(ir))
 319     {
 320         const int ndim = ir->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
 321         do_box_rel(ndim, ir->deform, state->box_rel, state->box, true);
 322     }
 323 }
 324
 325 void preserve_box_shape(const t_inputrec* ir, matrix box_rel, matrix box)
 326 {
 327     if (inputrecPreserveShape(ir))
 328     {
 329         const int ndim = ir->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
 330         do_box_rel(ndim, ir->deform, box_rel, box, false);
 331     }
 332 }
 333
 334 void printLambdaStateToLog(FILE* fplog, gmx::ArrayRef<const real> lambda, const bool isInitialOutput)
 335 {
 336     if (fplog != nullptr)
 337     {
 338         fprintf(fplog, "%s vector of lambda components:[ ", isInitialOutput ? "Initial" : "Current");
 339         for (const auto& l : lambda)
 340         {
 341             fprintf(fplog, "%10.4f ", l);
 342         }
 343         fprintf(fplog, "]\n%s", isInitialOutput ? "" : "\n");
 344     }
 345 }
 346
 347 void initialize_lambdas(FILE*                            fplog,
 348                         const FreeEnergyPerturbationType freeEnergyPerturbationType,
 349                         const bool                       haveSimulatedTempering,
 350                         const t_lambda&                  fep,
 351                         gmx::ArrayRef<const real>        simulatedTemperingTemps,
 352                         gmx::ArrayRef<real>              ref_t,
 353                         bool                             isMaster,
 354                         int*                             fep_state,
 355                         gmx::ArrayRef<real>              lambda)
 356 {
 357     /* TODO: Clean up initialization of fep_state and lambda in
 358        t_state.  This function works, but could probably use a logic
 359        rewrite to keep all the different types of efep straight. */
 360
 361     if ((freeEnergyPerturbationType == FreeEnergyPerturbationType::No) && (!haveSimulatedTempering))
 362     {
 363         return;
 364     }
 365
 366     if (isMaster)
 367     {
 368         *fep_state = fep.init_fep_state; /* this might overwrite the checkpoint
 369                                              if checkpoint is set -- a kludge is in for now
 370                                              to prevent this.*/
 371     }
 372
 373     for (int i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
 374     {
 375         double thisLambda;
 376         /* overwrite lambda state with init_lambda for now for backwards compatibility */
 377         if (fep.init_lambda >= 0) /* if it's -1, it was never initialized */
 378         {
 379             thisLambda = fep.init_lambda;
 380         }
 381         else
 382         {
 383             thisLambda = fep.all_lambda[i][fep.init_fep_state];
 384         }
 385         if (isMaster)
 386         {
 387             lambda[i] = thisLambda;
 388         }
 389     }
 390     if (haveSimulatedTempering)
 391     {
 392         /* need to rescale control temperatures to match current state */
 393         for (int i = 0; i < ref_t.ssize(); i++)
 394         {
 395             if (ref_t[i] > 0)
 396             {
 397                 ref_t[i] = simulatedTemperingTemps[fep.init_fep_state];
 398             }
 399         }
 400     }
 401
 402     /* Send to the log the information on the current lambdas */
 403     const bool isInitialOutput = true;
 404     printLambdaStateToLog(fplog, lambda, isInitialOutput);
 405 }