Feat/hybrid center xyz file

changes/user/enhancement.print-inter-region-center-xyz.md

@@ -0,0 +1 @@
+- Print the inner region center of hybrid type calculations (e.g. QM/MM) to a separate output file with the file ending .center.xyz

docs/sphinx/src/featureList/featureList.rst

@@ -18,7 +18,9 @@ The supported :ref:`jobtype <jobtype>` values are:
 * ``mm-md`` - molecular-mechanics molecular dynamics
 * ``qm-md`` - quantum-mechanics molecular dynamics
 * ``qm-rpmd`` - quantum-mechanics ring-polymer molecular dynamics
+* ``qmmm-md`` - hybrid quantum mechanics - molecular mechanics molecular dynamics
 * ``mm-opt`` - molecular-mechanics geometry optimization
+* ``mm-hessian`` - molecular mechanics Hessian calculation
 
 *******************
 Molecular Mechanics
@@ -113,7 +115,6 @@ Planned Items
 
 The following items are planned for future releases:
 
-* hybrid QM/MM job type
 * MM-RPMD
 * Verlet-list force evaluation
 * Ewald summation

docs/sphinx/src/userGuide/inputFile.rst

@@ -97,7 +97,7 @@ Possible options are:
 
    3. **qm-rpmd** - Represents a full quantum mechanics ring polymer molecular dynamics simulation. For more information see the :ref:`ringpolymermdKeys` section
 
-   4. **qmmm-md** - Represents a hybrid quantum mechanics - molecular mechanics molecular dynamics simulation. (Not implemented yet)
+   4. **qmmm-md** - Represents a hybrid quantum mechanics - molecular mechanics molecular dynamics simulation.
 
    5. **mm-opt** - represents a geometry optimization calculation using molecular mechanics.
 
@@ -476,6 +476,24 @@ The ``traj_file`` keyword sets the name for the :ref:`trajectoryFile`, which sto
 
 .. centered:: *default value* = "default.xyz"
 
+.. _hybridcenterfilekey:
+
+Hybrid Center File
+==================
+
+.. admonition:: Key
+    :class: tip
+
+    hybrid_center_file = {file} -> "default.center.xyz"
+
+The ``hybrid_center_file`` keyword sets the name for the :ref:`hybridCenterFile`, which stores the position of the inner region center in a hybrid type MD simulation.
+
+.. centered:: *default value* = "default.center.xyz"
+
+.. Note::
+
+    The ``hybrid_center_file`` keyword is ignored for non-hybrid-type calculations.
+
 .. _velocityfilekey:
 
 Velocity File

docs/sphinx/src/userGuide/outputFiles.rst

@@ -143,6 +143,24 @@ Hessian Info File
 
 Stores metadata for a generated Hessian, including the matrix filename, Hessian builder, optimization flag, finite-difference displacement, sign convention, unit and matrix dimensions. The file starts with ``format = pq-hessian-info-v1`` so downstream tools can identify the metadata format.
 
+.. _hybridCenterFile:
+
+******************
+Hybrid Center File
+******************
+
+**File Type:** ``.center.xyz``
+
+Stores the coordinates (*x*, *y*, *z*) of the inner region center in ``.xyz`` file format:
+
+    | **line 1:** 1
+    | **line 2:** empty or "# step = step_number"
+    | **line 3:** X *x* *y* *z*
+
+There is one dummy atom used to denote the position of the inner region center.
+The atom has the name ``X`` and its coordinates are denoted by *x* *y* *z* in Å.
+The file is only written for :ref:`hybrid type calculations <HybridCalculationKeys>`.
+
 .. _instantEnergyFile:
 
 ********************

include/config/defaults.hpp

