Adding energy output field for density fitting
authorChristian Blau <cblau@gwdg.de>
Mon, 2 Sep 2019 13:47:24 +0000 (15:47 +0200)
committerErik Lindahl <erik.lindahl@gmail.com>
Fri, 6 Sep 2019 18:56:36 +0000 (20:56 +0200)
Added a new energy output field for energies that stem from the density
fitting code.

refs #2282

Change-Id: I78c64213f4b958f00a938f0a67c586c2f4e9033e

16 files changed:
src/gromacs/applied_forces/densityfitting.cpp
src/gromacs/applied_forces/densityfittingforceprovider.cpp
src/gromacs/fileio/tpxio.cpp
src/gromacs/listed_forces/bonded.cpp
src/gromacs/mdlib/energyoutput.cpp
src/gromacs/mdlib/energyoutput.h
src/gromacs/mdlib/tests/energyoutput.cpp
src/gromacs/mdrun/md.cpp
src/gromacs/mdrun/mimic.cpp
src/gromacs/mdrun/minimize.cpp
src/gromacs/mdrun/rerun.cpp
src/gromacs/modularsimulator/energyelement.cpp
src/gromacs/modularsimulator/energyelement.h
src/gromacs/topology/ifunc.cpp
src/gromacs/topology/ifunc.h
src/gromacs/utility/mdmodulenotification.h

index 962b2e8aa69df9c43b819f39925c9ce62b5e8645..0d122c2e1f2c3b17a263e946d1ad745819fd96f5 100644 (file)
@@ -249,6 +249,13 @@ class DensityFitting final : public IMDModule
                     this->densityFittingSimulationParameters_.setPeriodicBoundaryConditionType(pbc);
                 };
             notifier->notifier_.subscribe(setPeriodicBoundaryContionsFunction);
+
+            // adding output to energy file
+            const auto requestEnergyOutput
+                = [this](MdModulesEnergyOutputToDensityFittingRequestChecker *energyOutputRequest) {
+                        this->setEnergyOutputRequest(energyOutputRequest);
+                    };
+            notifier->notifier_.subscribe(requestEnergyOutput);
         }
 
         //! From IMDModule
@@ -283,8 +290,19 @@ class DensityFitting final : public IMDModule
          */
         void constructLocalAtomSet(LocalAtomSetManager * localAtomSetManager)
         {
-            LocalAtomSet atomSet = localAtomSetManager->add(densityFittingOptions_.buildParameters().indices_);
-            densityFittingSimulationParameters_.setLocalAtomSet(atomSet);
+            if (densityFittingOptions_.active())
+            {
+                LocalAtomSet atomSet
+                    = localAtomSetManager->add(densityFittingOptions_.buildParameters().indices_);
+                densityFittingSimulationParameters_.setLocalAtomSet(atomSet);
+            }
+        }
+
+        /*! \brief Request energy output to energy file during simulation.
+         */
+        void setEnergyOutputRequest(MdModulesEnergyOutputToDensityFittingRequestChecker *energyOutputRequest)
+        {
+            energyOutputRequest->energyOutputToDensityFitting_ = densityFittingOptions_.active();
         }
 
     private:
index 43f556032edafcc63adce6a6cf1a155054138dfb..57ccd175cc86dbb0523b07346f7be15c248f05df 100644 (file)
@@ -224,7 +224,7 @@ void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput
     // calculate corresponding potential energy
     const float similarity  = measure_.similarity(gaussTransform_.constView());
     const real  energy      = -similarity * parameters_.forceConstant_;
