diff --git a/src/module/module.cpp b/src/module/module.cpp index d4953db1cd..387cda17b2 100644 --- a/src/module/module.cpp +++ b/src/module/module.cpp @@ -26,6 +26,7 @@ EnumOptions 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"}, diff --git a/src/module/types.h b/src/module/types.h index 2ce1d453b6..8c1a4c9c0f 100644 --- a/src/module/types.h +++ b/src/module/types.h @@ -20,6 +20,7 @@ enum ModuleType Compare, DAngle, DataTest, + DrivenMD, Energy, EPSR, EPSRManager, diff --git a/src/modules/drivenMD/CMakeLists.txt b/src/modules/drivenMD/CMakeLists.txt new file mode 100644 index 0000000000..ad792d62d1 --- /dev/null +++ b/src/modules/drivenMD/CMakeLists.txt @@ -0,0 +1 @@ +dissolve_add_module(drivenMD.h drivenMD) diff --git a/src/modules/drivenMD/drivenMD.cpp b/src/modules/drivenMD/drivenMD.cpp new file mode 100644 index 0000000000..015122eac4 --- /dev/null +++ b/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("Configuration", "Set target configuration for the module", targetConfiguration_) + ->setEditSignals({KeywordBase::ClearModuleData, KeywordBase::RecreateRenderables}); + keywords_.addTarget("Target", "Add specified Module (and it's Reference data) as a refinement target", + targets_, + std::vector{ModuleTypes::NeutronSQ, ModuleTypes::XRaySQ}); +} diff --git a/src/modules/drivenMD/drivenMD.h b/src/modules/drivenMD/drivenMD.h new file mode 100644 index 0000000000..c9fc9f1b6f --- /dev/null +++ b/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 + +// DrivenMD Module +class DrivenMDModule : public Module +{ + public: + DrivenMDModule(); + ~DrivenMDModule() override = default; + + /* + * Definition + */ + public: + private: + // Vector storing atom pairs and associated potentials + std::vector, std::shared_ptr, Data1D>> empiricalPotentials_; + // Target Modules containing data to refine against + std::vector targets_; + + public: + // Return list of target Modules / data for refinement + const std::vector &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; +}; diff --git a/src/modules/drivenMD/functions.cpp b/src/modules/drivenMD/functions.cpp new file mode 100644 index 0000000000..a75a4c1c39 --- /dev/null +++ b/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); +} diff --git a/src/modules/drivenMD/gui/CMakeLists.txt b/src/modules/drivenMD/gui/CMakeLists.txt new file mode 100644 index 0000000000..96a541db4b --- /dev/null +++ b/src/modules/drivenMD/gui/CMakeLists.txt @@ -0,0 +1 @@ +dissolve_add_module_gui(drivenMD) diff --git a/src/modules/drivenMD/gui/drivenMDFuncs.cpp b/src/modules/drivenMD/gui/drivenMDFuncs.cpp new file mode 100644 index 0000000000..bd5e8c8c12 --- /dev/null +++ b/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; +} diff --git a/src/modules/drivenMD/gui/drivenMDWidget.h b/src/modules/drivenMD/gui/drivenMDWidget.h new file mode 100644 index 0000000000..44d03e2c73 --- /dev/null +++ b/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_; +}; diff --git a/src/modules/drivenMD/gui/drivenMDWidget.ui b/src/modules/drivenMD/gui/drivenMDWidget.ui new file mode 100644 index 0000000000..159430921f --- /dev/null +++ b/src/modules/drivenMD/gui/drivenMDWidget.ui @@ -0,0 +1,47 @@ + + + DrivenMDModuleWidget + + + + 0 + 0 + 473 + 395 + + + + + 0 + 0 + + + + + 8 + + + + AnalyseModule Controls + + + + 4 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + diff --git a/src/modules/drivenMD/process.cpp b/src/modules/drivenMD/process.cpp new file mode 100644 index 0000000000..b0e27817ef --- /dev/null +++ b/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 + +// 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> forces; + + // Does a PartialSet already exist for this Configuration? + auto originalGRObject = processingData.realiseIf( + fmt::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_->atomTypes(), 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("ReferenceData", module->name()); + + for (auto &i : atoms) + { + fmt::print("X coord = {}, Y coord = {}, Z coord = {}\n", i.x(), i.y(), i.z()); + auto &f = forces.emplace_back(Vec3(0.0, 0.0, 0.0)); + fmt::print("Atom {}\n", i.globalIndex()); + for (auto n = 0; n < 3; ++n) + { + switch (n) + { + // Move x + case 0: + // Get position, change via delta then set + fmt::print("X coord = {}\n", i.x()); + i.translateCoordinates(-delta, 0, 0); + fmt::print("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; + fmt::print("X+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error); + // Set new position (needs 2x to get +x from original) + i.translateCoordinates(2 * delta, 0, 0); + fmt::print("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; + fmt::print("X+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error); + // Reset position + i.translateCoordinates(-delta, 0, 0); + fmt::print("X coord = {}\n", i.x()); + break; + // Move y and repeat + case 1: + fmt::print("Y coord = {}\n", i.y()); + i.translateCoordinates(0, -delta, 0); + fmt::print("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; + fmt::print("Y+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error); + i.translateCoordinates(0, 2 * delta, 0); + fmt::print("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; + fmt::print("Y+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error); + i.translateCoordinates(0, -delta, 0); + fmt::print("Y coord = {}\n", i.y()); + break; + // Move z + case 2: + fmt::print("Z coord = {}\n", i.z()); + i.translateCoordinates(0, 0, -delta); + fmt::print("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; + fmt::print("Z+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error); + i.translateCoordinates(0, 0, 2 * delta); + fmt::print("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; + fmt::print("Z+ = {}\n", Error::rmse(originalgr.total(), originalReferenceData).error); + i.translateCoordinates(0, 0, -delta); + fmt::print("Z coord = {}\n", i.z()); + break; + } + } + } + } + + for (auto &&[i, f] : zip(atoms, forces)) + { + fmt::print("X coord = {}, Y coord = {}, Z coord = {}\n", i.x(), i.y(), i.z()); + i.translateCoordinates(f); + f.print(); + printf("We've moved\n"); + fmt::print("X coord = {}, Y coord = {}, Z coord = {}\n", i.x(), i.y(), i.z()); + } + + targetConfiguration_->incrementContentsVersion(); + targetConfiguration_->updateObjectRelationships(); + + return ExecutionResult::Success; +} diff --git a/src/modules/registry.cpp b/src/modules/registry.cpp index 4739903a38..9916cc74a5 100644 --- a/src/modules/registry.cpp +++ b/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" @@ -53,6 +54,7 @@ ModuleRegistry::ModuleRegistry() registerProducer(ModuleTypes::Bragg, "Calculate Bragg intensities", "Correlation Functions"); registerProducer(ModuleTypes::Compare, "Compare data sets and calculate errors", "Checks & Tests"); registerProducer(ModuleTypes::DAngle, "Calculate distance/angle maps", "Analysis"); + registerProducer(ModuleTypes::DrivenMD, "SAXS/SANS driven molecular dynamics test", "Forcefield"); registerProducer(ModuleTypes::Energy, "Calculate the total energy of a Configuration", "Forcefield"); registerProducer(ModuleTypes::EPSR, "Refine interatomic potentials in a manner consistent with EPSR", "Forcefield");