@@ -32,23 +32,24 @@
  */
 struct DefaultFiles
 {
-    static constexpr auto restartFile  = "default.rst";
-    static constexpr auto energyFile   = "default.en";
-    static constexpr auto instEnFile   = "default.instant_en";
-    static constexpr auto momentumFile = "default.mom";
-    static constexpr auto trajFile     = "default.xyz";
-    static constexpr auto velFile      = "default.vel";
-    static constexpr auto forceFile    = "default.force";
-    static constexpr auto chargeFile   = "default.charge";
-    static constexpr auto logFile      = "default.log";
-    static constexpr auto stdoutFile   = "default.stdout";
-    static constexpr auto refFile      = "default.ref";
-    static constexpr auto infoFile     = "default.info";
-    static constexpr auto virialFile   = "default.vir";
-    static constexpr auto stressFile   = "default.stress";
-    static constexpr auto boxFile      = "default.box";
-    static constexpr auto optFile      = "default.opt";
-    static constexpr auto timingsFile  = "default.timings";
+    static constexpr auto restartFile      = "default.rst";
+    static constexpr auto energyFile       = "default.en";
+    static constexpr auto instEnFile       = "default.instant_en";
+    static constexpr auto momentumFile     = "default.mom";
+    static constexpr auto trajFile         = "default.xyz";
+    static constexpr auto hybridCenterFile = "default.center.xyz";
+    static constexpr auto velFile          = "default.vel";
+    static constexpr auto forceFile        = "default.force";
+    static constexpr auto chargeFile       = "default.charge";
+    static constexpr auto logFile          = "default.log";
+    static constexpr auto stdoutFile       = "default.stdout";
+    static constexpr auto refFile          = "default.ref";
+    static constexpr auto infoFile         = "default.info";
+    static constexpr auto virialFile       = "default.vir";
+    static constexpr auto stressFile       = "default.stress";
+    static constexpr auto boxFile          = "default.box";
+    static constexpr auto optFile          = "default.opt";
+    static constexpr auto timingsFile      = "default.timings";
 
     static constexpr auto rpmdRstFile    = "default.rpmd.rst";
     static constexpr auto rpmdTrajFile   = "default.rpmd.xyz";
@@ -83,6 +84,8 @@ namespace defaults
     static constexpr bool   HESSIAN_OPTIMIZE_DEFAULT     = true;
     static constexpr auto*  HESSIAN_BUILDER_DEFAULT      = "central";
 
+    static constexpr char INNER_REGION_CENTER_ATOM_NAME = 'X';
+
     static constexpr double COULOMB_CUT_OFF_DEFAULT  = 12.5;   // in Angstrom
     static constexpr double SCALE_14_COULOMB_DEFAULT = 1.0;
     static constexpr double SCALE_14_VAN_DER_WAALS_DEFAULT = 1.0;

include/engine/engineOutput.hpp

@@ -45,6 +45,11 @@
 #include "trajectoryOutput.hpp"
 #include "virialOutput.hpp"
 
+namespace configurator
+{
+    class HybridConfigurator;   // forward declaration
+}   // namespace configurator
+
 namespace engine
 {
     /**
@@ -61,6 +66,7 @@ namespace engine
         std::unique_ptr<output::InfoOutput>   _infoOutput;
 
         std::unique_ptr<output::TrajectoryOutput> _xyzOutput;
+        std::unique_ptr<output::TrajectoryOutput> _xyzHybridCenterOutput;
         std::unique_ptr<output::TrajectoryOutput> _velOutput;
         std::unique_ptr<output::TrajectoryOutput> _forceOutput;
         std::unique_ptr<output::TrajectoryOutput> _chargeOutput;
@@ -97,7 +103,12 @@ namespace engine
             const size_t step,
             const physicalData::PhysicalData &
         );
+
         void writeXyzFile(simulationBox::SimulationBox &, const size_t);
+        void writeHybridCenterXyzFile(
+            const configurator::HybridConfigurator &,
+            const size_t
+        );
         void writeVelFile(simulationBox::SimulationBox &, const size_t);
         void writeForceFile(simulationBox::SimulationBox &, const size_t);
         void writeChargeFile(simulationBox::SimulationBox &, const size_t);
@@ -154,6 +165,7 @@ namespace engine
         [[nodiscard]] output::EnergyOutput     &getEnergyOutput();
         [[nodiscard]] output::EnergyOutput     &getInstantEnergyOutput();
         [[nodiscard]] output::TrajectoryOutput &getXyzOutput();
+        [[nodiscard]] output::TrajectoryOutput &getXyzHybridCenterOutput();
         [[nodiscard]] output::TrajectoryOutput &getVelOutput();
         [[nodiscard]] output::TrajectoryOutput &getForceOutput();
         [[nodiscard]] output::TrajectoryOutput &getChargeOutput();

include/engine/hybridMDEngine.hpp

@@ -24,7 +24,6 @@
 
 #define _HYBRID_MD_ENGINE_HPP_
 
-#include "hybridConfigurator.hpp"
 #include "mdEngine.hpp"
 #include "qmCapableEngine.hpp"
 
@@ -41,8 +40,6 @@ namespace engine
     class HybridMDEngine : virtual public MDEngine, public QMCapableEngine
     {
        protected:
-        configurator::HybridConfigurator _configurator{};
-
         void combineInnerOuterForces();
 
         void addCurrentForcesToInnerAndReset(

include/engine/mdEngine.hpp

@@ -49,7 +49,8 @@ namespace engine
     class MDEngine : public Engine
     {
        protected:
-        resetKinetics::ResetKinetics _resetKinetics;
+        resetKinetics::ResetKinetics     _resetKinetics;
+        configurator::HybridConfigurator _configurator{};
 
         std::unique_ptr<integrator::Integrator> _integrator;
         std::unique_ptr<thermostat::Thermostat> _thermostat;
@@ -79,6 +80,7 @@ namespace engine
         [[nodiscard]] manostat::Manostat           &getManostat();
         [[nodiscard]] output::EnergyOutput         &getInstantEnergyOutput();
         [[nodiscard]] output::MomentumOutput       &getMomentumOutput();
+        [[nodiscard]] output::TrajectoryOutput     &getXyzHybridCenterOutput();
         [[nodiscard]] output::TrajectoryOutput     &getVelOutput();
         [[nodiscard]] output::TrajectoryOutput     &getChargeOutput();
         [[nodiscard]] output::VirialOutput         &getVirialOutput();

include/input/inputFileParser/outputInputParser.hpp

@@ -76,6 +76,11 @@ namespace input
             const size_t
         );
 
+        void parseHybridCenterFilename(
+            const std::vector<std::string> &,
+            const size_t
+        );
+
         void parseVelocityFilename(
             const std::vector<std::string> &,
             const size_t

include/output/trajectoryOutput.hpp

@@ -26,13 +26,19 @@
 
 #include <cstddef>   // for size_t
 
+#include "hybridConfigurator.hpp"
 #include "output.hpp"   // for Output
 
 namespace simulationBox
 {
     class SimulationBox;   // forward declaration
 }   // namespace simulationBox
 
+namespace configurator
+{
+    class HybridConfigurator;   // forward declaration
+}   // namespace configurator
+
 namespace output
 {
     /**
@@ -48,6 +54,10 @@ namespace output
 
         void writeHeader(const simulationBox::SimulationBox &);
         void writeXyz(simulationBox::SimulationBox &, const size_t);
+        void writeHybridCenterXyz(
+            const configurator::HybridConfigurator &configurator,
+            const size_t                            step
+        );
         void writeVelocities(simulationBox::SimulationBox &, const size_t);
         void writeForces(simulationBox::SimulationBox &, const size_t);
         void writeCharges(simulationBox::SimulationBox &, const size_t);

include/settings/outputFileSettings.hpp

@@ -48,17 +48,20 @@ namespace settings
         static inline bool        _filePrefixSet         = false;
         static inline std::string _filePrefix;
 
-        static inline std::string _energyFile = DefaultFiles::energyFile;
-        static inline std::string _instEnFile = DefaultFiles::instEnFile;
-        static inline std::string _rstFile    = DefaultFiles::restartFile;
-        static inline std::string _momFile    = DefaultFiles::momentumFile;
-        static inline std::string _trajFile   = DefaultFiles::trajFile;
-        static inline std::string _velFile    = DefaultFiles::velFile;
-        static inline std::string _forceFile  = DefaultFiles::forceFile;
-        static inline std::string _chargeFile = DefaultFiles::chargeFile;
-        static inline std::string _logFile    = DefaultFiles::logFile;
-        static inline std::string _refFile    = DefaultFiles::refFile;
-        static inline std::string _infoFile   = DefaultFiles::infoFile;
+        // clang-format off
+        static inline std::string _energyFile       = DefaultFiles::energyFile;
+        static inline std::string _instEnFile       = DefaultFiles::instEnFile;
+        static inline std::string _rstFile          = DefaultFiles::restartFile;
+        static inline std::string _momFile          = DefaultFiles::momentumFile;
+        static inline std::string _trajFile         = DefaultFiles::trajFile;
+        static inline std::string _hybridCenterFile = DefaultFiles::hybridCenterFile;
+        static inline std::string _velFile          = DefaultFiles::velFile;
+        static inline std::string _forceFile        = DefaultFiles::forceFile;
+        static inline std::string _chargeFile       = DefaultFiles::chargeFile;
+        static inline std::string _logFile          = DefaultFiles::logFile;
+        static inline std::string _refFile          = DefaultFiles::refFile;
+        static inline std::string _infoFile         = DefaultFiles::infoFile;
+        // clang-format on
 
         static inline std::string _virialFile = DefaultFiles::virialFile;
         static inline std::string _stressFile = DefaultFiles::stressFile;
@@ -96,6 +99,7 @@ namespace settings
         static void setInstantEnergyFileName(const std::string_view);
         static void setMomentumFileName(const std::string_view);
         static void setTrajectoryFileName(const std::string_view);
+        static void setHybridCenterFileName(const std::string_view);
         static void setVelocityFileName(const std::string_view);
         static void setForceFileName(const std::string_view);
         static void setChargeFileName(const std::string_view);
@@ -135,6 +139,7 @@ namespace settings
         [[nodiscard]] static std::string getInstantEnergyFileName();
         [[nodiscard]] static std::string getMomentumFileName();
         [[nodiscard]] static std::string getTrajectoryFileName();
+        [[nodiscard]] static std::string getHybridCenterFileName();
         [[nodiscard]] static std::string getVelocityFileName();
         [[nodiscard]] static std::string getForceFileName();
         [[nodiscard]] static std::string getChargeFileName();

src/engine/engineOutput.cpp

@@ -27,6 +27,7 @@
 #include "boxOutput.hpp"
 #include "defaults.hpp"
 #include "energyOutput.hpp"
+#include "hybridConfigurator.hpp"
 #include "infoOutput.hpp"
 #include "logOutput.hpp"
 #include "momentumOutput.hpp"
@@ -55,6 +56,7 @@ using namespace engine;
 using namespace simulationBox;
 using namespace physicalData;
 using namespace thermostat;
+using namespace configurator;
 
 using output::BoxFileOutput;
 using output::EnergyOutput;
@@ -81,6 +83,7 @@ EngineOutput::EngineOutput()
       _instantEnergyOutput(std::make_unique<EnergyOutput>(DefaultFiles::instEnFile)),
       _infoOutput(std::make_unique<InfoOutput>(DefaultFiles::infoFile)),
       _xyzOutput(std::make_unique<TrajectoryOutput>(DefaultFiles::trajFile)),
+      _xyzHybridCenterOutput(std::make_unique<TrajectoryOutput>(DefaultFiles::hybridCenterFile)),
       _velOutput(std::make_unique<TrajectoryOutput>(DefaultFiles::velFile)),
       _forceOutput(std::make_unique<TrajectoryOutput>(DefaultFiles::forceFile)),
       _chargeOutput(std::make_unique<TrajectoryOutput>(DefaultFiles::chargeFile)),
@@ -161,6 +164,21 @@ void EngineOutput::writeXyzFile(SimulationBox &simulationBox, const size_t step)
     _xyzOutput->writeXyz(simulationBox, step);
 }
 
+/**
+ * @brief wrapper for hybrid center xyz file output function
+ *
+ * @param configurator
+ * @param step
+ */
+void EngineOutput::writeHybridCenterXyzFile(
+    const HybridConfigurator &configurator,
+    const size_t              step
+)
+{
+    auto _ = scoped("TrajectoryOutput");
+    _xyzHybridCenterOutput->writeHybridCenterXyz(configurator, step);
+}
+
 /**
  * @brief wrapper for velocity file output function
  *
@@ -443,6 +461,16 @@ MomentumOutput &EngineOutput::getMomentumOutput() { return *_momentumOutput; }
  */
 TrajectoryOutput &EngineOutput::getXyzOutput() { return *_xyzOutput; }
 
+/**
+ * @brief getter for hybrid center xyz output
+ *
+ * @return TrajectoryOutput
+ */
+TrajectoryOutput &EngineOutput::getXyzHybridCenterOutput()
+{
+    return *_xyzHybridCenterOutput;
+}
+
 /**
  * @brief getter for velocity output
  *

src/engine/mdEngine.cpp

@@ -251,6 +251,9 @@ void MDEngine::writeOutput()
         );   // use physicalData instead of averagePhysicalData
 
         _engineOutput.writeBoxFile(effStep, _simulationBox->getBox());
+
+        if (Settings::isHybridJobtype())
+            _engineOutput.writeHybridCenterXyzFile(_configurator, effStep);
     }
 
     // NOTE:
@@ -335,6 +338,16 @@ output::MomentumOutput &MDEngine::getMomentumOutput()
     return _engineOutput.getMomentumOutput();
 }
 
+/**
+ * @brief get the reference to the xyz hybrid center output
+ *
+ * @return output::TrajectoryOutput&
+ */
+output::TrajectoryOutput &MDEngine::getXyzHybridCenterOutput()
+{
+    return _engineOutput.getXyzHybridCenterOutput();
+}
+
 /**
  * @brief get the reference to the vel output
  *