-    forceProviderOutput->enerd_.term[F_COM_PULL] += energy;
+    forceProviderOutput->enerd_.term[F_DENSITYFITTING] += energy;
 }
 
 /********************************************************************
index 6a4ea30d3aee473443b1f7b1a921fc0e87799216..c58d3749023341dbfcec14615b7051ead3af2d49 100644 (file)
@@ -206,6 +206,7 @@ static const t_ftupd ftupd[] = {
     { 79, F_DVDL_BONDED,      },
     { 79, F_DVDL_RESTRAINT    },
     { 79, F_DVDL_TEMPERATURE  },
+    { tpxv_GenericInternalParameters, F_DENSITYFITTING },
 };
 #define NFTUPD asize(ftupd)
 
index ae28c7c775816f2df10eb251c332aaca252e7cdd..902453320fe0ae8e50c8342202cc65077adc6021 100644 (file)
@@ -3850,6 +3850,7 @@ const std::array<BondedInteractions, F_NRE> c_bondedInteractionFunctions
     BondedInteractions {unimplemented, -1 },                              // F_VSITE4FDN
     BondedInteractions {unimplemented, -1 },                              // F_VSITEN
     BondedInteractions {unimplemented, -1 },                              // F_COM_PULL
+    BondedInteractions {unimplemented, -1 },                              // F_DENSITYFITTING
     BondedInteractions {unimplemented, -1 },                              // F_EQM
     BondedInteractions {unimplemented, -1 },                              // F_EPOT
     BondedInteractions {unimplemented, -1 },                              // F_EKIN
index 5a8b55a1d040f26a10761d2fcad9df38bac7562f..6628445044847125f8d61b64613c427242c241da 100644 (file)
@@ -79,6 +79,7 @@
 #include "gromacs/trajectory/energyframe.h"
 #include "gromacs/utility/arraysize.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/mdmodulenotification.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/stringutil.h"
 
@@ -124,12 +125,13 @@ namespace gmx
  * \todo Remove GMX_CONSTRAINTVIR
  * \todo Write free-energy output also to energy file (after adding more tests)
  */
