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 namespace
 177 {
 178 /*!
 179  * \brief Enum describing the contents df_history_t writes to modular checkpoint
 180  *
 181  * When changing the checkpoint content, add a new element just above Count, and adjust the
 182  * checkpoint functionality.
 183  */
 184 enum class DFHistoryCheckpointVersion
 185 {
 186     Base, //!< First version of modular checkpointing
 187     Count //!< Number of entries. Add new versions right above this!
 188 };
 189 constexpr auto c_dfHistoryCurrentVersion =
 190         DFHistoryCheckpointVersion(int(DFHistoryCheckpointVersion::Count) - 1);
 191 } // namespace
 192
 193 template<gmx::CheckpointDataOperation operation>
 194 void df_history_t::doCheckpoint(gmx::CheckpointData<operation> checkpointData, LambdaWeightCalculation elamstats)
 195 {
 196     gmx::checkpointVersion(&checkpointData, "df_history_t version", c_dfHistoryCurrentVersion);
 197
 198     auto numLambdas = nlambda;
 199     checkpointData.scalar("nlambda", &numLambdas);
 200     if (operation == gmx::CheckpointDataOperation::Read)
 201     {
 202         // If this isn't matching, we haven't allocated the right amount of data
 203         GMX_RELEASE_ASSERT(numLambdas == nlambda,
 204                            "df_history_t checkpoint reading: Lambda vectors size mismatch.");
 205     }
 206
 207     checkpointData.scalar("bEquil", &bEquil);
 208     checkpointData.arrayRef("n_at_lam", gmx::makeCheckpointArrayRefFromArray<operation>(n_at_lam, nlambda));
 209
 210     checkpointData.arrayRef("sum_weights",
 211                             gmx::makeCheckpointArrayRefFromArray<operation>(sum_weights, nlambda));
 212     checkpointData.arrayRef("sum_dg", gmx::makeCheckpointArrayRefFromArray<operation>(sum_dg, nlambda));
 213     for (int idx = 0; idx < nlambda; idx++)
 214     {
 215         checkpointData.arrayRef(gmx::formatString("Tij[%d]", idx),
 216                                 gmx::makeCheckpointArrayRefFromArray<operation>(Tij[idx], nlambda));
 217         checkpointData.arrayRef(
 218                 gmx::formatString("Tij_empirical[%d]", idx),
 219                 gmx::makeCheckpointArrayRefFromArray<operation>(Tij_empirical[idx], nlambda));
 220     }
 221
 222     if (EWL(elamstats))
 223     {
 224         checkpointData.arrayRef("wl_histo",
 225                                 gmx::makeCheckpointArrayRefFromArray<operation>(wl_histo, nlambda));
 226         checkpointData.scalar("wl_delta", &wl_delta);
 227     }
 228     if ((elamstats == LambdaWeightCalculation::Minvar) || (elamstats == LambdaWeightCalculation::Barker)
 229         || (elamstats == LambdaWeightCalculation::Metropolis))
 230     {
 231         checkpointData.arrayRef(
 232                 "sum_minvar", gmx::makeCheckpointArrayRefFromArray<operation>(sum_minvar, nlambda));
 233         checkpointData.arrayRef("sum_variance",
 234                                 gmx::makeCheckpointArrayRefFromArray<operation>(sum_variance, nlambda));
 235         for (int idx = 0; idx < nlambda; idx++)
 236         {
 237             checkpointData.arrayRef(gmx::formatString("accum_p[%d]", idx),
 238                                     gmx::makeCheckpointArrayRefFromArray<operation>(accum_p[idx], nlambda));
 239             checkpointData.arrayRef(gmx::formatString("accum_m[%d]", idx),
 240                                     gmx::makeCheckpointArrayRefFromArray<operation>(accum_m[idx], nlambda));
 241             checkpointData.arrayRef(
 242                     gmx::formatString("accum_p2[%d]", idx),
 243                     gmx::makeCheckpointArrayRefFromArray<operation>(accum_p2[idx], nlambda));
 244             checkpointData.arrayRef(
 245                     gmx::formatString("accum_m2[%d]", idx),
 246                     gmx::makeCheckpointArrayRefFromArray<operation>(accum_m2[idx], nlambda));
 247         }
 248     }
 249 }
 250
 251 // explicit template instantiation
 252 template void df_history_t::doCheckpoint(gmx::CheckpointData<gmx::CheckpointDataOperation::Read> checkpointData,
 253                                          LambdaWeightCalculation elamstats);
 254 template void df_history_t::doCheckpoint(gmx::CheckpointData<gmx::CheckpointDataOperation::Write> checkpointData,
 255                                          LambdaWeightCalculation elamstats);
 256
 257 void init_dfhist_state(t_state* state, int dfhistNumLambda)
 258 {
 259     if (dfhistNumLambda > 0)
 260     {
 261         snew(state->dfhist, 1);
 262         init_df_history(state->dfhist, dfhistNumLambda);
 263     }
 264     else
 265     {
 266         state->dfhist = nullptr;
 267     }
 268 }
 269
 270 void comp_state(const t_state* st1, const t_state* st2, gmx_bool bRMSD, real ftol, real abstol)
 271 {
 272     int i, j, nc;
 273
 274     fprintf(stdout, "comparing flags\n");
 275     cmp_int(stdout, "flags", -1, st1->flags, st2->flags);
 276     fprintf(stdout, "comparing box\n");
 277     cmp_rvecs(stdout, "box", DIM, st1->box, st2->box, FALSE, ftol, abstol);
 278     fprintf(stdout, "comparing box_rel\n");
 279     cmp_rvecs(stdout, "box_rel", DIM, st1->box_rel, st2->box_rel, FALSE, ftol, abstol);
 280     fprintf(stdout, "comparing boxv\n");
 281     cmp_rvecs(stdout, "boxv", DIM, st1->boxv, st2->boxv, FALSE, ftol, abstol);
 282     if (st1->flags & enumValueToBitMask(StateEntry::SVirPrev))
 283     {
 284         fprintf(stdout, "comparing shake vir_prev\n");
 285         cmp_rvecs(stdout, "svir_prev", DIM, st1->svir_prev, st2->svir_prev, FALSE, ftol, abstol);
 286     }
 287     if (st1->flags & enumValueToBitMask(StateEntry::FVirPrev))
 288     {
 289         fprintf(stdout, "comparing force vir_prev\n");
 290         cmp_rvecs(stdout, "fvir_prev", DIM, st1->fvir_prev, st2->fvir_prev, FALSE, ftol, abstol);
 291     }
 292     if (st1->flags & enumValueToBitMask(StateEntry::PressurePrevious))
 293     {
 294         fprintf(stdout, "comparing prev_pres\n");
 295         cmp_rvecs(stdout, "pres_prev", DIM, st1->pres_prev, st2->pres_prev, FALSE, ftol, abstol);
 296     }
 297     cmp_int(stdout, "ngtc", -1, st1->ngtc, st2->ngtc);
 298     cmp_int(stdout, "nhchainlength", -1, st1->nhchainlength, st2->nhchainlength);
 299     if (st1->ngtc == st2->ngtc && st1->nhchainlength == st2->nhchainlength)
 300     {
 301         for (i = 0; i < st1->ngtc; i++)
 302         {
 303             nc = i * st1->nhchainlength;
 304             for (j = 0; j < nc; j++)
 305             {
 306                 cmp_real(stdout, "nosehoover_xi", i, st1->nosehoover_xi[nc + j], st2->nosehoover_xi[nc + j], ftol, abstol);
 307             }
 308         }
 309     }
 310     cmp_int(stdout, "nnhpres", -1, st1->nnhpres, st2->nnhpres);
 311     if (st1->nnhpres == st2->nnhpres && st1->nhchainlength == st2->nhchainlength)
 312     {
 313         for (i = 0; i < st1->nnhpres; i++)
 314         {
 315             nc = i * st1->nhchainlength;
 316             for (j = 0; j < nc; j++)
 317             {
 318                 cmp_real(stdout, "nosehoover_xi", i, st1->nhpres_xi[nc + j], st2->nhpres_xi[nc + j], ftol, abstol);
 319             }
 320         }
 321     }
 322
 323     cmp_int(stdout, "natoms", -1, st1->natoms, st2->natoms);
 324     if (st1->natoms == st2->natoms)
 325     {
 326         if ((st1->flags & enumValueToBitMask(StateEntry::X))
 327             && (st2->flags & enumValueToBitMask(StateEntry::X)))
 328         {
 329             fprintf(stdout, "comparing x\n");
 330             cmp_rvecs(stdout, "x", st1->natoms, st1->x.rvec_array(), st2->x.rvec_array(), bRMSD, ftol, abstol);
 331         }
 332         if ((st1->flags & enumValueToBitMask(StateEntry::V))
 333             && (st2->flags & enumValueToBitMask(StateEntry::V)))
 334         {
 335             fprintf(stdout, "comparing v\n");
 336             cmp_rvecs(stdout, "v", st1->natoms, st1->v.rvec_array(), st2->v.rvec_array(), bRMSD, ftol, abstol);
 337         }
 338     }
 339 }
 340
 341 rvec* makeRvecArray(gmx::ArrayRef<const gmx::RVec> v, gmx::index n)
 342 {
 343     GMX_ASSERT(v.ssize() >= n, "We can't copy more elements than the vector size");
 344
 345     rvec* dest;
 346
 347     snew(dest, n);
 348
 349     const rvec* vPtr = as_rvec_array(v.data());
 350     for (gmx::index i = 0; i < n; i++)
 351     {
 352         copy_rvec(vPtr[i], dest[i]);
 353     }
 354
 355     return dest;
 356 }
 357
 358 t_state::t_state() :
 359     natoms(0),
 360     ngtc(0),
 361     nnhpres(0),
 362     nhchainlength(0),
 363     flags(0),
 364     fep_state(0),
 365     lambda(),
 366
 367     baros_integral(0),
 368     veta(0),
 369     vol0(0),
 370
 371     ekinstate(),
 372     hist(),
 373     dfhist(nullptr),
 374     awhHistory(nullptr),
 375     ddp_count(0),
 376     ddp_count_cg_gl(0)
 377
 378 {
 379     // It would be nicer to initialize these with {} or {{0}} in the
 380     // above initialization list, but uncrustify doesn't understand
 381     // that.
 382     // TODO Fix this if we switch to clang-format some time.
 383     lambda = { { 0 } };
 384     clear_mat(box);
 385     clear_mat(box_rel);
 386     clear_mat(boxv);
 387     clear_mat(pres_prev);
 388     clear_mat(svir_prev);
 389     clear_mat(fvir_prev);
 390 }
 391
 392 void set_box_rel(const t_inputrec* ir, t_state* state)
 393 {
 394     /* Make sure the box obeys the restrictions before we fix the ratios */
 395     correct_box(nullptr, 0, state->box);
 396
 397     clear_mat(state->box_rel);
 398
 399     if (inputrecPreserveShape(ir))
 400     {
 401         const int ndim = ir->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
 402         do_box_rel(ndim, ir->deform, state->box_rel, state->box, true);
 403     }
 404 }
 405
 406 void preserve_box_shape(const t_inputrec* ir, matrix box_rel, matrix box)
 407 {
 408     if (inputrecPreserveShape(ir))
 409     {
 410         const int ndim = ir->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
 411         do_box_rel(ndim, ir->deform, box_rel, box, false);
 412     }
 413 }
 414
 415 void printLambdaStateToLog(FILE* fplog, gmx::ArrayRef<const real> lambda, const bool isInitialOutput)
 416 {
 417     if (fplog != nullptr)
 418     {
 419         fprintf(fplog, "%s vector of lambda components:[ ", isInitialOutput ? "Initial" : "Current");
 420         for (const auto& l : lambda)
 421         {
 422             fprintf(fplog, "%10.4f ", l);
 423         }
 424         fprintf(fplog, "]\n%s", isInitialOutput ? "" : "\n");
 425     }
 426 }
 427
 428 void initialize_lambdas(FILE*                            fplog,
 429                         const FreeEnergyPerturbationType freeEnergyPerturbationType,
 430                         const bool                       haveSimulatedTempering,
 431                         const t_lambda&                  fep,
 432                         gmx::ArrayRef<const real>        simulatedTemperingTemps,
 433                         gmx::ArrayRef<real>              ref_t,
 434                         bool                             isMaster,
 435                         int*                             fep_state,
 436                         gmx::ArrayRef<real>              lambda)
 437 {
 438     /* TODO: Clean up initialization of fep_state and lambda in
 439        t_state.  This function works, but could probably use a logic
 440        rewrite to keep all the different types of efep straight. */
 441
 442     if ((freeEnergyPerturbationType == FreeEnergyPerturbationType::No) && (!haveSimulatedTempering))
 443     {
 444         return;
 445     }
 446
 447     if (isMaster)
 448     {
 449         *fep_state = fep.init_fep_state; /* this might overwrite the checkpoint
 450                                              if checkpoint is set -- a kludge is in for now
 451                                              to prevent this.*/
 452     }
 453
 454     for (int i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
 455     {
 456         double thisLambda;
 457         /* overwrite lambda state with init_lambda for now for backwards compatibility */
 458         if (fep.init_lambda >= 0) /* if it's -1, it was never initialized */
 459         {
 460             thisLambda = fep.init_lambda;
 461         }
 462         else
 463         {
 464             thisLambda = fep.all_lambda[i][fep.init_fep_state];
 465         }
 466         if (isMaster)
 467         {
 468             lambda[i] = thisLambda;
 469         }
 470     }
 471     if (haveSimulatedTempering)
 472     {
 473         /* need to rescale control temperatures to match current state */
 474         for (int i = 0; i < ref_t.ssize(); i++)
 475         {
 476             if (ref_t[i] > 0)
 477             {
 478                 ref_t[i] = simulatedTemperingTemps[fep.init_fep_state];
 479             }
 480         }
 481     }
 482
 483     /* Send to the log the information on the current lambdas */
 484     const bool isInitialOutput = true;
 485     printLambdaStateToLog(fplog, lambda, isInitialOutput);
 486 }