src/gromacs/modularsimulator/statepropagatordata.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief Defines the state for the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "statepropagatordata.h"
  45
  46 #include "gromacs/commandline/filenm.h"
  47 #include "gromacs/domdec/collect.h"
  48 #include "gromacs/domdec/domdec.h"
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  52 #include "gromacs/mdlib/mdatoms.h"
  53 #include "gromacs/mdlib/mdoutf.h"
  54 #include "gromacs/mdlib/stat.h"
  55 #include "gromacs/mdlib/update.h"
  56 #include "gromacs/mdtypes/checkpointdata.h"
  57 #include "gromacs/mdtypes/commrec.h"
  58 #include "gromacs/mdtypes/forcebuffers.h"
  59 #include "gromacs/mdtypes/forcerec.h"
  60 #include "gromacs/mdtypes/inputrec.h"
  61 #include "gromacs/mdtypes/mdatom.h"
  62 #include "gromacs/mdtypes/mdrunoptions.h"
  63 #include "gromacs/mdtypes/state.h"
  64 #include "gromacs/pbcutil/pbc.h"
  65 #include "gromacs/topology/atoms.h"
  66 #include "gromacs/topology/topology.h"
  67
  68 #include "freeenergyperturbationdata.h"
  69 #include "modularsimulator.h"
  70 #include "simulatoralgorithm.h"
  71
  72 namespace gmx
  73 {
  74 StatePropagatorData::StatePropagatorData(int                numAtoms,
  75                                          FILE*              fplog,
  76                                          const t_commrec*   cr,
  77                                          t_state*           globalState,
  78                                          bool               useGPU,
  79                                          bool               canMoleculesBeDistributedOverPBC,
  80                                          bool               writeFinalConfiguration,
  81                                          const std::string& finalConfigurationFilename,
  82                                          const t_inputrec*  inputrec,
  83                                          const t_mdatoms*   mdatoms,
  84                                          const gmx_mtop_t*  globalTop) :
  85     totalNumAtoms_(numAtoms),
  86     localNAtoms_(0),
  87     box_{ { 0 } },
  88     previousBox_{ { 0 } },
  89     ddpCount_(0),
  90     element_(std::make_unique<Element>(this,
  91                                        fplog,
  92                                        cr,
  93                                        inputrec->nstxout,
  94                                        inputrec->nstvout,
  95                                        inputrec->nstfout,
  96                                        inputrec->nstxout_compressed,
  97                                        canMoleculesBeDistributedOverPBC,
  98                                        writeFinalConfiguration,
  99                                        finalConfigurationFilename,
 100                                        inputrec,
 101                                        globalTop)),
 102     vvResetVelocities_(false),
 103     isRegularSimulationEnd_(false),
 104     lastStep_(-1),
 105     globalState_(globalState)
 106 {
 107     bool stateHasVelocities;
 108     // Local state only becomes valid now.
 109     if (DOMAINDECOMP(cr))
 110     {
 111         auto localState = std::make_unique<t_state>();
 112         dd_init_local_state(cr->dd, globalState, localState.get());
 113         stateHasVelocities = ((static_cast<unsigned int>(localState->flags) & (1U << estV)) != 0U);
 114         setLocalState(std::move(localState));
 115     }
 116     else
 117     {
 118         state_change_natoms(globalState, globalState->natoms);
 119         f_.resize(globalState->natoms);
 120         localNAtoms_ = globalState->natoms;
 121         x_           = globalState->x;
 122         v_           = globalState->v;
 123         copy_mat(globalState->box, box_);
 124         stateHasVelocities = ((static_cast<unsigned int>(globalState->flags) & (1U << estV)) != 0U);
 125         previousX_.resizeWithPadding(localNAtoms_);
 126         ddpCount_ = globalState->ddp_count;
 127         copyPosition();
 128     }
 129     if (useGPU)
 130     {
 131         changePinningPolicy(&x_, gmx::PinningPolicy::PinnedIfSupported);
 132     }
 133
 134     if (DOMAINDECOMP(cr) && MASTER(cr))
 135     {
 136         xGlobal_.resizeWithPadding(totalNumAtoms_);
 137         previousXGlobal_.resizeWithPadding(totalNumAtoms_);
 138         vGlobal_.resizeWithPadding(totalNumAtoms_);
 139         fGlobal_.resizeWithPadding(totalNumAtoms_);
 140     }
 141
 142     if (!inputrec->bContinuation)
 143     {
 144         if (stateHasVelocities)
 145         {
 146             auto v = velocitiesView().paddedArrayRef();
 147             // Set the velocities of vsites, shells and frozen atoms to zero
 148             for (int i = 0; i < mdatoms->homenr; i++)
 149             {
 150                 if (mdatoms->ptype[i] == eptVSite || mdatoms->ptype[i] == eptShell)
 151                 {
 152                     clear_rvec(v[i]);
 153                 }
 154                 else if (mdatoms->cFREEZE)
 155                 {
 156                     for (int m = 0; m < DIM; m++)
 157                     {
 158                         if (inputrec->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 159                         {
 160                             v[i][m] = 0;
 161                         }
 162                     }
 163                 }
 164             }
 165         }
 166         if (inputrec->eI == eiVV)
 167         {
 168             vvResetVelocities_ = true;
 169         }
 170     }
 171 }
 172
 173 StatePropagatorData::Element* StatePropagatorData::element()
 174 {
 175     return element_.get();
 176 }
 177
 178 void StatePropagatorData::setup()
 179 {
 180     if (element_)
 181     {
 182         element_->elementSetup();
 183     }
 184 }
 185
 186 ArrayRefWithPadding<RVec> StatePropagatorData::positionsView()
 187 {
 188     return x_.arrayRefWithPadding();
 189 }
 190
 191 ArrayRefWithPadding<const RVec> StatePropagatorData::constPositionsView() const
 192 {
 193     return x_.constArrayRefWithPadding();
 194 }
 195
 196 ArrayRefWithPadding<RVec> StatePropagatorData::previousPositionsView()
 197 {
 198     return previousX_.arrayRefWithPadding();
 199 }
 200
 201 ArrayRefWithPadding<const RVec> StatePropagatorData::constPreviousPositionsView() const
 202 {
 203     return previousX_.constArrayRefWithPadding();
 204 }
 205
 206 ArrayRefWithPadding<RVec> StatePropagatorData::velocitiesView()
 207 {
 208     return v_.arrayRefWithPadding();
 209 }
 210
 211 ArrayRefWithPadding<const RVec> StatePropagatorData::constVelocitiesView() const
 212 {
 213     return v_.constArrayRefWithPadding();
 214 }
 215
 216 ForceBuffersView& StatePropagatorData::forcesView()
 217 {
 218     return f_.view();
 219 }
 220
 221 const ForceBuffersView& StatePropagatorData::constForcesView() const
 222 {
 223     return f_.view();
 224 }
 225
 226 rvec* StatePropagatorData::box()
 227 {
 228     return box_;
 229 }
 230
 231 const rvec* StatePropagatorData::constBox() const
 232 {
 233     return box_;
 234 }
 235
 236 rvec* StatePropagatorData::previousBox()
 237 {
 238     return previousBox_;
 239 }
 240
 241 const rvec* StatePropagatorData::constPreviousBox() const
 242 {
 243     return previousBox_;
 244 }
 245
 246 int StatePropagatorData::localNumAtoms() const
 247 {
 248     return localNAtoms_;
 249 }
 250
 251 int StatePropagatorData::totalNumAtoms() const
 252 {
 253     return totalNumAtoms_;
 254 }
 255
 256 std::unique_ptr<t_state> StatePropagatorData::localState()
 257 {
 258     auto state   = std::make_unique<t_state>();
 259     state->flags = (1U << estX) | (1U << estV) | (1U << estBOX);
 260     state_change_natoms(state.get(), localNAtoms_);
 261     state->x = x_;
 262     state->v = v_;
 263     copy_mat(box_, state->box);
 264     state->ddp_count       = ddpCount_;
 265     state->ddp_count_cg_gl = ddpCountCgGl_;
 266     state->cg_gl           = cgGl_;
 267     return state;
 268 }
 269
 270 void StatePropagatorData::setLocalState(std::unique_ptr<t_state> state)
 271 {
 272     localNAtoms_ = state->natoms;
 273     x_.resizeWithPadding(localNAtoms_);
 274     previousX_.resizeWithPadding(localNAtoms_);
 275     v_.resizeWithPadding(localNAtoms_);
 276     x_ = state->x;
 277     v_ = state->v;
 278     copy_mat(state->box, box_);
 279     copyPosition();
 280     ddpCount_     = state->ddp_count;
 281     ddpCountCgGl_ = state->ddp_count_cg_gl;
 282     cgGl_         = state->cg_gl;
 283
 284     if (vvResetVelocities_)
 285     {
 286         /* DomDec runs twice early in the simulation, once at setup time, and once before the first
 287          * step. Every time DD runs, it sets a new local state here. We are saving a backup during
 288          * setup time (ok for non-DD cases), so we need to update our backup to the DD state before
 289          * the first step here to avoid resetting to an earlier DD state. This is done before any
 290          * propagation that needs to be reset, so it's not very safe but correct for now.
 291          * TODO: Get rid of this once input is assumed to be at half steps
 292          */
 293         velocityBackup_ = v_;
 294     }
 295 }
 296
 297 t_state* StatePropagatorData::globalState()
 298 {
 299     return globalState_;
 300 }
 301
 302 ForceBuffers* StatePropagatorData::forcePointer()
 303 {
 304     return &f_;
 305 }
 306
 307 void StatePropagatorData::copyPosition()
 308 {
 309     int nth = gmx_omp_nthreads_get(emntUpdate);
 310
 311 #pragma omp parallel for num_threads(nth) schedule(static) default(none) shared(nth)
 312     for (int th = 0; th < nth; th++)
 313     {
 314         int start_th, end_th;
 315         getThreadAtomRange(nth, th, localNAtoms_, &start_th, &end_th);
 316         copyPosition(start_th, end_th);
 317     }
 318
 319     /* Box is changed in update() when we do pressure coupling,
 320      * but we should still use the old box for energy corrections and when
 321      * writing it to the energy file, so it matches the trajectory files for
 322      * the same timestep above. Make a copy in a separate array.
 323      */
 324     copy_mat(box_, previousBox_);
 325 }
 326
 327 void StatePropagatorData::copyPosition(int start, int end)
 328 {
 329     for (int i = start; i < end; ++i)
 330     {
 331         previousX_[i] = x_[i];
 332     }
 333 }
 334
 335 void StatePropagatorData::Element::scheduleTask(Step step,
 336                                                 Time gmx_unused            time,
 337                                                 const RegisterRunFunction& registerRunFunction)
 338 {
 339     if (statePropagatorData_->vvResetVelocities_)
 340     {
 341         statePropagatorData_->vvResetVelocities_ = false;
 342         registerRunFunction([this]() { statePropagatorData_->resetVelocities(); });
 343     }
 344     // copy x -> previousX
 345     registerRunFunction([this]() { statePropagatorData_->copyPosition(); });
 346     // if it's a write out step, keep a copy for writeout
 347     if (step == writeOutStep_ || (step == lastStep_ && writeFinalConfiguration_))
 348     {
 349         registerRunFunction([this]() { saveState(); });
 350     }
 351 }
 352
 353 void StatePropagatorData::Element::saveState()
 354 {
 355     GMX_ASSERT(!localStateBackup_, "Save state called again before previous state was written.");
 356     localStateBackup_ = statePropagatorData_->localState();
 357     if (freeEnergyPerturbationData_)
 358     {
 359         localStateBackup_->fep_state = freeEnergyPerturbationData_->currentFEPState();
 360         for (unsigned long i = 0; i < localStateBackup_->lambda.size(); ++i)
 361         {
 362             localStateBackup_->lambda[i] = freeEnergyPerturbationData_->constLambdaView()[i];
 363         }
 364         localStateBackup_->flags |= (1U << estLAMBDA) | (1U << estFEPSTATE);
 365     }
 366 }
 367
 368 std::optional<SignallerCallback> StatePropagatorData::Element::registerTrajectorySignallerCallback(TrajectoryEvent event)
 369 {
 370     if (event == TrajectoryEvent::StateWritingStep)
 371     {
 372         return [this](Step step, Time /*unused*/) { this->writeOutStep_ = step; };
 373     }
 374     return std::nullopt;
 375 }
 376
 377 std::optional<ITrajectoryWriterCallback>
 378 StatePropagatorData::Element::registerTrajectoryWriterCallback(TrajectoryEvent event)
 379 {
 380     if (event == TrajectoryEvent::StateWritingStep)
 381     {
 382         return [this](gmx_mdoutf* outf, Step step, Time time, bool writeTrajectory, bool gmx_unused writeLog) {
 383             if (writeTrajectory)
 384             {
 385                 write(outf, step, time);
 386             }
 387         };
 388     }
 389     return std::nullopt;
 390 }
 391
 392 void StatePropagatorData::Element::write(gmx_mdoutf_t outf, Step currentStep, Time currentTime)
 393 {
 394     wallcycle_start(mdoutf_get_wcycle(outf), ewcTRAJ);
 395     unsigned int mdof_flags = 0;
 396     if (do_per_step(currentStep, nstxout_))
 397     {
 398         mdof_flags |= MDOF_X;
 399     }
 400     if (do_per_step(currentStep, nstvout_))
 401     {
 402         mdof_flags |= MDOF_V;
 403     }
 404     if (do_per_step(currentStep, nstfout_))
 405     {
 406         mdof_flags |= MDOF_F;
 407     }
 408     if (do_per_step(currentStep, nstxout_compressed_))
 409     {
 410         mdof_flags |= MDOF_X_COMPRESSED;
 411     }
 412     if (do_per_step(currentStep, mdoutf_get_tng_box_output_interval(outf)))
 413     {
 414         mdof_flags |= MDOF_BOX;
 415     }
 416     if (do_per_step(currentStep, mdoutf_get_tng_lambda_output_interval(outf)))
 417     {
 418         mdof_flags |= MDOF_LAMBDA;
 419     }
 420     if (do_per_step(currentStep, mdoutf_get_tng_compressed_box_output_interval(outf)))
 421     {
 422         mdof_flags |= MDOF_BOX_COMPRESSED;
 423     }
 424     if (do_per_step(currentStep, mdoutf_get_tng_compressed_lambda_output_interval(outf)))
 425     {
 426         mdof_flags |= MDOF_LAMBDA_COMPRESSED;
 427     }
 428
 429     if (mdof_flags == 0)
 430     {
 431         wallcycle_stop(mdoutf_get_wcycle(outf), ewcTRAJ);
 432         return;
 433     }
 434     GMX_ASSERT(localStateBackup_, "Trajectory writing called, but no state saved.");
 435
 436     // TODO: This is only used for CPT - needs to be filled when we turn CPT back on
 437     ObservablesHistory* observablesHistory = nullptr;
 438
 439     mdoutf_write_to_trajectory_files(
 440             fplog_, cr_, outf, static_cast<int>(mdof_flags), statePropagatorData_->totalNumAtoms_,
 441             currentStep, currentTime, localStateBackup_.get(), statePropagatorData_->globalState_,
 442             observablesHistory, statePropagatorData_->f_.view().force(), &dummyCheckpointDataHolder_);
 443
 444     if (currentStep != lastStep_ || !isRegularSimulationEnd_)
 445     {
 446         localStateBackup_.reset();
 447     }
 448     wallcycle_stop(mdoutf_get_wcycle(outf), ewcTRAJ);
 449 }
 450
 451 void StatePropagatorData::Element::elementSetup()
 452 {
 453     if (statePropagatorData_->vvResetVelocities_)
 454     {
 455         // MD-VV does the first velocity half-step only to calculate the constraint virial,
 456         // then resets the velocities since the input is assumed to be positions and velocities
 457         // at full time step. TODO: Change this to have input at half time steps.
 458         statePropagatorData_->velocityBackup_ = statePropagatorData_->v_;
 459     }
 460 }
 461
 462 void StatePropagatorData::resetVelocities()
 463 {
 464     v_ = velocityBackup_;
 465 }
 466
 467 namespace
 468 {
 469 /*!
 470  * \brief Enum describing the contents StatePropagatorData::Element writes to modular checkpoint
 471  *
 472  * When changing the checkpoint content, add a new element just above Count, and adjust the
 473  * checkpoint functionality.
 474  */
 475 enum class CheckpointVersion
 476 {
 477     Base, //!< First version of modular checkpointing
 478     Count //!< Number of entries. Add new versions right above this!
 479 };
 480 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
 481 } // namespace
 482
 483 template<CheckpointDataOperation operation>
 484 void StatePropagatorData::Element::doCheckpointData(CheckpointData<operation>* checkpointData,
 485                                                     const t_commrec*           cr)
 486 {
 487     ArrayRef<RVec> xGlobalRef;
 488     ArrayRef<RVec> vGlobalRef;
 489     if (DOMAINDECOMP(cr))
 490     {
 491         if (MASTER(cr))
 492         {
 493             xGlobalRef = statePropagatorData_->xGlobal_;
 494             vGlobalRef = statePropagatorData_->vGlobal_;
 495         }
 496         if (operation == CheckpointDataOperation::Write)
 497         {
 498             dd_collect_vec(cr->dd, statePropagatorData_->ddpCount_, statePropagatorData_->ddpCountCgGl_,
 499                            statePropagatorData_->cgGl_, statePropagatorData_->x_, xGlobalRef);
 500             dd_collect_vec(cr->dd, statePropagatorData_->ddpCount_, statePropagatorData_->ddpCountCgGl_,
 501                            statePropagatorData_->cgGl_, statePropagatorData_->v_, vGlobalRef);
 502         }
 503     }
 504     else
 505     {
 506         xGlobalRef = statePropagatorData_->x_;
 507         vGlobalRef = statePropagatorData_->v_;
 508     }
 509     if (MASTER(cr))
 510     {
 511         GMX_ASSERT(checkpointData, "Master needs a valid pointer to a CheckpointData object");
 512         checkpointVersion(checkpointData, "StatePropagatorData version", c_currentVersion);
 513
 514         checkpointData->arrayRef("positions", makeCheckpointArrayRef<operation>(xGlobalRef));
 515         checkpointData->arrayRef("velocities", makeCheckpointArrayRef<operation>(vGlobalRef));
 516         checkpointData->tensor("box", statePropagatorData_->box_);
 517         checkpointData->scalar("ddpCount", &statePropagatorData_->ddpCount_);
 518         checkpointData->scalar("ddpCountCgGl", &statePropagatorData_->ddpCountCgGl_);
 519         checkpointData->arrayRef("cgGl", makeCheckpointArrayRef<operation>(statePropagatorData_->cgGl_));
 520     }
 521 }
 522
 523 void StatePropagatorData::Element::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
 524                                                        const t_commrec*                   cr)
 525 {
 526     doCheckpointData<CheckpointDataOperation::Write>(
 527             checkpointData ? &checkpointData.value() : nullptr, cr);
 528 }
 529
 530 /*!
 531  * \brief Update the legacy global state
 532  *
 533  * When restoring from checkpoint, data will be distributed during domain decomposition at setup stage.
 534  * Domain decomposition still uses the legacy global t_state object so make sure it's up-to-date.
 535  */
 536 static void updateGlobalState(t_state*                      globalState,
 537                               const PaddedHostVector<RVec>& x,
 538                               const PaddedHostVector<RVec>& v,
 539                               const tensor                  box,
 540                               int                           ddpCount,
 541                               int                           ddpCountCgGl,
 542                               const std::vector<int>&       cgGl)
 543 {
 544     globalState->x = x;
 545     globalState->v = v;
 546     copy_mat(box, globalState->box);
 547     globalState->ddp_count       = ddpCount;
 548     globalState->ddp_count_cg_gl = ddpCountCgGl;
 549     globalState->cg_gl           = cgGl;
 550 }
 551
 552 void StatePropagatorData::Element::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
 553                                                           const t_commrec*                  cr)
 554 {
 555     doCheckpointData<CheckpointDataOperation::Read>(checkpointData ? &checkpointData.value() : nullptr, cr);
 556
 557     // Copy data to global state to be distributed by DD at setup stage
 558     if (DOMAINDECOMP(cr) && MASTER(cr))
 559     {
 560         updateGlobalState(statePropagatorData_->globalState_, statePropagatorData_->xGlobal_,
 561                           statePropagatorData_->vGlobal_, statePropagatorData_->box_,
 562                           statePropagatorData_->ddpCount_, statePropagatorData_->ddpCountCgGl_,
 563                           statePropagatorData_->cgGl_);
 564     }
 565 }
 566
 567 const std::string& StatePropagatorData::Element::clientID()
 568 {
 569     return identifier_;
 570 }
 571
 572 void StatePropagatorData::Element::trajectoryWriterTeardown(gmx_mdoutf* gmx_unused outf)
 573 {
 574     // Note that part of this code is duplicated in do_md_trajectory_writing.
 575     // This duplication is needed while both legacy and modular code paths are in use.
 576     // TODO: Remove duplication asap, make sure to keep in sync in the meantime.
 577     if (!writeFinalConfiguration_ || !isRegularSimulationEnd_)
 578     {
 579         return;
 580     }
 581
 582     GMX_ASSERT(localStateBackup_, "Final trajectory writing called, but no state saved.");
 583
 584     wallcycle_start(mdoutf_get_wcycle(outf), ewcTRAJ);
 585     if (DOMAINDECOMP(cr_))
 586     {
 587         auto globalXRef =
 588                 MASTER(cr_) ? statePropagatorData_->globalState_->x : gmx::ArrayRef<gmx::RVec>();
 589         dd_collect_vec(cr_->dd, localStateBackup_->ddp_count, localStateBackup_->ddp_count_cg_gl,
 590                        localStateBackup_->cg_gl, localStateBackup_->x, globalXRef);
 591         auto globalVRef =
 592                 MASTER(cr_) ? statePropagatorData_->globalState_->v : gmx::ArrayRef<gmx::RVec>();
 593         dd_collect_vec(cr_->dd, localStateBackup_->ddp_count, localStateBackup_->ddp_count_cg_gl,
 594                        localStateBackup_->cg_gl, localStateBackup_->v, globalVRef);
 595     }
 596     else
 597     {
 598         // We have the whole state locally: copy the local state pointer
 599         statePropagatorData_->globalState_ = localStateBackup_.get();
 600     }
 601
 602     if (MASTER(cr_))
 603     {
 604         fprintf(stderr, "\nWriting final coordinates.\n");
 605         if (canMoleculesBeDistributedOverPBC_ && !systemHasPeriodicMolecules_)
 606         {
 607             // Make molecules whole only for confout writing
 608             do_pbc_mtop(pbcType_, localStateBackup_->box, top_global_,
 609                         statePropagatorData_->globalState_->x.rvec_array());
 610         }
 611         write_sto_conf_mtop(finalConfigurationFilename_.c_str(), *top_global_->name, top_global_,
 612                             statePropagatorData_->globalState_->x.rvec_array(),
 613                             statePropagatorData_->globalState_->v.rvec_array(), pbcType_,
 614                             localStateBackup_->box);
 615     }
 616     wallcycle_stop(mdoutf_get_wcycle(outf), ewcTRAJ);
 617 }
 618
 619 std::optional<SignallerCallback> StatePropagatorData::Element::registerLastStepCallback()
 620 {
 621     return [this](Step step, Time /*time*/) {
 622         lastStep_               = step;
 623         isRegularSimulationEnd_ = (step == lastPlannedStep_);
 624     };
 625 }
 626
 627 StatePropagatorData::Element::Element(StatePropagatorData* statePropagatorData,
 628                                       FILE*                fplog,
 629                                       const t_commrec*     cr,
 630                                       int                  nstxout,
 631                                       int                  nstvout,
 632                                       int                  nstfout,
 633                                       int                  nstxout_compressed,
 634                                       bool                 canMoleculesBeDistributedOverPBC,
 635                                       bool                 writeFinalConfiguration,
 636                                       std::string          finalConfigurationFilename,
 637                                       const t_inputrec*    inputrec,
 638                                       const gmx_mtop_t*    globalTop) :
 639     statePropagatorData_(statePropagatorData),
 640     nstxout_(nstxout),
 641     nstvout_(nstvout),
 642     nstfout_(nstfout),
 643     nstxout_compressed_(nstxout_compressed),
 644     writeOutStep_(-1),
 645     freeEnergyPerturbationData_(nullptr),
 646     isRegularSimulationEnd_(false),
 647     lastStep_(-1),
 648     canMoleculesBeDistributedOverPBC_(canMoleculesBeDistributedOverPBC),
 649     systemHasPeriodicMolecules_(inputrec->bPeriodicMols),
 650     pbcType_(inputrec->pbcType),
 651     lastPlannedStep_(inputrec->nsteps + inputrec->init_step),
 652     writeFinalConfiguration_(writeFinalConfiguration),
 653     finalConfigurationFilename_(std::move(finalConfigurationFilename)),
 654     fplog_(fplog),
 655     cr_(cr),
 656     top_global_(globalTop)
 657 {
 658 }
 659 void StatePropagatorData::Element::setFreeEnergyPerturbationData(FreeEnergyPerturbationData* freeEnergyPerturbationData)
 660 {
 661     freeEnergyPerturbationData_ = freeEnergyPerturbationData;
 662 }
 663
 664 ISimulatorElement* StatePropagatorData::Element::getElementPointerImpl(
 665         LegacySimulatorData gmx_unused*        legacySimulatorData,
 666         ModularSimulatorAlgorithmBuilderHelper gmx_unused* builderHelper,
 667         StatePropagatorData*                               statePropagatorData,
 668         EnergyData gmx_unused*      energyData,
 669         FreeEnergyPerturbationData* freeEnergyPerturbationData,
 670         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper)
 671 {
 672     statePropagatorData->element()->setFreeEnergyPerturbationData(freeEnergyPerturbationData);
 673     return statePropagatorData->element();
 674 }
 675
 676 } // namespace gmx