-EnergyOutput::EnergyOutput(ener_file        *fp_ene,
-                           const gmx_mtop_t *mtop,
-                           const t_inputrec *ir,
-                           const pull_t     *pull_work,
-                           FILE             *fp_dhdl,
-                           bool              isRerun)
+EnergyOutput::EnergyOutput(ener_file *              fp_ene,
+                           const gmx_mtop_t *       mtop,
+                           const t_inputrec *       ir,
+                           const pull_t *           pull_work,
+                           FILE *                   fp_dhdl,
+                           bool                     isRerun,
+                           const MdModulesNotifier &mdModulesNotifier)
 {
     const char                *ener_nm[F_NRE];
     static const char         *vir_nm[] = {
@@ -264,6 +266,12 @@ EnergyOutput::EnergyOutput(ener_file        *fp_ene,
     bEner_[F_ORIRESDEV]      = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
     bEner_[F_COM_PULL]       = ((ir->bPull && pull_have_potential(pull_work)) || ir->bRot);
 
+    MdModulesEnergyOutputToDensityFittingRequestChecker mdModulesAddOutputToDensityFittingFieldRequest;
+    mdModulesNotifier.notifier_.notify(&mdModulesAddOutputToDensityFittingFieldRequest);
+
+    bEner_[F_DENSITYFITTING] = mdModulesAddOutputToDensityFittingFieldRequest.energyOutputToDensityFitting_;
+
+
     // Counting the energy terms that will be printed and saving their names
     f_nre_ = 0;
     for (i = 0; i < F_NRE; i++)
index 0a6c8e4ecd7fc141310f5ccfe29b9f735543c3a3..4ecba82ce3f10cedf7d9102a5c44541555891615 100644 (file)
@@ -71,6 +71,7 @@ namespace gmx
 {
 class Awh;
 class Constraints;
+struct MdModulesNotifier;
 }
 
 //! \brief Printed names for intergroup energies
@@ -117,13 +118,15 @@ class EnergyOutput
          * \param[in] pull_work  Pulling simulations data
          * \param[in] fp_dhdl    FEP file.
          * \param[in] isRerun    Is this is a rerun instead of the simulations.
+         * \param[in] mdModulesNotifier Notifications to MD modules.
          */
-        EnergyOutput(ener_file        *fp_ene,
-                     const gmx_mtop_t *mtop,
-                     const t_inputrec *ir,
-                     const pull_t     *pull_work,
-                     FILE             *fp_dhdl,
-                     bool              isRerun);
+        EnergyOutput(ener_file *              fp_ene,
+                     const gmx_mtop_t *       mtop,
+                     const t_inputrec *       ir,
+                     const pull_t *           pull_work,
+                     FILE *                   fp_dhdl,
+                     bool                     isRerun,
+                     const MdModulesNotifier &mdModulesNotifier);
 
         ~EnergyOutput();
 
index a2f05a886af48a8661333a7499b04c914587f914..6ae860c0b1482c87e16314125b6c8b485f147fa2 100644 (file)
@@ -67,6 +67,7 @@
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/topology/topology.h"
 #include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/mdmodulenotification.h"
 #include "gromacs/utility/stringutil.h"
 #include "gromacs/utility/textreader.h"
 #include "gromacs/utility/unique_cptr.h"
@@ -541,7 +542,6 @@ class EnergyOutputTest : public ::testing::TestWithParam<EnergyOutputTestParamet
                 state_.nhpres_xi[k]  = (*testValue += 0.1);
                 state_.nhpres_vxi[k] = (*testValue += 0.1);
             }
-
         }
 
         /*! \brief Check if the contents of the .edr file correspond to the reference data.
@@ -623,7 +623,8 @@ TEST_P(EnergyOutputTest, CheckOutput)
         inputrec_.ref_p[YY][XX]   = 1.0;
     }
 
-    std::unique_ptr<EnergyOutput> energyOutput = std::make_unique<EnergyOutput>(energyFile_, &mtop_, &inputrec_, nullptr, nullptr, parameters.isRerun);
+    MdModulesNotifier             mdModulesNotifier;
+    std::unique_ptr<EnergyOutput> energyOutput = std::make_unique<EnergyOutput>(energyFile_, &mtop_, &inputrec_, nullptr, nullptr, parameters.isRerun, mdModulesNotifier);
 
     // Add synthetic data for a single step
     double testValue = 10.0;
index 9e33cc02ddf7b68e8005185f3e94344f863a7738..622887eaf455c4ea363c747ebb178b048fa3f02d 100644 (file)
@@ -257,7 +257,7 @@ void gmx::LegacySimulator::do_md()
     }
     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, ir, top_global, oenv, wcycle,
                                          StartingBehavior::NewSimulation);
-    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), false);
+    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), false, mdModulesNotifier);
 
     gstat = global_stat_init(ir);
 
index bdc0da71758bbe5c6b90d4e7100935e10780a03a..b10bed8bd03715b2601f758e97b3f3940cb9f581 100644 (file)
@@ -228,7 +228,7 @@ void gmx::LegacySimulator::do_mimic()
 
     gmx_mdoutf        *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, ir, top_global, oenv, wcycle,
                                           StartingBehavior::NewSimulation);
-    gmx::EnergyOutput  energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), true);
+    gmx::EnergyOutput  energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), true, mdModulesNotifier);
 
     gstat = global_stat_init(ir);
 
index c83f1658063f7e85cb6c89383cb47aa28a509785..31cccde41923411733479a704104eafb63bf3840 100644 (file)
@@ -1115,7 +1115,7 @@ LegacySimulator::do_cg()
             vsite, constr, nullptr);
     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec, top_global, nullptr, wcycle,
                                          StartingBehavior::NewSimulation);
-    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr, false);
+    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr, false, mdModulesNotifier);
 
     /* Print to log file */
     print_em_start(fplog, cr, walltime_accounting, wcycle, CG);
@@ -1754,7 +1754,7 @@ LegacySimulator::do_lbfgs()
             vsite, constr, nullptr);
     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec, top_global, nullptr, wcycle,
                                          StartingBehavior::NewSimulation);
