feat: drivenMD setup and testing

src/module/module.cpp

@@ -26,6 +26,7 @@ EnumOptions<ModuleTypes::ModuleType> moduleTypes_("ModuleType", {{ModuleTypes::A
                                                                  {ModuleTypes::Compare, "Compare"},
                                                                  {ModuleTypes::DAngle, "DAngle"},
                                                                  {ModuleTypes::DataTest, "DataTest"},
+                                                                 {ModuleTypes::DrivenMD, "DrivenMD"},
                                                                  {ModuleTypes::Energy, "Energy"},
                                                                  {ModuleTypes::EPSR, "EPSR"},
                                                                  {ModuleTypes::EPSRManager, "EPSRManager"},

src/module/types.h

@@ -20,6 +20,7 @@ enum ModuleType
     Compare,
     DAngle,
     DataTest,
+    DrivenMD,
     Energy,
     EPSR,
     EPSRManager,

src/modules/drivenMD/CMakeLists.txt

@@ -0,0 +1 @@
+dissolve_add_module(drivenMD.h drivenMD)

src/modules/drivenMD/drivenMD.cpp

@@ -0,0 +1,18 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "modules/drivenMD/drivenMD.h"
+#include "keywords/configuration.h"
+#include "keywords/fileAndFormat.h"
+#include "keywords/moduleVector.h"
+#include "keywords/speciesSiteVector.h"
+#include "keywords/weightedModuleVector.h"
+
+DrivenMDModule::DrivenMDModule() : Module(ModuleTypes::DrivenMD)
+{
+    keywords_.addTarget<ConfigurationKeyword>("Configuration", "Set target configuration for the module", targetConfiguration_)
+        ->setEditSignals({KeywordBase::ClearModuleData, KeywordBase::RecreateRenderables});
+    keywords_.addTarget<ModuleVectorKeyword>("Target", "Add specified Module (and it's Reference data) as a refinement target",
+                                             targets_,
+                                             std::vector<ModuleTypes::ModuleType>{ModuleTypes::NeutronSQ, ModuleTypes::XRaySQ});
+}

src/modules/drivenMD/drivenMD.h

@@ -0,0 +1,52 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#pragma once
+
+#include "base/enumOptions.h"
+#include "classes/scatteringMatrix.h"
+#include "math/data1D.h"
+#include "math/range.h"
+#include "module/groups.h"
+#include "module/module.h"
+#include "templates/array3D.h"
+#include <tuple>
+
+// DrivenMD Module
+class DrivenMDModule : public Module
+{
+    public:
+    DrivenMDModule();
+    ~DrivenMDModule() override = default;
+
+    /*
+     * Definition
+     */
+    public:
+    private:
+    // Vector storing atom pairs and associated potentials
+    std::vector<std::tuple<std::shared_ptr<AtomType>, std::shared_ptr<AtomType>, Data1D>> empiricalPotentials_;
+    // Target Modules containing data to refine against
+    std::vector<Module *> targets_;
+
+    public:
+    // Return list of target Modules / data for refinement
+    const std::vector<Module *> &targets() const;
+    // Set whether to apply this module's generated potentials to the global pair potentials
+    void setApplyPotentials(bool b);
+
+    /*
+     * Functions
+     */
+    private:
+    // Target Configuration (determined from target modules)
+    Configuration *targetConfiguration_{nullptr};
+    // Calculate partial g(r) in serial with simple double-loop
+    void calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet);
+    /*
+     * Processing
+     */
+    private:
+    // Run main processing
+    Module::ExecutionResult process(ModuleContext &moduleContext) override;
+};

src/modules/drivenMD/functions.cpp

@@ -0,0 +1,28 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "main/dissolve.h"
+#include "module/context.h"
+#include "modules/drivenMD/drivenMD.h"
+#include "templates/algorithms.h"
+
+// Calculate partial g(r) in serial with simple double-loop
+void DrivenMDModule::calculateGRTestSerial(Configuration *cfg, PartialSet &partialSet)
+{
+    // Calculate radial distribution functions with a simple double loop, in serial then return the partial set!!!!!!!
+    const auto *box = cfg->box();
+
+    dissolve::for_each_pair(
+        ParallelPolicies::seq, cfg->atoms().begin(), cfg->atoms().end(),
+        [box, &partialSet](auto i, auto &ii, auto j, auto &jj)
+        {
+            if (&ii != &jj)
+                partialSet.fullHistogram(ii.localTypeIndex(), jj.localTypeIndex()).bin(box->minimumDistance(ii.r(), jj.r()));
+        });
+
+    // Transform histogram data into radial distribution functions
+    partialSet.formPartials(box->volume());
+
+    // Sum total functions
+    partialSet.formTotals(true);
+}

src/modules/drivenMD/gui/CMakeLists.txt

@@ -0,0 +1 @@
+dissolve_add_module_gui(drivenMD)

src/modules/drivenMD/gui/drivenMDFuncs.cpp

@@ -0,0 +1,15 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "main/dissolve.h"
+#include "modules/drivenMD/drivenMD.h"
+#include "modules/drivenMD/gui/drivenMDWidget.h"
+
+DrivenMDModuleWidget::DrivenMDModuleWidget(QWidget *parent, DrivenMDModule *module, Dissolve &dissolve)
+    : ModuleWidget(parent, dissolve), module_(module)
+{
+    // Set up user interface
+    ui_.setupUi(this);
+
+    refreshing_ = false;
+}

src/modules/drivenMD/gui/drivenMDWidget.h

@@ -0,0 +1,34 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#pragma once
+
+#include "modules/drivenMD/gui/ui_drivenMDWidget.h"
+#include "modules/widget.h"
+
+// Forward Declarations
+class DrivenMDModule;
+class Dissolve;
+class Module;
+
+// Module Widget
+class DrivenMDModuleWidget : public ModuleWidget
+{
+    // All Qt declarations derived from QObject must include this macro
+    Q_OBJECT
+
+    public:
+    DrivenMDModuleWidget(QWidget *parent, DrivenMDModule *module, Dissolve &dissolve);
+    ~DrivenMDModuleWidget() override = default;
+
+    private:
+    // Associated Module
+    DrivenMDModule *module_;
+
+    /*
+     * UI
+     */
+    private:
+    // Main form declaration
+    Ui::DrivenMDModuleWidget ui_;
+};

src/modules/drivenMD/gui/drivenMDWidget.ui

@@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
+ <class>DrivenMDModuleWidget</class>
+ <widget class="QWidget" name="DrivenMDModuleWidget">
+  <property name="geometry">
+   <rect>
+    <x>0</x>
+    <y>0</y>
+    <width>473</width>
+    <height>395</height>
+   </rect>
+  </property>
+  <property name="sizePolicy">
+   <sizepolicy hsizetype="Expanding" vsizetype="Expanding">
+    <horstretch>0</horstretch>
+    <verstretch>0</verstretch>
+   </sizepolicy>
+  </property>
+  <property name="font">
+   <font>
+    <pointsize>8</pointsize>
+   </font>
+  </property>
+  <property name="windowTitle">
+   <string>AnalyseModule Controls</string>
+  </property>
+  <layout class="QVBoxLayout" name="verticalLayout">
+   <property name="spacing">
+    <number>4</number>
+   </property>
+   <property name="leftMargin">
+    <number>0</number>
+   </property>
+   <property name="topMargin">
+    <number>0</number>
+   </property>
+   <property name="rightMargin">
+    <number>0</number>
+   </property>
+   <property name="bottomMargin">
+    <number>0</number>
+   </property>
+  </layout>
+ </widget>
+ <resources/>
+ <connections/>
+</ui>

src/modules/drivenMD/process.cpp

@@ -0,0 +1,159 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "base/sysFunc.h"
+#include "keywords/module.h"
+#include "main/dissolve.h"
+#include "math/error.h"
+#include "math/ft.h"
+#include "module/context.h"
+#include "module/group.h"
+#include "modules/drivenMD/drivenMD.h"
+#include "modules/gr/gr.h"
+#include <functional>
+
+// Run main processing
+Module::ExecutionResult DrivenMDModule::process(ModuleContext &moduleContext)
+{
+    auto &processingData = moduleContext.dissolve().processingModuleData();
+
+    /*
+     * 1) Loop over a site
+     * 2) Move in the x, y and z direction
+     * 3) Calculate gr at each step
+     * 4) FT gr-1 to F(Q) and weight by scattering length
+     * 5) Take the negative gradient of the difference between the gr and the F(Q) which corresponds to the derived force
+     * that we then apply to MD simulation
+     */
+
+    std::vector<Vec3<double>> forces;
+
+    // Does a PartialSet already exist for this Configuration?
+    auto originalGRObject = processingData.realiseIf<PartialSet>(
+        std::format("{}//OriginalGR", targetConfiguration_->niceName()), name_, GenericItem::InRestartFileFlag);
+    auto &originalgr = originalGRObject.first;
+    auto rdfRange = targetConfiguration_->box()->inscribedSphereRadius();
+    auto binWidth = 0.025;
+    rdfRange = int(rdfRange / binWidth) * binWidth;
+    if (originalGRObject.second == GenericItem::ItemStatus::Created)
+        originalgr.setUp(targetConfiguration_->atomTypePopulations(), rdfRange, binWidth);
+
+    originalgr.setUpHistograms(rdfRange, binWidth);
+    originalgr.reset();
+
+    double delta{2};
+    auto atoms = targetConfiguration_->atoms();
+    forces.reserve(atoms.size());
+
+    for (auto *module : targets_)
+    {
+        // Retrieve the ReferenceData
+        if (!processingData.contains("ReferenceData", module->name()))
+        {
+            Messenger::error("Reference data not found for target '{}'.\n", module->name());
+            return ExecutionResult::Failed;
+        }
+        const auto &originalReferenceData = processingData.value<Data1D>("ReferenceData", module->name());
+
+        for (auto &i : atoms)
+        {
+            std::format("X coord = {}, Y coord = {}, Z coord = {}\n", i.x(), i.y(), i.z());
+            auto &f = forces.emplace_back(Vec3<double>(0.0, 0.0, 0.0));
+            std::format("Atom {}\n", i.globalIndex());
+            for (auto n = 0; n < 3; ++n)
+            {
+                switch (n)
+                {
+                    // Move x
+                    case 0:
+                        // Get position, change via delta then set
+                        std::format("X coord = {}\n", i.x());
+                        i.translateCoordinates(-delta, 0, 0);
+                        std::format("X coord = {}\n", i.x());
+                        // Reset GR
+                        originalgr.reset();
+                        // Calculate GR
+                        calculateGRTestSerial(targetConfiguration_, originalgr);
+                        // FT to structure factor
+                        Fourier::sineFT(originalgr.total(), 1.0 / (2.0 * PI * PI * 1.39), 0.05, 0.05, 30.0,
+                                        WindowFunction::Form::Lorch0);
+                        // Store the error in a vector
+                        f.x = Error::rmse(originalgr.total(), originalReferenceData).error;
+                        std::format("X+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error);
+                        // Set new position (needs 2x to get +x from original)
+                        i.translateCoordinates(2 * delta, 0, 0);
+                        std::format("X coord = {}\n", i.x());
+                        originalgr.reset();
+                        calculateGRTestSerial(targetConfiguration_, originalgr);
+                        Fourier::sineFT(originalgr.total(), 1.0 / (2.0 * PI * PI * 1.39), 0.05, 0.05, 30.0,
+                                        WindowFunction::Form::Lorch0);
+                        // Calculate new error
+                        f.x -= Error::rmse(originalgr.total(), originalReferenceData).error;
+                        std::format("X+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error);
+                        // Reset position
+                        i.translateCoordinates(-delta, 0, 0);
+                        std::format("X coord = {}\n", i.x());
+                        break;
+                    // Move y and repeat
+                    case 1:
+                        std::format("Y coord = {}\n", i.y());
+                        i.translateCoordinates(0, -delta, 0);
+                        std::format("Y coord = {}\n", i.y());
+                        originalgr.reset();
+                        calculateGRTestSerial(targetConfiguration_, originalgr);
+                        Fourier::sineFT(originalgr.total(), 1.0 / (2.0 * PI * PI * 1.39), 0.05, 0.05, 30.0,
+                                        WindowFunction::Form::Lorch0);
+                        f.y = Error::rmse(originalgr.total(), originalReferenceData).error;
+                        std::format("Y+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error);
+                        i.translateCoordinates(0, 2 * delta, 0);
+                        std::format("Y coord = {}\n", i.y());
+                        originalgr.reset();
+                        calculateGRTestSerial(targetConfiguration_, originalgr);
+                        Fourier::sineFT(originalgr.total(), 1.0 / (2.0 * PI * PI * 1.39), 0.05, 0.05, 30.0,
+                                        WindowFunction::Form::Lorch0);
+                        f.y -= Error::rmse(originalgr.total(), originalReferenceData).error;
+                        std::format("Y+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error);
+                        i.translateCoordinates(0, -delta, 0);
+                        std::format("Y coord = {}\n", i.y());
+                        break;
+                    // Move z
+                    case 2:
+                        std::format("Z coord = {}\n", i.z());
+                        i.translateCoordinates(0, 0, -delta);
+                        std::format("Z coord = {}\n", i.z());
+                        originalgr.reset();
+                        calculateGRTestSerial(targetConfiguration_, originalgr);
+                        Fourier::sineFT(originalgr.total(), 1.0 / (2.0 * PI * PI * 1.39), 0.05, 0.05, 30.0,
+                                        WindowFunction::Form::Lorch0);
+                        f.z = Error::rmse(originalgr.total(), originalReferenceData).error;
+                        std::format("Z+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error);
+                        i.translateCoordinates(0, 0, 2 * delta);
+                        std::format("Z coord = {}\n", i.z());
+                        originalgr.reset();
+                        calculateGRTestSerial(targetConfiguration_, originalgr);
+                        Fourier::sineFT(originalgr.total(), 1.0 / (2.0 * PI * PI * 1.39), 0.05, 0.05, 30.0,
+                                        WindowFunction::Form::Lorch0);
+                        f.z -= Error::rmse(originalgr.total(), originalReferenceData).error;
+                        std::format("Z+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error);
+                        i.translateCoordinates(0, 0, -delta);
+                        std::format("Z coord = {}\n", i.z());
+                        break;
+                }
+            }
+        }
+    }
+
+    for (auto &&[i, f] : zip(atoms, forces))
+    {
+        std::format("X coord = {}, Y coord = {}, Z coord = {}\n", i.x(), i.y(), i.z());
+        i.translateCoordinates(f);
+        f.print();
+        std::format("We've moved\n");
+        std::format("X coord = {}, Y coord = {}, Z coord = {}\n", i.x(), i.y(), i.z());
+    }
+
+    targetConfiguration_->incrementContentsVersion();
+    targetConfiguration_->updateObjectRelationships();
+
+    return ExecutionResult::Success;
+}

src/modules/registry.cpp

@@ -11,6 +11,7 @@
 #include "modules/bragg/bragg.h"
 #include "modules/compare/compare.h"
 #include "modules/dAngle/dAngle.h"
+#include "modules/drivenMD/drivenMD.h"
 #include "modules/energy/energy.h"
 #include "modules/epsr/epsr.h"
 #include "modules/epsrManager/epsrManager.h"
@@ -54,6 +55,7 @@ ModuleRegistry::ModuleRegistry()
     registerProducer<BraggModule>(ModuleTypes::Bragg, "Calculate Bragg intensities", "Correlation Functions");
     registerProducer<CompareModule>(ModuleTypes::Compare, "Compare data sets and calculate errors", "Checks & Tests");
     registerProducer<DAngleModule>(ModuleTypes::DAngle, "Calculate distance/angle maps", "Analysis");
+    registerProducer<DrivenMDModule>(ModuleTypes::DrivenMD, "SAXS/SANS driven molecular dynamics test", "Forcefield");
     registerProducer<EnergyModule>(ModuleTypes::Energy, "Calculate the total energy of a Configuration", "Forcefield");
     registerProducer<EPSRModule>(ModuleTypes::EPSR, "Refine interatomic potentials in a manner consistent with EPSR",
                                  "Forcefield");