-    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr, false);
+    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr, false, mdModulesNotifier);
 
     start = 0;
     end   = mdatoms->homenr;
@@ -2449,7 +2449,7 @@ LegacySimulator::do_steep()
             vsite, constr, nullptr);
     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec, top_global, nullptr, wcycle,
                                          StartingBehavior::NewSimulation);
-    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr, false);
+    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, inputrec, pull_work, nullptr, false, mdModulesNotifier);
 
     /* Print to log file  */
     print_em_start(fplog, cr, walltime_accounting, wcycle, SD);
index e722f39511cedc13bf049eba7692535715aa1e89..6c7519252bd3c05e6d90db00021dc3d489428cc4 100644 (file)
@@ -300,7 +300,7 @@ void gmx::LegacySimulator::do_rerun()
     initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
     gmx_mdoutf       *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, ir, top_global, oenv, wcycle,
                                          StartingBehavior::NewSimulation);
-    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), true);
+    gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), true, mdModulesNotifier);
 
     gstat = global_stat_init(ir);
 
index d35d731a15388b9b0683ea880fb63d4d7472e4f1..589988d17701637d9a41ce159b24f3046a050329 100644 (file)
@@ -64,16 +64,17 @@ namespace gmx
 class Awh;
 
 EnergyElement::EnergyElement(
-        StatePropagatorData *statePropagatorData,
-        const gmx_mtop_t    *globalTopology,
-        const t_inputrec    *inputrec,
-        const MDAtoms       *mdAtoms,
-        gmx_enerdata_t      *enerd,
-        gmx_ekindata_t      *ekind,
-        const Constraints   *constr,
-        FILE                *fplog,
-        t_fcdata            *fcd,
-        bool                 isMaster) :
+        StatePropagatorData     *statePropagatorData,
+        const gmx_mtop_t        *globalTopology,
+        const t_inputrec        *inputrec,
+        const MDAtoms           *mdAtoms,
+        gmx_enerdata_t          *enerd,
+        gmx_ekindata_t          *ekind,
+        const Constraints       *constr,
+        FILE                    *fplog,
+        t_fcdata                *fcd,
+        const MdModulesNotifier &mdModulesNotifier,
+        bool                     isMaster) :
     isMaster_(isMaster),
     energyWritingStep_(-1),
     energyCalculationStep_(-1),
@@ -94,6 +95,7 @@ EnergyElement::EnergyElement(
     constr_(constr),
     fplog_(fplog),
     fcd_(fcd),
+    mdModulesNotifier_(mdModulesNotifier),
     groups_(&globalTopology->groups)
 {
     clear_mat(forceVirial_);
@@ -144,7 +146,8 @@ void EnergyElement::trajectoryWriterSetup(gmx_mdoutf *outf)
     pull_t *pull_work = nullptr;
     energyOutput_ = std::make_unique<EnergyOutput>(
                 mdoutf_get_fp_ene(outf), top_global_,
-                inputrec_, pull_work, mdoutf_get_fp_dhdl(outf), false);
+                inputrec_, pull_work, mdoutf_get_fp_dhdl(outf), false,
+                mdModulesNotifier_);
 
     if (!isMaster_)
     {
index ebb3d8d40df0ff275ccf953a79b584a614262556..3819d72a6a8b409aa08ddc62399c5ca089752959 100644 (file)
@@ -59,6 +59,7 @@ class Constraints;
 class EnergyOutput;
 class MDAtoms;
 class StatePropagatorData;
+struct MdModulesNotifier;
 
 //! \addtogroup module_modularsimulator
 //! \{
@@ -102,6 +103,7 @@ class EnergyElement final :
             const Constraints *constr,
             FILE              *fplog,
             t_fcdata          *fcd,
+            const MdModulesNotifier &mdModulesNotifier,
             bool               isMaster);
 
         /*! \brief Register run function for step / time
@@ -273,25 +275,27 @@ class EnergyElement final :
          * Pointers to Simulator data
          */
         //! The microstate
-        StatePropagatorData   statePropagatorData_;
+        StatePropagatorData     *statePropagatorData_;
         //! Contains user input mdp options.
-        const t_inputrec       *inputrec_;
+        const t_inputrec        *inputrec_;
         //! Full system topology.
-        const gmx_mtop_t       *top_global_;
+        const gmx_mtop_t        *top_global_;
         //! Atom parameters for this domain.
-        const MDAtoms          *mdAtoms_;
+        const MDAtoms           *mdAtoms_;
         //! Energy data structure
-        gmx_enerdata_t         *enerd_;
+        gmx_enerdata_t          *enerd_;
         //! Kinetic energy data
-        gmx_ekindata_t         *ekind_;
+        gmx_ekindata_t          *ekind_;
         //! Handles constraints.
-        const Constraints      *constr_;
+        const Constraints       *constr_;
         //! Handles logging.
-        FILE                   *fplog_;
+        FILE                    *fplog_;
         //! Helper struct for force calculations.
-        t_fcdata               *fcd_;
+        t_fcdata                *fcd_;
+        //! Notification to MD modules
+        const MdModulesNotifier &mdModulesNotifier_;
         //! Global topology groups
-        const SimulationGroups *groups_;
+        const SimulationGroups  *groups_;
 };
 
 //! /}
index fed56d5c3112739bd3bda1fd6cbca60e4ab8de54..56346e70e9589e643cddcdbec47b82833ed99c89 100644 (file)
@@ -155,6 +155,7 @@ const t_interaction_function interaction_function[F_NRE] =
     def_vsite   ("VSITE4FDN", "Virtual site 4fdn", 5, 3                               ),
     def_vsite   ("VSITEN",   "Virtual site N",   2, 2                               ),
     def_nofc    ("COM_PULL", "COM Pull En."     ),
+    def_nofc    ("DENSITYFIT", "Density fitting"),
     def_nofc    ("EQM",      "Quantum En."      ),
     def_nofc    ("EPOT",     "Potential"        ),
     def_nofc    ("EKIN",     "Kinetic En."      ),
index b57e7368672c2634c89fddd3135256ca41348609..c3c1fffa43db9e02d26532bab51dfa89038cce1c 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2014,2015,2016,2018, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015,2016,2018,2019, 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.
@@ -117,7 +117,7 @@ struct t_interaction_function // NOLINT (clang-analyzer-optin.performance.Paddin
 #define IS_TABULATED(ftype) (interaction_function[(ftype)].flags & IF_TABULATED)
 
 /* this MUST correspond to the
-   t_interaction_function[F_NRE] in gmxlib/ifunc.c */
+   t_interaction_function[F_NRE] in gmxlib/ifunc.cpp */
 enum
 {
     F_BONDS,
@@ -194,6 +194,7 @@ enum
     F_VSITE4FDN,
     F_VSITEN,
     F_COM_PULL,
+    F_DENSITYFITTING,
     F_EQM,
     F_EPOT,
     F_EKIN,
index bec57d4a22fb4b103e45ec863f8973764e85a039..7f669aa32c7642ce0edb8262fefc48fa2c5f536c 100644 (file)
@@ -204,6 +204,11 @@ struct PeriodicBoundaryConditionType
     int pbcType;
 };
 
+struct MdModulesEnergyOutputToDensityFittingRequestChecker
+{
+    bool energyOutputToDensityFitting_ = false;
+};
+
 struct MdModulesNotifier
 {
 //! Register callback function types for MdModule
@@ -213,6 +218,7 @@ struct MdModulesNotifier
         KeyValueTreeObjectBuilder,
         const KeyValueTreeObject &,
         LocalAtomSetManager *,
+        MdModulesEnergyOutputToDensityFittingRequestChecker *,
         MdModulesCheckpointReadingDataOnMaster,
         MdModulesCheckpointReadingBroadcast,
         MdModulesWriteCheckpointData,