From 2343dde3d152150cbae0d6ebb5bfe7f8abbe1135 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Mon, 13 Nov 2023 15:37:54 +0000 Subject: [PATCH 01/21] module added to GUI --- CMakeLists.txt | 38 +-- flake.nix | 4 +- src/module/module.cpp | 3 +- src/module/module.h | 3 +- src/modules/dlPoly/CMakeLists.txt | 1 + src/modules/dlPoly/dlPoly.cpp | 65 +++++ src/modules/dlPoly/dlPoly.h | 87 +++++++ src/modules/dlPoly/functions.cpp | 183 ++++++++++++++ src/modules/dlPoly/process.cpp | 402 ++++++++++++++++++++++++++++++ src/modules/registry.cpp | 2 + 10 files changed, 766 insertions(+), 22 deletions(-) create mode 100644 src/modules/dlPoly/CMakeLists.txt create mode 100644 src/modules/dlPoly/dlPoly.cpp create mode 100644 src/modules/dlPoly/dlPoly.h create mode 100644 src/modules/dlPoly/functions.cpp create mode 100644 src/modules/dlPoly/process.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 0084bc6092..e060626900 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -270,26 +270,26 @@ else() endif() # Enable build of test suite? -option(BUILD_TESTS "Build test suite" OFF) -if(BUILD_TESTS) - enable_testing() - add_subdirectory(examples) - add_subdirectory(tests) -endif(BUILD_TESTS) - -if(BUILD_BENCHMARKS) +#option(BUILD_TESTS "Build test suite" OFF) +#if(BUILD_TESTS) +# enable_testing() +# add_subdirectory(examples) +# add_subdirectory(tests) +#endif(BUILD_TESTS) + +#if(BUILD_BENCHMARKS) # The conan package of goooglebenchmark segfaults - so lets use fetch content - include(FetchContent) - set(BENCHMARK_ENABLE_TESTING "OFF") - fetchcontent_declare( - googlebenchmark - GIT_REPOSITORY https://github.com/google/benchmark.git - GIT_TAG v1.5.4 - ) - - fetchcontent_makeavailable(googlebenchmark) - add_subdirectory(benchmark) -endif(BUILD_BENCHMARKS) +# include(FetchContent) +# set(BENCHMARK_ENABLE_TESTING "OFF") +# fetchcontent_declare( +# googlebenchmark +# GIT_REPOSITORY https://github.com/google/benchmark.git +# GIT_TAG v1.5.4 +# ) +# +# fetchcontent_makeavailable(googlebenchmark) +# add_subdirectory(benchmark) +#endif(BUILD_BENCHMARKS) # Build CLI target executable target_link_libraries( diff --git a/flake.nix b/flake.nix index 7849687464..a5e542e29a 100644 --- a/flake.nix +++ b/flake.nix @@ -22,6 +22,7 @@ version = "1.4.0"; base_libs = pkgs: with pkgs; [ + gfortran antlr4 antlr4.runtime.cpp antlr4.runtime.cpp.dev @@ -176,6 +177,7 @@ export __EGL_VENDOR_LIBRARY_FILENAMES=${pkgs.mesa.drivers}/share/glvnd/egl_vendor.d/50_mesa.json export LD_LIBRARY_PATH=${pkgs.lib.makeLibraryPath [pkgs.mesa.drivers]}:${pkgs.lib.makeSearchPathOutput "lib" "lib/vdpau" [pkgs.libvdpau]}:${pkgs.lib.makeLibraryPath [pkgs.libglvnd]}"''${LD_LIBRARY_PATH:+:$LD_LIBRARY_PATH}" export QT_PLUGIN_PATH="${qt-idaaas.packages.${system}.qtsvg}/lib/qt-6/plugins:$QT_PLUGIN_PATH" + export FC=${pkgs.lib.makeSearchPathOutput "gfortran" "bin/gfortran" [pkgs.gfortran]} ''; @@ -278,4 +280,4 @@ user-env = import ./nix/user-env.nix; }; }); -} +} \ No newline at end of file diff --git a/src/module/module.cpp b/src/module/module.cpp index 66118ca1e6..cc5270a02d 100644 --- a/src/module/module.cpp +++ b/src/module/module.cpp @@ -49,7 +49,8 @@ EnumOptions moduleTypes_("ModuleType", {{ModuleTypes::A {ModuleTypes::SiteRDF, "SiteRDF"}, {ModuleTypes::SQ, "SQ"}, {ModuleTypes::Test, "Test"}, - {ModuleTypes::XRaySQ, "XRaySQ"}}); + {ModuleTypes::XRaySQ, "XRaySQ"}, + {ModuleTypes::DlPoly, "DlPoly"}}); // Return module type string for specified type enumeration std::string moduleType(ModuleTypes::ModuleType type) { return moduleTypes_.keyword(type); } diff --git a/src/module/module.h b/src/module/module.h index 14f0a095bc..7d52bfc6e7 100644 --- a/src/module/module.h +++ b/src/module/module.h @@ -57,7 +57,8 @@ enum ModuleType Skeleton, SQ, Test, - XRaySQ + XRaySQ, + DlPoly }; // Return module type string for specified type enumeration std::string moduleType(ModuleTypes::ModuleType type); diff --git a/src/modules/dlPoly/CMakeLists.txt b/src/modules/dlPoly/CMakeLists.txt new file mode 100644 index 0000000000..4959c98ca8 --- /dev/null +++ b/src/modules/dlPoly/CMakeLists.txt @@ -0,0 +1 @@ +dissolve_add_module(dlPoly.h dlPoly) diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp new file mode 100644 index 0000000000..9e0e86d350 --- /dev/null +++ b/src/modules/dlPoly/dlPoly.cpp @@ -0,0 +1,65 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#include "modules/dlPoly/dlPoly.h" +#include "keywords/bool.h" +#include "keywords/configuration.h" +#include "keywords/double.h" +#include "keywords/integer.h" +#include "keywords/optionalDouble.h" +#include "keywords/optionalInt.h" +#include "keywords/speciesVector.h" + +DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) +{ + keywords_.addTarget("Configuration", "Set target configuration for the module", targetConfiguration_); + + keywords_.setOrganisation("Options", "Simulation"); + keywords_.add("NSteps", "Number of DlPoly steps to perform", nSteps_, 1); + keywords_.add>("Timestep", "Timestep type to use in calculation", timestepType_, + DlPolyModule::timestepType()); + keywords_.add("DeltaT", "Fixed timestep (ps) to use in DlPoly simulation", fixedTimestep_, 0.0); + keywords_.add("RandomVelocities", + "Whether random velocities should always be assigned before beginning DlPoly simulation", + randomVelocities_); + + keywords_.setOrganisation("Options", "Control"); + keywords_ + .add("RestrictToSpecies", "Restrict the calculation to the specified Species", restrictToSpecies_) + ->setEditSignals({KeywordBase::KeywordSignal::ClearModuleData}); + keywords_.add("OnlyWhenEnergyStable", "Only run DlPoly when target Configuration energies are stable", + onlyWhenEnergyStable_); + + keywords_.setOrganisation("Options", "Output"); + keywords_.add("EnergyFrequency", "Frequency at which to calculate total system energy", + energyFrequency_, 0, std::nullopt, 5, "Off"); + keywords_.add("OutputFrequency", "Frequency at which to output step information", outputFrequency_, + 0, std::nullopt, 5, "Off"); + keywords_.add("TrajectoryFrequency", "Write frequency for trajectory file", trajectoryFrequency_, 0, + std::nullopt, 5, "Off"); + + keywords_.setOrganisation("Advanced"); + keywords_.add("CapForces", "Control whether atomic forces are capped every step", capForces_); + keywords_.add("CapForcesAt", "Set cap on allowable force (kJ/mol) per atom", capForcesAt_, 0.0); + keywords_.add( + "CutoffDistance", "Interatomic cutoff distance to use for energy calculation (0.0 to use pair potential range)", + cutoffDistance_, 0.0, std::nullopt, 0.1, "Use PairPotential Range"); + keywords_.add( + "IntraOnly", + "Only forces arising from intramolecular terms (including pair potential contributions) will be calculated", + intramolecularForcesOnly_); + + // Deprecated + static bool deprecatedBool_{false}; + keywords_.addDeprecated("VariableTimestep", + "Whether a variable timestep should be used, determined from the maximal force vector", + deprecatedBool_); +} + +// Return enum options for TimestepType +EnumOptions DlPolyModule::timestepType() +{ + return EnumOptions( + "TimestepType", + {{TimestepType::Fixed, "Fixed"}, {TimestepType::Variable, "Variable"}, {TimestepType::Automatic, "Auto"}}); +} diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h new file mode 100644 index 0000000000..952a86f088 --- /dev/null +++ b/src/modules/dlPoly/dlPoly.h @@ -0,0 +1,87 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#pragma once + +#include "base/enumOptions.h" +#include "module/module.h" + +// Forward Declarations +class PotentialMap; +class Species; + +// Molecular Dynamics Module +class DlPolyModule : public Module +{ + public: + DlPolyModule(); + ~DlPolyModule() override = default; + + /* + * Definition + */ + public: + // Timestep Type + enum class TimestepType + { + Fixed, + Variable, + Automatic + }; + // Return enum options for TimestepType + static EnumOptions timestepType(); + + private: + // Target configurations + Configuration *targetConfiguration_{nullptr}; + // Control whether atomic forces are capped every step + bool capForces_{false}; + // Set cap on allowable force (kJ/mol) per atom + double capForcesAt_{1.0e7}; + // Interatomic cutoff distance to employ + std::optional cutoffDistance_; + // Timestep type to employ + TimestepType timestepType_{TimestepType::Automatic}; + // Fixed timestep (ps) to use in MD simulation + double fixedTimestep_{5.0e-4}; + // Frequency at which to calculate total system energy + std::optional energyFrequency_{10}; + // Whether to restrict force calculation to intramolecular contributions only + bool intramolecularForcesOnly_{false}; + // Number of steps to perform + int nSteps_{50}; + // Only run MD when target Configuration energies are stable + bool onlyWhenEnergyStable_{true}; + // Frequency at which to output step information + std::optional outputFrequency_{5}; + // Whether random velocities should always be assigned before beginning MD simulation + bool randomVelocities_{false}; + // Species to restrict calculation to + std::vector restrictToSpecies_; + // Write frequency for trajectory file + std::optional trajectoryFrequency_; + + /* + * Functions + */ + private: + // Cap forces in Configuration + static int capForces(double maxForceSq, std::vector> &fInter, std::vector> &fIntra); + // Determine timestep to use + static std::optional determineTimeStep(TimestepType timestepType, double requestedTimeStep, + const std::vector> &fInter, + const std::vector> &fIntra); + + public: + // Evolve Species coordinates, returning new coordinates + static std::vector> evolve(const ProcessPool &procPool, const PotentialMap &potentialMap, const Species *sp, + double temperature, int nSteps, double maxDeltaT, + const std::vector> &rInit, std::vector> &velocities); + + /* + * Processing + */ + private: + // Run main processing + Module::ExecutionResult process(ModuleContext &moduleContext) override; +}; diff --git a/src/modules/dlPoly/functions.cpp b/src/modules/dlPoly/functions.cpp new file mode 100644 index 0000000000..93e9e1b262 --- /dev/null +++ b/src/modules/dlPoly/functions.cpp @@ -0,0 +1,183 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#include "classes/configuration.h" +#include "classes/species.h" +#include "data/atomicMasses.h" +#include "modules/forces/forces.h" +#include "modules/dlPoly/dlPoly.h" + +// Cap forces in Configuration +int DlPolyModule::capForces(double maxForce, std::vector> &fInter, std::vector> &fIntra) +{ + double fMag; + const auto maxForceSq = maxForce * maxForce; + auto nCapped = 0; + for (auto &&[inter, intra] : zip(fInter, fIntra)) + { + fMag = (inter + intra).magnitudeSq(); + if (fMag < maxForceSq) + continue; + + fMag = maxForce / sqrt(fMag); + inter *= fMag; + intra *= fMag; + + ++nCapped; + } + + return nCapped; +} + +// Determine timestep to use +std::optional DlPolyModule::determineTimeStep(TimestepType timestepType, double requestedTimeStep, + const std::vector> &fInter, + const std::vector> &fIntra) +{ + if (timestepType == TimestepType::Fixed) + return requestedTimeStep; + + // Simple variable timestep + if (timestepType == TimestepType::Variable) + { + auto absFMax = 0.0; + for (auto &&[inter, intra] : zip(fInter, fIntra)) + absFMax = std::max(absFMax, (inter + intra).absMax()); + + return 1.0 / absFMax; + } + + // Automatic timestep determination, using maximal interatomic force to guide the timestep up to the current fixed timestep + // value + auto absFMaxInter = + std::max_element(fInter.begin(), fInter.end(), [](auto &left, auto &right) { return left.absMax() < right.absMax(); }) + ->absMax(); + + auto deltaT = 100.0 / absFMaxInter; + if (deltaT < (requestedTimeStep / 100.0)) + return {}; + return deltaT > requestedTimeStep ? requestedTimeStep : deltaT; +} + +// Evolve Species coordinates, returning new coordinates +std::vector> DlPolyModule::evolve(const ProcessPool &procPool, const PotentialMap &potentialMap, const Species *sp, + double temperature, int nSteps, double maxDeltaT, + const std::vector> &rInit, std::vector> &velocities) +{ + assert(sp); + assert(sp->nAtoms() == velocities.size()); + + // Create arrays + std::vector mass(sp->nAtoms(), 0.0); + std::vector> fInter(sp->nAtoms()), fIntra(sp->nAtoms()), accelerations(sp->nAtoms()); + + // Variables + auto &atoms = sp->atoms(); + double tInstant, ke, tScale; + + // Units + // J = kg m2 s-2 --> 10 J = g Ang2 ps-2 + // If ke is in units of [g mol-1 Angstroms2 ps-2] then must use kb in units of 10 J mol-1 K-1 (= 0.8314462) + const auto kb = 0.8314462; + + // Store atomic masses for future use + std::transform(atoms.begin(), atoms.end(), mass.begin(), [](const auto &atom) { return AtomicMass::mass(atom.Z()); }); + + // Calculate total velocity and mass over all atoms + auto massSum = std::accumulate(mass.begin(), mass.end(), 0.0); + auto vCom = std::transform_reduce(velocities.begin(), velocities.end(), mass.begin(), Vec3()); + + // Remove any velocity shift + vCom /= massSum; + std::transform(velocities.begin(), velocities.end(), velocities.begin(), [vCom](auto vel) { return vel - vCom; }); + + // Calculate instantaneous temperature + ke = std::transform_reduce(mass.begin(), mass.end(), velocities.begin(), 0.0, std::plus<>(), + [](const auto m, const auto &v) { return 0.5 * m * v.dp(v); }); + tInstant = ke * 2.0 / (3.0 * atoms.size() * kb); + + // Rescale velocities for desired temperature + tScale = sqrt(temperature / tInstant); + std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto v) { return v * tScale; }); + + // Zero force arrays + std::fill(fInter.begin(), fInter.end(), Vec3()); + std::fill(fIntra.begin(), fIntra.end(), Vec3()); + + ForcesModule::totalForces(procPool, sp, potentialMap, ForcesModule::ForceCalculationType::Full, fInter, fIntra, rInit); + + // Must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) + std::transform(fInter.begin(), fInter.end(), fInter.begin(), [](auto f) { return f * 100.0; }); + std::transform(fIntra.begin(), fIntra.end(), fIntra.begin(), [](auto f) { return f * 100.0; }); + + // Check for suitable timestep + if (!determineTimeStep(TimestepType::Automatic, maxDeltaT, fInter, fIntra)) + { + Messenger::print("Forces are currently too high for species MD to proceed. Try decreasing the maximum timestep.\n"); + return rInit; + } + + // Copy coordinates ready for propagation + auto rNew = rInit; + + // Ready to do MD propagation of the species + for (auto step = 1; step <= nSteps; ++step) + { + // Get timestep + auto optDT = determineTimeStep(TimestepType::Automatic, maxDeltaT, fInter, fIntra); + if (!optDT) + { + Messenger::warn("A reasonable timestep could not be determined. Stopping evolution.\n"); + break; + } + auto dT = *optDT; + auto deltaTSq = dT * dT; + + // Velocity Verlet first stage (A) + // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 + // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt + // B: a(t+dt) = F(t+dt)/m + // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt + for (auto &&[r, v, a] : zip(rNew, velocities, accelerations)) + { + // Propagate positions (by whole step)... + r += v * dT + a * 0.5 * deltaTSq; + + // ...velocities (by half step)... + v += a * 0.5 * dT; + } + + // Zero force arrays + std::fill(fInter.begin(), fInter.end(), Vec3()); + std::fill(fIntra.begin(), fIntra.end(), Vec3()); + + // Calculate forces - must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) + ForcesModule::totalForces(procPool, sp, potentialMap, ForcesModule::ForceCalculationType::Full, fInter, fIntra, rNew); + std::transform(fInter.begin(), fInter.end(), fInter.begin(), [](auto f) { return f * 100.0; }); + std::transform(fIntra.begin(), fIntra.end(), fIntra.begin(), [](auto f) { return f * 100.0; }); + + // Velocity Verlet second stage (B) and velocity scaling + // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 + // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt + // B: a(t+dt) = F(t+dt)/m + // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt + ke = 0.0; + for (auto &&[f1, f2, v, a, m] : zip(fInter, fIntra, velocities, accelerations, mass)) + { + // Determine new accelerations + a = (f1 + f2) / m; + + // ..and finally velocities again (by second half-step) + v += a * 0.5 * dT; + + ke += 0.5 * m * v.dp(v); + } + + // Rescale velocities for desired temperature + tInstant = ke * 2.0 / (3.0 * sp->nAtoms() * kb); + tScale = sqrt(temperature / tInstant); + std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto &v) { return v * tScale; }); + } + + return rNew; +} diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp new file mode 100644 index 0000000000..2525a89629 --- /dev/null +++ b/src/modules/dlPoly/process.cpp @@ -0,0 +1,402 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#include "base/lineParser.h" +#include "base/randomBuffer.h" +#include "base/timer.h" +#include "data/atomicMasses.h" +#include "main/dissolve.h" +#include "module/context.h" +#include "modules/energy/energy.h" +#include "modules/forces/forces.h" +#include "modules/dlPoly/dlPoly.h" + +// Run main processing +Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) +{ + // Check for zero Configuration targets + if (!targetConfiguration_) + { + Messenger::error("No configuration target set for module '{}'.\n", name()); + return ExecutionResult::Failed; + } + + // Get control parameters + const auto maxForce = capForcesAt_ * 100.0; // To convert from kJ/mol to 10 J/mol + auto rCut = cutoffDistance_.value_or(moduleContext.dissolve().pairPotentialRange()); + + // Units + // J = kg m2 s-2 --> 10 J = g Ang2 ps-2 + // If ke is in units of [g mol-1 Angstroms2 ps-2] then must use kb in units of 10 J mol-1 K-1 (= 0.8314462) + const auto kb = 0.8314462; + + // Print argument/parameter summary + Messenger::print("MD: Cutoff distance is {}\n", rCut); + Messenger::print("MD: Number of steps = {}\n", nSteps_); + Messenger::print("MD: Timestep type is '{}'\n", timestepType().keyword(timestepType_)); + if (onlyWhenEnergyStable_) + Messenger::print("MD: Only perform MD if target Configuration energies are stable.\n"); + if (trajectoryFrequency_.value_or(0) > 0) + Messenger::print("MD: Trajectory file will be appended every {} step(s).\n", trajectoryFrequency_.value()); + else + Messenger::print("MD: Trajectory file off.\n"); + if (capForces_) + Messenger::print("MD: Forces will be capped to {:10.3e} kJ/mol per atom per axis.\n", maxForce / 100.0); + if (energyFrequency_.value_or(0) > 0) + Messenger::print("MD: Energy will be calculated every {} step(s).\n", energyFrequency_.value()); + else + Messenger::print("MD: Energy will be not be calculated.\n"); + if (outputFrequency_.value_or(0) > 0) + Messenger::print("MD: Summary will be written every {} step(s).\n", outputFrequency_.value()); + else + Messenger::print("MD: Summary will not be written.\n"); + if (!restrictToSpecies_.empty()) + Messenger::print("MD: Calculation will be restricted to species: {}\n", + joinStrings(restrictToSpecies_, " ", [](const auto &sp) { return sp->name(); })); + Messenger::print("\n"); + + if (onlyWhenEnergyStable_) + { + auto stabilityResult = + EnergyModule::checkStability(moduleContext.dissolve().processingModuleData(), targetConfiguration_); + if (stabilityResult == EnergyModule::NotAssessable) + { + return ExecutionResult::Failed; + } + else if (stabilityResult == EnergyModule::EnergyUnstable) + { + Messenger::print("Skipping MD for Configuration '{}'.\n", targetConfiguration_->niceName()); + return ExecutionResult::NotExecuted; + } + } + + // Get temperature from Configuration + const auto temperature = targetConfiguration_->temperature(); + + // Create arrays + std::vector mass(targetConfiguration_->nAtoms(), 0.0); + std::vector> fBound(targetConfiguration_->nAtoms()), fUnbound(targetConfiguration_->nAtoms()), + accelerations(targetConfiguration_->nAtoms()); + + // Variables + auto nCapped = 0; + auto &atoms = targetConfiguration_->atoms(); + double tInstant, ke, tScale, peInter, peIntra; + + // Determine target molecules from the restrictedSpecies vector (if any) + std::vector targetMolecules; + std::vector free(targetConfiguration_->nAtoms(), 0); + if (restrictToSpecies_.empty()) + { + std::fill(free.begin(), free.end(), 1); + } + else + for (const auto &mol : targetConfiguration_->molecules()) + if (std::find(restrictToSpecies_.begin(), restrictToSpecies_.end(), mol->species()) != restrictToSpecies_.end()) + { + targetMolecules.push_back(mol.get()); + auto offset = mol->globalAtomOffset(); + std::fill(free.begin() + offset, free.begin() + offset + mol->atoms().size(), 1); + } + + /* + * Calculation Begins + */ + + // Initialise the random number buffer for all processes + RandomBuffer randomBuffer(moduleContext.processPool(), ProcessPool::PoolProcessesCommunicator); + + // Read in or assign random velocities + auto [velocities, status] = moduleContext.dissolve().processingModuleData().realiseIf>>( + fmt::format("{}//Velocities", targetConfiguration_->niceName()), name(), GenericItem::InRestartFileFlag); + if ((status == GenericItem::ItemStatus::Created || randomVelocities_ || + velocities.size() != targetConfiguration_->nAtoms()) && + !intramolecularForcesOnly_) + { + // Show warning message on array size mismatch + if (status != GenericItem::ItemStatus::Created && velocities.size() != targetConfiguration_->nAtoms()) + Messenger::warn( + "Size of existing velocities array doesn't match the current configuration size - they will be ignored."); + + Messenger::print("Random initial velocities will be assigned.\n"); + velocities.resize(targetConfiguration_->nAtoms(), Vec3()); + for (auto &&[v, iFree] : zip(velocities, free)) + { + if (iFree) + v.set(exp(randomBuffer.random() - 0.5), exp(randomBuffer.random() - 0.5), exp(randomBuffer.random() - 0.5)); + else + v.zero(); + v /= sqrt(TWOPI); + } + } + else if (intramolecularForcesOnly_) + { + Messenger::print("Only intramolecular forces will be calculated, so velocities will be zeroes.\n"); + velocities.resize(targetConfiguration_->nAtoms(), Vec3()); + std::fill(velocities.begin(), velocities.end(), Vec3()); + } + else + { + Messenger::print("Existing velocities will be used.\n"); + } + + Messenger::print("\n"); + + // Store atomic masses for future use + for (auto &&[i, m] : zip(atoms, mass)) + m = AtomicMass::mass(i.speciesAtom()->Z()); + + // Calculate total velocity and mass over all atoms + Vec3 vCom; + auto massSum = 0.0; + for (auto &&[v, m, iFree] : zip(velocities, mass, free)) + { + if (!iFree) + continue; + vCom += v * m; + massSum += m; + } + + // Finalise initial velocities (unless considering intramolecular forces only) + if (!intramolecularForcesOnly_) + { + // Remove any velocity shift, and re-zero velocities on fixed atoms + vCom /= massSum; + std::transform(velocities.begin(), velocities.end(), velocities.begin(), [vCom](auto vel) { return vel - vCom; }); + for (auto &&[v, iFree] : zip(velocities, free)) + if (!iFree) + v.zero(); + + // Calculate instantaneous temperature + ke = 0.0; + for (auto &&[m, v] : zip(mass, velocities)) + ke += 0.5 * m * v.dp(v); + tInstant = ke * 2.0 / (3.0 * atoms.size() * kb); + + // Rescale velocities for desired temperature + tScale = sqrt(temperature / tInstant); + std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto v) { return v * tScale; }); + } + + // Open trajectory file (if requested) + LineParser trajParser; + if (trajectoryFrequency_.value_or(0) > 0) + { + std::string trajectoryFile = fmt::format("{}.md.xyz", targetConfiguration_->name()); + if (moduleContext.processPool().isMaster()) + { + if ((!trajParser.appendOutput(trajectoryFile)) || (!trajParser.isFileGoodForWriting())) + { + Messenger::error("Failed to open MD trajectory output file '{}'.\n", trajectoryFile); + moduleContext.processPool().decideFalse(); + return ExecutionResult::Failed; + } + moduleContext.processPool().decideTrue(); + } + else if (!moduleContext.processPool().decision()) + { + return ExecutionResult::Failed; + } + } + + // Write header + if (outputFrequency_.value_or(0) > 0) + { + Messenger::print(" Energies (kJ/mol)\n"); + Messenger::print(" Step T(K) Kinetic Inter Intra Total " + "deltaT(ps)\n"); + } + + // Start a timer + Timer timer, commsTimer(false); + + // If we're not using a fixed timestep the forces need to be available immediately + if (timestepType_ != TimestepType::Fixed) + { + // Zero force arrays + std::fill(fUnbound.begin(), fUnbound.end(), Vec3()); + std::fill(fBound.begin(), fBound.end(), Vec3()); + + if (targetMolecules.empty()) + ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, + moduleContext.dissolve().potentialMap(), + intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull + : ForcesModule::ForceCalculationType::Full, + fUnbound, fBound, commsTimer); + else + ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, targetMolecules, + moduleContext.dissolve().potentialMap(), + intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull + : ForcesModule::ForceCalculationType::Full, + fUnbound, fBound, commsTimer); + + // Must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) + std::transform(fUnbound.begin(), fUnbound.end(), fUnbound.begin(), [](auto f) { return f * 100.0; }); + std::transform(fBound.begin(), fBound.end(), fBound.begin(), [](auto f) { return f * 100.0; }); + + // Check for suitable timestep + if (!determineTimeStep(timestepType_, fixedTimestep_, fUnbound, fBound)) + { + Messenger::print("Forces are currently too high for MD to proceed. Skipping this run.\n"); + return ExecutionResult::NotExecuted; + } + } + + // Ready to do MD propagation of system + auto step = 1; + for (step = 1; step <= nSteps_; ++step) + { + // Get timestep + auto optDT = determineTimeStep(timestepType_, fixedTimestep_, fUnbound, fBound); + if (!optDT) + { + Messenger::warn("A reasonable timestep could not be determined. Stopping evolution.\n"); + break; + } + auto dT = *optDT; + auto deltaTSq = dT * dT; + + // Velocity Verlet first stage (A) + // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 + // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt + // B: a(t+dt) = F(t+dt)/m + // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt + for (auto &&[i, v, a] : zip(atoms, velocities, accelerations)) + { + // Propagate positions (by whole step)... + i.translateCoordinates(v * dT + a * 0.5 * deltaTSq); + + // ...velocities (by half step)... + v += a * 0.5 * dT; + } + + // Update Atom locations + targetConfiguration_->updateAtomLocations(); + + // Zero force arrays + std::fill(fUnbound.begin(), fUnbound.end(), Vec3()); + std::fill(fBound.begin(), fBound.end(), Vec3()); + + // Calculate forces - must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) + if (targetMolecules.empty()) + ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, + moduleContext.dissolve().potentialMap(), + intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull + : ForcesModule::ForceCalculationType::Full, + fUnbound, fBound, commsTimer); + else + ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, targetMolecules, + moduleContext.dissolve().potentialMap(), + intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull + : ForcesModule::ForceCalculationType::Full, + fUnbound, fBound, commsTimer); + std::transform(fUnbound.begin(), fUnbound.end(), fUnbound.begin(), [](auto f) { return f * 100.0; }); + std::transform(fBound.begin(), fBound.end(), fBound.begin(), [](auto f) { return f * 100.0; }); + + // Cap forces + if (capForces_) + nCapped = capForces(maxForce, fUnbound, fBound); + + // Velocity Verlet second stage (B) and velocity scaling + // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 + // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt + // B: a(t+dt) = F(t+dt)/m + // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt + ke = 0.0; + for (auto &&[f1, f2, v, a, m] : zip(fUnbound, fBound, velocities, accelerations, mass)) + { + // Determine new accelerations + a = (f1 + f2) / m; + + // ..and finally velocities again (by second half-step) + v += a * 0.5 * dT; + + ke += 0.5 * m * v.dp(v); + } + + // Rescale velocities for desired temperature + tInstant = ke * 2.0 / (3.0 * targetConfiguration_->nAtoms() * kb); + tScale = sqrt(temperature / tInstant); + std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto &v) { return v * tScale; }); + + // Convert ke from 10J/mol to kJ/mol + ke *= 0.01; + + // Write step summary? + if (outputFrequency_ && (step == 1 || (step % outputFrequency_.value() == 0))) + { + // Include total energy term? + if (energyFrequency_ && (step % energyFrequency_.value() == 0)) + { + peInter = EnergyModule::interAtomicEnergy(moduleContext.processPool(), targetConfiguration_, + moduleContext.dissolve().potentialMap()); + peIntra = EnergyModule::intraMolecularEnergy(moduleContext.processPool(), targetConfiguration_, + moduleContext.dissolve().potentialMap()); + Messenger::print(" {:<10d} {:10.3e} {:10.3e} {:10.3e} {:10.3e} {:10.3e} {:10.3e}\n", step, + tInstant, ke, peInter, peIntra, ke + peIntra + peInter, dT); + } + else + Messenger::print(" {:<10d} {:10.3e} {:10.3e} {:10.3e}\n", step, + tInstant, ke, dT); + } + + // Save trajectory frame + if (trajectoryFrequency_ && (step % trajectoryFrequency_.value() == 0)) + { + if (moduleContext.processPool().isMaster()) + { + // Write number of atoms + trajParser.writeLineF("{}\n", targetConfiguration_->nAtoms()); + + // Construct and write header + std::string header = fmt::format("Step {} of {}, T = {:10.3e}, ke = {:10.3e}", step, nSteps_, tInstant, ke); + if (energyFrequency_ && (step % energyFrequency_.value() == 0)) + header += fmt::format(", inter = {:10.3e}, intra = {:10.3e}, tot = {:10.3e}", peInter, peIntra, + ke + peInter + peIntra); + if (!trajParser.writeLine(header)) + { + moduleContext.processPool().decideFalse(); + return ExecutionResult::Failed; + } + + // Write Atoms + for (const auto &i : atoms) + { + if (!trajParser.writeLineF("{:<3} {:10.3f} {:10.3f} {:10.3f}\n", Elements::symbol(i.speciesAtom()->Z()), + i.r().x, i.r().y, i.r().z)) + { + moduleContext.processPool().decideFalse(); + return ExecutionResult::Failed; + } + } + + moduleContext.processPool().decideTrue(); + } + else if (!moduleContext.processPool().decision()) + { + return ExecutionResult::Failed; + } + } + } + timer.stop(); + + // Close trajectory file + if (trajectoryFrequency_.value_or(0) > 0 && moduleContext.processPool().isMaster()) + trajParser.closeFiles(); + + if (capForces_) + Messenger::print("A total of {} forces were capped over the course of the dynamics ({:9.3e} per step).\n", nCapped, + double(nCapped) / nSteps_); + Messenger::print("{} steps performed ({} work, {} comms)\n", step - 1, timer.totalTimeString(), + commsTimer.totalTimeString()); + + // Increment configuration changeCount + if (step > 1) + targetConfiguration_->incrementContentsVersion(); + + /* + * Calculation End + */ + + return ExecutionResult::Success; +} diff --git a/src/modules/registry.cpp b/src/modules/registry.cpp index a4ca7925a0..7f64c1ae17 100644 --- a/src/modules/registry.cpp +++ b/src/modules/registry.cpp @@ -35,6 +35,7 @@ #include "modules/sq/sq.h" #include "modules/test/test.h" #include "modules/xRaySQ/xRaySQ.h" +#include "modules/dlPoly/dlPoly.h" ModuleRegistry::ModuleRegistry() { @@ -85,6 +86,7 @@ ModuleRegistry::ModuleRegistry() registerProducer(ModuleTypes::SQ, "Transform g(r) into unweighted S(Q)", "Correlation Functions"); registerProducer(ModuleTypes::Test, "Development Module"); registerProducer(ModuleTypes::XRaySQ, "Calculate x-ray-weighted S(Q)", "Correlation Functions"); + registerProducer(ModuleTypes::DlPoly, "Dl poly", "Evolution"); } /* From 52a7a3e687f8a919f5a83d7256acdc56f8078a93 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 15 Nov 2023 17:10:04 +0000 Subject: [PATCH 02/21] export CONTROL file classes implemented --- src/io/export/CMakeLists.txt | 2 ++ src/modules/dlPoly/dlPoly.cpp | 15 ++++----------- src/modules/dlPoly/dlPoly.h | 3 +++ src/modules/dlPoly/process.cpp | 18 ++++++++++++++++++ 4 files changed, 27 insertions(+), 11 deletions(-) diff --git a/src/io/export/CMakeLists.txt b/src/io/export/CMakeLists.txt index 2bed6dfd51..36c40c1311 100644 --- a/src/io/export/CMakeLists.txt +++ b/src/io/export/CMakeLists.txt @@ -9,6 +9,7 @@ add_library( species.cpp trajectory.cpp values.cpp + dlPolyControl.cpp coordinates.h data1D.h data2D.h @@ -18,6 +19,7 @@ add_library( species.h trajectory.h values.h + dlPolyControl.h ) target_include_directories(export PRIVATE ${PROJECT_SOURCE_DIR}/src ${PROJECT_BINARY_DIR}/src) diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 9e0e86d350..311421bd5d 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -9,10 +9,14 @@ #include "keywords/optionalDouble.h" #include "keywords/optionalInt.h" #include "keywords/speciesVector.h" +#include "keywords/fileAndFormat.h" DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) { keywords_.addTarget("Configuration", "Set target configuration for the module", targetConfiguration_); + + keywords_.setOrganisation("Options", "File"); + keywords_.add("Format", "File / format for coordinates", dlPolyControlFormat_, "EndFormat"); keywords_.setOrganisation("Options", "Simulation"); keywords_.add("NSteps", "Number of DlPoly steps to perform", nSteps_, 1); @@ -23,13 +27,6 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) "Whether random velocities should always be assigned before beginning DlPoly simulation", randomVelocities_); - keywords_.setOrganisation("Options", "Control"); - keywords_ - .add("RestrictToSpecies", "Restrict the calculation to the specified Species", restrictToSpecies_) - ->setEditSignals({KeywordBase::KeywordSignal::ClearModuleData}); - keywords_.add("OnlyWhenEnergyStable", "Only run DlPoly when target Configuration energies are stable", - onlyWhenEnergyStable_); - keywords_.setOrganisation("Options", "Output"); keywords_.add("EnergyFrequency", "Frequency at which to calculate total system energy", energyFrequency_, 0, std::nullopt, 5, "Off"); @@ -44,10 +41,6 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.add( "CutoffDistance", "Interatomic cutoff distance to use for energy calculation (0.0 to use pair potential range)", cutoffDistance_, 0.0, std::nullopt, 0.1, "Use PairPotential Range"); - keywords_.add( - "IntraOnly", - "Only forces arising from intramolecular terms (including pair potential contributions) will be calculated", - intramolecularForcesOnly_); // Deprecated static bool deprecatedBool_{false}; diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 952a86f088..6f5eecf92d 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -3,6 +3,7 @@ #pragma once +#include "io/export/dlPolyControl.h" #include "base/enumOptions.h" #include "module/module.h" @@ -60,6 +61,8 @@ class DlPolyModule : public Module std::vector restrictToSpecies_; // Write frequency for trajectory file std::optional trajectoryFrequency_; + // Filename and format for CONTROL export + DlPolyControlExportFileFormat dlPolyControlFormat_; /* * Functions diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 2525a89629..473b6a2661 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -21,6 +21,24 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } + // Only the pool master saves the data + if (moduleContext.processPool().isMaster()) + { + Messenger::print("Export: Writing TESTTEST file ({}) for Configuration '{}'...\n", + dlPolyControlFormat_.formatDescription(), targetConfiguration_->name()); + + if (!dlPolyControlFormat_.exportData(targetConfiguration_)) + { + Messenger::print("Export: Failed to export coordinates file '{}'.\n", dlPolyControlFormat_.filename()); + moduleContext.processPool().decideFalse(); + return ExecutionResult::Failed; + } + + moduleContext.processPool().decideTrue(); + } + else if (!moduleContext.processPool().decision()) + return ExecutionResult::Failed; + // Get control parameters const auto maxForce = capForcesAt_ * 100.0; // To convert from kJ/mol to 10 J/mol auto rCut = cutoffDistance_.value_or(moduleContext.dissolve().pairPotentialRange()); From f8ac94b173aed723631cdeacfc75d6d5b0a0dbf4 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 15 Nov 2023 17:13:53 +0000 Subject: [PATCH 03/21] export CONTROL file classes implemented --- src/io/export/dlPolyControl.cpp | 56 +++++++++++++++++++++++++++++++++ src/io/export/dlPolyControl.h | 47 +++++++++++++++++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 src/io/export/dlPolyControl.cpp create mode 100644 src/io/export/dlPolyControl.h diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp new file mode 100644 index 0000000000..09ded8222d --- /dev/null +++ b/src/io/export/dlPolyControl.cpp @@ -0,0 +1,56 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#include "io/export/dlPolyControl.h" +#include "base/lineParser.h" +#include "base/sysFunc.h" +#include "classes/atomType.h" +#include "classes/box.h" +#include "classes/configuration.h" +#include "classes/speciesAtom.h" +#include "data/atomicMasses.h" + +DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view filename, DlPolyControlExportFormat format) + : FileAndFormat(formats_, filename, (int)format) +{ + formats_ = EnumOptions( + "CoordinateExportFileFormat", {{DlPolyControlExportFormat::DLPOLY, "dlpoly", "DL_POLY CONFIG File"}}); +} + + +// Export DlPolyControl as CONTROL +bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg) +{ + // Export title + if (!parser.writeLineF("{} @ {}\n", cfg->name(), cfg->contentsVersion())) + return false; + + + return true; +} + +// Export DlPolyControl using current filename and format +bool DlPolyControlExportFileFormat::exportData(Configuration *cfg) +{ + // Open the file + LineParser parser; + if (!parser.openOutput(filename_)) + { + parser.closeFiles(); + return false; + } + + // Write data + auto result = false; + switch (formats_.enumerationByIndex(*formatIndex_)) + { + case (DlPolyControlExportFormat::DLPOLY): + result = exportDLPOLY(parser, cfg); + break; + default: + throw(std::runtime_error(fmt::format("DlPolyControl format '{}' export has not been implemented.\n", + formats_.keywordByIndex(*formatIndex_)))); + } + + return true; +} diff --git a/src/io/export/dlPolyControl.h b/src/io/export/dlPolyControl.h new file mode 100644 index 0000000000..73c6422def --- /dev/null +++ b/src/io/export/dlPolyControl.h @@ -0,0 +1,47 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#pragma once + +#include "io/fileAndFormat.h" + +// Forward Declarations +class Configuration; + +// Coordinate Export Formats +class DlPolyControlExportFileFormat : public FileAndFormat +{ + public: + // Available DlPolyControl formats + enum class DlPolyControlExportFormat + { + DLPOLY + }; + DlPolyControlExportFileFormat(std::string_view filename = "", DlPolyControlExportFormat format = DlPolyControlExportFormat::DLPOLY); + ~DlPolyControlExportFileFormat() override = default; + + /* + * Formats + */ + private: + // Format enum options + EnumOptions formats_; + + /* + * Filename / Basename + */ + public: + // Return whether the file must exist + bool fileMustExist() const override { return false; } + + /* + * Export Functions + */ + private: + // Export Configuration as DL_POLY CONTROL + bool exportDLPOLY(LineParser &parser, Configuration *cfg); + + public: + // Export Configuration using current filename and format + bool exportData(Configuration *cfg); +}; From f14bdbc3e8c925a13d9e8859eda84da21ea039ae Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 15 Nov 2023 20:54:12 +0000 Subject: [PATCH 04/21] add all props to control file exporter --- src/io/export/dlPolyControl.cpp | 18 +- src/io/export/dlPolyControl.h | 7 +- src/modules/dlPoly/dlPoly.h | 14 +- src/modules/dlPoly/process.cpp | 384 +------------------------------- 4 files changed, 36 insertions(+), 387 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 09ded8222d..3a041b4cb0 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -19,18 +19,20 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi // Export DlPolyControl as CONTROL -bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg) +bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) { // Export title if (!parser.writeLineF("{} @ {}\n", cfg->name(), cfg->contentsVersion())) return false; + if (!parser.writeLineF("{}\n", nSteps)) + return false; return true; } // Export DlPolyControl using current filename and format -bool DlPolyControlExportFileFormat::exportData(Configuration *cfg) +bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) { // Open the file LineParser parser; @@ -45,7 +47,17 @@ bool DlPolyControlExportFileFormat::exportData(Configuration *cfg) switch (formats_.enumerationByIndex(*formatIndex_)) { case (DlPolyControlExportFormat::DLPOLY): - result = exportDLPOLY(parser, cfg); + result = exportDLPOLY(parser, + cfg, + capForces, + capForcesAt, + cutoffDistance, + fixedTimestep, + energyFrequency, + nSteps, + outputFrequency, + randomVelocities, + trajectoryFrequency); break; default: throw(std::runtime_error(fmt::format("DlPolyControl format '{}' export has not been implemented.\n", diff --git a/src/io/export/dlPolyControl.h b/src/io/export/dlPolyControl.h index 73c6422def..8b41bda0fb 100644 --- a/src/io/export/dlPolyControl.h +++ b/src/io/export/dlPolyControl.h @@ -39,9 +39,12 @@ class DlPolyControlExportFileFormat : public FileAndFormat */ private: // Export Configuration as DL_POLY CONTROL - bool exportDLPOLY(LineParser &parser, Configuration *cfg); + bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); public: // Export Configuration using current filename and format - bool exportData(Configuration *cfg); + bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); + }; + + diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 6f5eecf92d..4c67127a7e 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -31,6 +31,10 @@ class DlPolyModule : public Module }; // Return enum options for TimestepType static EnumOptions timestepType(); + + // Filename and format for CONTROL export + DlPolyControlExportFileFormat dlPolyControlFormat_; + private: // Target configurations @@ -47,23 +51,15 @@ class DlPolyModule : public Module double fixedTimestep_{5.0e-4}; // Frequency at which to calculate total system energy std::optional energyFrequency_{10}; - // Whether to restrict force calculation to intramolecular contributions only - bool intramolecularForcesOnly_{false}; // Number of steps to perform int nSteps_{50}; - // Only run MD when target Configuration energies are stable - bool onlyWhenEnergyStable_{true}; // Frequency at which to output step information std::optional outputFrequency_{5}; // Whether random velocities should always be assigned before beginning MD simulation bool randomVelocities_{false}; - // Species to restrict calculation to - std::vector restrictToSpecies_; // Write frequency for trajectory file std::optional trajectoryFrequency_; - // Filename and format for CONTROL export - DlPolyControlExportFileFormat dlPolyControlFormat_; - + /* * Functions */ diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 473b6a2661..5d84d05bf9 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -21,13 +21,22 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } - // Only the pool master saves the data + // Save DlPoly CONTROL file with props if (moduleContext.processPool().isMaster()) { Messenger::print("Export: Writing TESTTEST file ({}) for Configuration '{}'...\n", dlPolyControlFormat_.formatDescription(), targetConfiguration_->name()); - if (!dlPolyControlFormat_.exportData(targetConfiguration_)) + if (!dlPolyControlFormat_.exportData(targetConfiguration_, + capForces_, + capForcesAt_, + cutoffDistance_, + fixedTimestep_, + energyFrequency_, + nSteps_, + outputFrequency_, + randomVelocities_, + trajectoryFrequency_)) { Messenger::print("Export: Failed to export coordinates file '{}'.\n", dlPolyControlFormat_.filename()); moduleContext.processPool().decideFalse(); @@ -39,378 +48,7 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) else if (!moduleContext.processPool().decision()) return ExecutionResult::Failed; - // Get control parameters - const auto maxForce = capForcesAt_ * 100.0; // To convert from kJ/mol to 10 J/mol - auto rCut = cutoffDistance_.value_or(moduleContext.dissolve().pairPotentialRange()); - // Units - // J = kg m2 s-2 --> 10 J = g Ang2 ps-2 - // If ke is in units of [g mol-1 Angstroms2 ps-2] then must use kb in units of 10 J mol-1 K-1 (= 0.8314462) - const auto kb = 0.8314462; - - // Print argument/parameter summary - Messenger::print("MD: Cutoff distance is {}\n", rCut); - Messenger::print("MD: Number of steps = {}\n", nSteps_); - Messenger::print("MD: Timestep type is '{}'\n", timestepType().keyword(timestepType_)); - if (onlyWhenEnergyStable_) - Messenger::print("MD: Only perform MD if target Configuration energies are stable.\n"); - if (trajectoryFrequency_.value_or(0) > 0) - Messenger::print("MD: Trajectory file will be appended every {} step(s).\n", trajectoryFrequency_.value()); - else - Messenger::print("MD: Trajectory file off.\n"); - if (capForces_) - Messenger::print("MD: Forces will be capped to {:10.3e} kJ/mol per atom per axis.\n", maxForce / 100.0); - if (energyFrequency_.value_or(0) > 0) - Messenger::print("MD: Energy will be calculated every {} step(s).\n", energyFrequency_.value()); - else - Messenger::print("MD: Energy will be not be calculated.\n"); - if (outputFrequency_.value_or(0) > 0) - Messenger::print("MD: Summary will be written every {} step(s).\n", outputFrequency_.value()); - else - Messenger::print("MD: Summary will not be written.\n"); - if (!restrictToSpecies_.empty()) - Messenger::print("MD: Calculation will be restricted to species: {}\n", - joinStrings(restrictToSpecies_, " ", [](const auto &sp) { return sp->name(); })); - Messenger::print("\n"); - - if (onlyWhenEnergyStable_) - { - auto stabilityResult = - EnergyModule::checkStability(moduleContext.dissolve().processingModuleData(), targetConfiguration_); - if (stabilityResult == EnergyModule::NotAssessable) - { - return ExecutionResult::Failed; - } - else if (stabilityResult == EnergyModule::EnergyUnstable) - { - Messenger::print("Skipping MD for Configuration '{}'.\n", targetConfiguration_->niceName()); - return ExecutionResult::NotExecuted; - } - } - - // Get temperature from Configuration - const auto temperature = targetConfiguration_->temperature(); - - // Create arrays - std::vector mass(targetConfiguration_->nAtoms(), 0.0); - std::vector> fBound(targetConfiguration_->nAtoms()), fUnbound(targetConfiguration_->nAtoms()), - accelerations(targetConfiguration_->nAtoms()); - - // Variables - auto nCapped = 0; - auto &atoms = targetConfiguration_->atoms(); - double tInstant, ke, tScale, peInter, peIntra; - - // Determine target molecules from the restrictedSpecies vector (if any) - std::vector targetMolecules; - std::vector free(targetConfiguration_->nAtoms(), 0); - if (restrictToSpecies_.empty()) - { - std::fill(free.begin(), free.end(), 1); - } - else - for (const auto &mol : targetConfiguration_->molecules()) - if (std::find(restrictToSpecies_.begin(), restrictToSpecies_.end(), mol->species()) != restrictToSpecies_.end()) - { - targetMolecules.push_back(mol.get()); - auto offset = mol->globalAtomOffset(); - std::fill(free.begin() + offset, free.begin() + offset + mol->atoms().size(), 1); - } - - /* - * Calculation Begins - */ - - // Initialise the random number buffer for all processes - RandomBuffer randomBuffer(moduleContext.processPool(), ProcessPool::PoolProcessesCommunicator); - - // Read in or assign random velocities - auto [velocities, status] = moduleContext.dissolve().processingModuleData().realiseIf>>( - fmt::format("{}//Velocities", targetConfiguration_->niceName()), name(), GenericItem::InRestartFileFlag); - if ((status == GenericItem::ItemStatus::Created || randomVelocities_ || - velocities.size() != targetConfiguration_->nAtoms()) && - !intramolecularForcesOnly_) - { - // Show warning message on array size mismatch - if (status != GenericItem::ItemStatus::Created && velocities.size() != targetConfiguration_->nAtoms()) - Messenger::warn( - "Size of existing velocities array doesn't match the current configuration size - they will be ignored."); - - Messenger::print("Random initial velocities will be assigned.\n"); - velocities.resize(targetConfiguration_->nAtoms(), Vec3()); - for (auto &&[v, iFree] : zip(velocities, free)) - { - if (iFree) - v.set(exp(randomBuffer.random() - 0.5), exp(randomBuffer.random() - 0.5), exp(randomBuffer.random() - 0.5)); - else - v.zero(); - v /= sqrt(TWOPI); - } - } - else if (intramolecularForcesOnly_) - { - Messenger::print("Only intramolecular forces will be calculated, so velocities will be zeroes.\n"); - velocities.resize(targetConfiguration_->nAtoms(), Vec3()); - std::fill(velocities.begin(), velocities.end(), Vec3()); - } - else - { - Messenger::print("Existing velocities will be used.\n"); - } - - Messenger::print("\n"); - - // Store atomic masses for future use - for (auto &&[i, m] : zip(atoms, mass)) - m = AtomicMass::mass(i.speciesAtom()->Z()); - - // Calculate total velocity and mass over all atoms - Vec3 vCom; - auto massSum = 0.0; - for (auto &&[v, m, iFree] : zip(velocities, mass, free)) - { - if (!iFree) - continue; - vCom += v * m; - massSum += m; - } - - // Finalise initial velocities (unless considering intramolecular forces only) - if (!intramolecularForcesOnly_) - { - // Remove any velocity shift, and re-zero velocities on fixed atoms - vCom /= massSum; - std::transform(velocities.begin(), velocities.end(), velocities.begin(), [vCom](auto vel) { return vel - vCom; }); - for (auto &&[v, iFree] : zip(velocities, free)) - if (!iFree) - v.zero(); - - // Calculate instantaneous temperature - ke = 0.0; - for (auto &&[m, v] : zip(mass, velocities)) - ke += 0.5 * m * v.dp(v); - tInstant = ke * 2.0 / (3.0 * atoms.size() * kb); - - // Rescale velocities for desired temperature - tScale = sqrt(temperature / tInstant); - std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto v) { return v * tScale; }); - } - - // Open trajectory file (if requested) - LineParser trajParser; - if (trajectoryFrequency_.value_or(0) > 0) - { - std::string trajectoryFile = fmt::format("{}.md.xyz", targetConfiguration_->name()); - if (moduleContext.processPool().isMaster()) - { - if ((!trajParser.appendOutput(trajectoryFile)) || (!trajParser.isFileGoodForWriting())) - { - Messenger::error("Failed to open MD trajectory output file '{}'.\n", trajectoryFile); - moduleContext.processPool().decideFalse(); - return ExecutionResult::Failed; - } - moduleContext.processPool().decideTrue(); - } - else if (!moduleContext.processPool().decision()) - { - return ExecutionResult::Failed; - } - } - - // Write header - if (outputFrequency_.value_or(0) > 0) - { - Messenger::print(" Energies (kJ/mol)\n"); - Messenger::print(" Step T(K) Kinetic Inter Intra Total " - "deltaT(ps)\n"); - } - - // Start a timer - Timer timer, commsTimer(false); - - // If we're not using a fixed timestep the forces need to be available immediately - if (timestepType_ != TimestepType::Fixed) - { - // Zero force arrays - std::fill(fUnbound.begin(), fUnbound.end(), Vec3()); - std::fill(fBound.begin(), fBound.end(), Vec3()); - - if (targetMolecules.empty()) - ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, - moduleContext.dissolve().potentialMap(), - intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull - : ForcesModule::ForceCalculationType::Full, - fUnbound, fBound, commsTimer); - else - ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, targetMolecules, - moduleContext.dissolve().potentialMap(), - intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull - : ForcesModule::ForceCalculationType::Full, - fUnbound, fBound, commsTimer); - - // Must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) - std::transform(fUnbound.begin(), fUnbound.end(), fUnbound.begin(), [](auto f) { return f * 100.0; }); - std::transform(fBound.begin(), fBound.end(), fBound.begin(), [](auto f) { return f * 100.0; }); - - // Check for suitable timestep - if (!determineTimeStep(timestepType_, fixedTimestep_, fUnbound, fBound)) - { - Messenger::print("Forces are currently too high for MD to proceed. Skipping this run.\n"); - return ExecutionResult::NotExecuted; - } - } - - // Ready to do MD propagation of system - auto step = 1; - for (step = 1; step <= nSteps_; ++step) - { - // Get timestep - auto optDT = determineTimeStep(timestepType_, fixedTimestep_, fUnbound, fBound); - if (!optDT) - { - Messenger::warn("A reasonable timestep could not be determined. Stopping evolution.\n"); - break; - } - auto dT = *optDT; - auto deltaTSq = dT * dT; - - // Velocity Verlet first stage (A) - // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 - // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt - // B: a(t+dt) = F(t+dt)/m - // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt - for (auto &&[i, v, a] : zip(atoms, velocities, accelerations)) - { - // Propagate positions (by whole step)... - i.translateCoordinates(v * dT + a * 0.5 * deltaTSq); - - // ...velocities (by half step)... - v += a * 0.5 * dT; - } - - // Update Atom locations - targetConfiguration_->updateAtomLocations(); - - // Zero force arrays - std::fill(fUnbound.begin(), fUnbound.end(), Vec3()); - std::fill(fBound.begin(), fBound.end(), Vec3()); - - // Calculate forces - must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) - if (targetMolecules.empty()) - ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, - moduleContext.dissolve().potentialMap(), - intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull - : ForcesModule::ForceCalculationType::Full, - fUnbound, fBound, commsTimer); - else - ForcesModule::totalForces(moduleContext.processPool(), targetConfiguration_, targetMolecules, - moduleContext.dissolve().potentialMap(), - intramolecularForcesOnly_ ? ForcesModule::ForceCalculationType::IntraMolecularFull - : ForcesModule::ForceCalculationType::Full, - fUnbound, fBound, commsTimer); - std::transform(fUnbound.begin(), fUnbound.end(), fUnbound.begin(), [](auto f) { return f * 100.0; }); - std::transform(fBound.begin(), fBound.end(), fBound.begin(), [](auto f) { return f * 100.0; }); - - // Cap forces - if (capForces_) - nCapped = capForces(maxForce, fUnbound, fBound); - - // Velocity Verlet second stage (B) and velocity scaling - // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 - // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt - // B: a(t+dt) = F(t+dt)/m - // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt - ke = 0.0; - for (auto &&[f1, f2, v, a, m] : zip(fUnbound, fBound, velocities, accelerations, mass)) - { - // Determine new accelerations - a = (f1 + f2) / m; - - // ..and finally velocities again (by second half-step) - v += a * 0.5 * dT; - - ke += 0.5 * m * v.dp(v); - } - - // Rescale velocities for desired temperature - tInstant = ke * 2.0 / (3.0 * targetConfiguration_->nAtoms() * kb); - tScale = sqrt(temperature / tInstant); - std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto &v) { return v * tScale; }); - - // Convert ke from 10J/mol to kJ/mol - ke *= 0.01; - - // Write step summary? - if (outputFrequency_ && (step == 1 || (step % outputFrequency_.value() == 0))) - { - // Include total energy term? - if (energyFrequency_ && (step % energyFrequency_.value() == 0)) - { - peInter = EnergyModule::interAtomicEnergy(moduleContext.processPool(), targetConfiguration_, - moduleContext.dissolve().potentialMap()); - peIntra = EnergyModule::intraMolecularEnergy(moduleContext.processPool(), targetConfiguration_, - moduleContext.dissolve().potentialMap()); - Messenger::print(" {:<10d} {:10.3e} {:10.3e} {:10.3e} {:10.3e} {:10.3e} {:10.3e}\n", step, - tInstant, ke, peInter, peIntra, ke + peIntra + peInter, dT); - } - else - Messenger::print(" {:<10d} {:10.3e} {:10.3e} {:10.3e}\n", step, - tInstant, ke, dT); - } - - // Save trajectory frame - if (trajectoryFrequency_ && (step % trajectoryFrequency_.value() == 0)) - { - if (moduleContext.processPool().isMaster()) - { - // Write number of atoms - trajParser.writeLineF("{}\n", targetConfiguration_->nAtoms()); - - // Construct and write header - std::string header = fmt::format("Step {} of {}, T = {:10.3e}, ke = {:10.3e}", step, nSteps_, tInstant, ke); - if (energyFrequency_ && (step % energyFrequency_.value() == 0)) - header += fmt::format(", inter = {:10.3e}, intra = {:10.3e}, tot = {:10.3e}", peInter, peIntra, - ke + peInter + peIntra); - if (!trajParser.writeLine(header)) - { - moduleContext.processPool().decideFalse(); - return ExecutionResult::Failed; - } - - // Write Atoms - for (const auto &i : atoms) - { - if (!trajParser.writeLineF("{:<3} {:10.3f} {:10.3f} {:10.3f}\n", Elements::symbol(i.speciesAtom()->Z()), - i.r().x, i.r().y, i.r().z)) - { - moduleContext.processPool().decideFalse(); - return ExecutionResult::Failed; - } - } - - moduleContext.processPool().decideTrue(); - } - else if (!moduleContext.processPool().decision()) - { - return ExecutionResult::Failed; - } - } - } - timer.stop(); - - // Close trajectory file - if (trajectoryFrequency_.value_or(0) > 0 && moduleContext.processPool().isMaster()) - trajParser.closeFiles(); - - if (capForces_) - Messenger::print("A total of {} forces were capped over the course of the dynamics ({:9.3e} per step).\n", nCapped, - double(nCapped) / nSteps_); - Messenger::print("{} steps performed ({} work, {} comms)\n", step - 1, timer.totalTimeString(), - commsTimer.totalTimeString()); - - // Increment configuration changeCount - if (step > 1) - targetConfiguration_->incrementContentsVersion(); /* * Calculation End From f92d0ae9c22e49b94b252d218787e3342fbe9993 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Thu, 16 Nov 2023 13:38:08 +0000 Subject: [PATCH 05/21] CONTROL file formatting implemented --- src/io/export/dlPolyControl.cpp | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 3a041b4cb0..243c561152 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -22,11 +22,36 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) { // Export title - if (!parser.writeLineF("{} @ {}\n", cfg->name(), cfg->contentsVersion())) + if (!parser.writeLineF("title {} @ {}\n\n", cfg->name(), cfg->contentsVersion())) return false; - if (!parser.writeLineF("{}\n", nSteps)) + if (!parser.writeLineF("io_file_config CONFIG\n")) + return false; + if (!parser.writeLineF("io_file_field FIELD\n")) + return false; + if (!parser.writeLineF("io_file_statis STATIS\n")) + return false; + if (!parser.writeLineF("io_file_revive REVIVE\n")) + return false; + if (!parser.writeLineF("io_file_revcon REVCON\n")) + return false; + if (!parser.writeLineF("temperature {} K\n", cfg->temperature())) + return false; + if (!parser.writeLineF("print_frequency {} steps\n", energyFrequency.value())) + return false; + if (!parser.writeLineF("stats_frequency {} steps\n", outputFrequency.value())) + return false; + if (!parser.writeLineF("cutoff {} ang\n", cutoffDistance.value())) + return false; + if (capForces && !parser.writeLineF("equilibration_force_cap {}\n", capForcesAt)) + return false; + if (!parser.writeLineF("time_run {} steps\n", nSteps)) + return false; + if (!parser.writeLineF("{}\n", fixedTimestep)) + return false; + if (!parser.writeLineF("{}\n", randomVelocities)) + return false; + if (!parser.writeLineF("{}\n", trajectoryFrequency.value())) return false; - return true; } From 3c9b38661120767f16e68eaeddd1cb5af45fec3a Mon Sep 17 00:00:00 2001 From: jounaidr Date: Fri, 17 Nov 2023 13:49:26 +0000 Subject: [PATCH 06/21] timestepType removed and boolean timestepVariable implemented --- src/io/export/dlPolyControl.cpp | 26 ++- src/io/export/dlPolyControl.h | 4 +- src/modules/dlPoly/dlPoly.cpp | 20 +- src/modules/dlPoly/dlPoly.h | 22 +- src/modules/dlPoly/functions.cpp | 360 +++++++++++++++---------------- src/modules/dlPoly/process.cpp | 1 + 6 files changed, 208 insertions(+), 225 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 243c561152..694848e07a 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -19,7 +19,7 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi // Export DlPolyControl as CONTROL -bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) +bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) { // Export title if (!parser.writeLineF("title {} @ {}\n\n", cfg->name(), cfg->contentsVersion())) @@ -46,18 +46,27 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("time_run {} steps\n", nSteps)) return false; - if (!parser.writeLineF("{}\n", fixedTimestep)) - return false; - if (!parser.writeLineF("{}\n", randomVelocities)) - return false; - if (!parser.writeLineF("{}\n", trajectoryFrequency.value())) - return false; + if (timestepVariable) + { + if (!parser.writeLineF("timestep_variable ON\n")) + return false; + } else { + if (!parser.writeLineF("timestep_variable OFF\n")) + return false; + } + + //if (!parser.writeLineF("{}\n", fixedTimestep)) + // return false; + //if (!parser.writeLineF("{}\n", randomVelocities)) + // return false; + //if (!parser.writeLineF("{}\n", trajectoryFrequency.value())) + // return false; return true; } // Export DlPolyControl using current filename and format -bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) +bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) { // Open the file LineParser parser; @@ -77,6 +86,7 @@ bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForce capForces, capForcesAt, cutoffDistance, + timestepVariable, fixedTimestep, energyFrequency, nSteps, diff --git a/src/io/export/dlPolyControl.h b/src/io/export/dlPolyControl.h index 8b41bda0fb..fbb41d7da6 100644 --- a/src/io/export/dlPolyControl.h +++ b/src/io/export/dlPolyControl.h @@ -39,11 +39,11 @@ class DlPolyControlExportFileFormat : public FileAndFormat */ private: // Export Configuration as DL_POLY CONTROL - bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); + bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); public: // Export Configuration using current filename and format - bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); + bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); }; diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 311421bd5d..7e5926026f 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -20,8 +20,9 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.setOrganisation("Options", "Simulation"); keywords_.add("NSteps", "Number of DlPoly steps to perform", nSteps_, 1); - keywords_.add>("Timestep", "Timestep type to use in calculation", timestepType_, - DlPolyModule::timestepType()); + keywords_.add("VariableTimestep", + "Whether a variable timestep should be used, determined from the maximal force vector", + timestepVariable_); keywords_.add("DeltaT", "Fixed timestep (ps) to use in DlPoly simulation", fixedTimestep_, 0.0); keywords_.add("RandomVelocities", "Whether random velocities should always be assigned before beginning DlPoly simulation", @@ -42,17 +43,4 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) "CutoffDistance", "Interatomic cutoff distance to use for energy calculation (0.0 to use pair potential range)", cutoffDistance_, 0.0, std::nullopt, 0.1, "Use PairPotential Range"); - // Deprecated - static bool deprecatedBool_{false}; - keywords_.addDeprecated("VariableTimestep", - "Whether a variable timestep should be used, determined from the maximal force vector", - deprecatedBool_); -} - -// Return enum options for TimestepType -EnumOptions DlPolyModule::timestepType() -{ - return EnumOptions( - "TimestepType", - {{TimestepType::Fixed, "Fixed"}, {TimestepType::Variable, "Variable"}, {TimestepType::Automatic, "Auto"}}); -} +} \ No newline at end of file diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 4c67127a7e..8e569ba756 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -21,22 +21,10 @@ class DlPolyModule : public Module /* * Definition */ - public: - // Timestep Type - enum class TimestepType - { - Fixed, - Variable, - Automatic - }; - // Return enum options for TimestepType - static EnumOptions timestepType(); - + + private: // Filename and format for CONTROL export DlPolyControlExportFileFormat dlPolyControlFormat_; - - - private: // Target configurations Configuration *targetConfiguration_{nullptr}; // Control whether atomic forces are capped every step @@ -46,7 +34,7 @@ class DlPolyModule : public Module // Interatomic cutoff distance to employ std::optional cutoffDistance_; // Timestep type to employ - TimestepType timestepType_{TimestepType::Automatic}; + bool timestepVariable_{false}; // Fixed timestep (ps) to use in MD simulation double fixedTimestep_{5.0e-4}; // Frequency at which to calculate total system energy @@ -66,10 +54,6 @@ class DlPolyModule : public Module private: // Cap forces in Configuration static int capForces(double maxForceSq, std::vector> &fInter, std::vector> &fIntra); - // Determine timestep to use - static std::optional determineTimeStep(TimestepType timestepType, double requestedTimeStep, - const std::vector> &fInter, - const std::vector> &fIntra); public: // Evolve Species coordinates, returning new coordinates diff --git a/src/modules/dlPoly/functions.cpp b/src/modules/dlPoly/functions.cpp index 93e9e1b262..95fd11cb3d 100644 --- a/src/modules/dlPoly/functions.cpp +++ b/src/modules/dlPoly/functions.cpp @@ -1,183 +1,183 @@ // SPDX-License-Identifier: GPL-3.0-or-later // Copyright (c) 2023 Team Dissolve and contributors -#include "classes/configuration.h" -#include "classes/species.h" -#include "data/atomicMasses.h" -#include "modules/forces/forces.h" -#include "modules/dlPoly/dlPoly.h" - -// Cap forces in Configuration -int DlPolyModule::capForces(double maxForce, std::vector> &fInter, std::vector> &fIntra) -{ - double fMag; - const auto maxForceSq = maxForce * maxForce; - auto nCapped = 0; - for (auto &&[inter, intra] : zip(fInter, fIntra)) - { - fMag = (inter + intra).magnitudeSq(); - if (fMag < maxForceSq) - continue; - - fMag = maxForce / sqrt(fMag); - inter *= fMag; - intra *= fMag; - - ++nCapped; - } - - return nCapped; -} - -// Determine timestep to use -std::optional DlPolyModule::determineTimeStep(TimestepType timestepType, double requestedTimeStep, - const std::vector> &fInter, - const std::vector> &fIntra) -{ - if (timestepType == TimestepType::Fixed) - return requestedTimeStep; - - // Simple variable timestep - if (timestepType == TimestepType::Variable) - { - auto absFMax = 0.0; - for (auto &&[inter, intra] : zip(fInter, fIntra)) - absFMax = std::max(absFMax, (inter + intra).absMax()); - - return 1.0 / absFMax; - } - - // Automatic timestep determination, using maximal interatomic force to guide the timestep up to the current fixed timestep - // value - auto absFMaxInter = - std::max_element(fInter.begin(), fInter.end(), [](auto &left, auto &right) { return left.absMax() < right.absMax(); }) - ->absMax(); - - auto deltaT = 100.0 / absFMaxInter; - if (deltaT < (requestedTimeStep / 100.0)) - return {}; - return deltaT > requestedTimeStep ? requestedTimeStep : deltaT; -} - -// Evolve Species coordinates, returning new coordinates -std::vector> DlPolyModule::evolve(const ProcessPool &procPool, const PotentialMap &potentialMap, const Species *sp, - double temperature, int nSteps, double maxDeltaT, - const std::vector> &rInit, std::vector> &velocities) -{ - assert(sp); - assert(sp->nAtoms() == velocities.size()); - - // Create arrays - std::vector mass(sp->nAtoms(), 0.0); - std::vector> fInter(sp->nAtoms()), fIntra(sp->nAtoms()), accelerations(sp->nAtoms()); - - // Variables - auto &atoms = sp->atoms(); - double tInstant, ke, tScale; - - // Units - // J = kg m2 s-2 --> 10 J = g Ang2 ps-2 - // If ke is in units of [g mol-1 Angstroms2 ps-2] then must use kb in units of 10 J mol-1 K-1 (= 0.8314462) - const auto kb = 0.8314462; - - // Store atomic masses for future use - std::transform(atoms.begin(), atoms.end(), mass.begin(), [](const auto &atom) { return AtomicMass::mass(atom.Z()); }); - - // Calculate total velocity and mass over all atoms - auto massSum = std::accumulate(mass.begin(), mass.end(), 0.0); - auto vCom = std::transform_reduce(velocities.begin(), velocities.end(), mass.begin(), Vec3()); - - // Remove any velocity shift - vCom /= massSum; - std::transform(velocities.begin(), velocities.end(), velocities.begin(), [vCom](auto vel) { return vel - vCom; }); - - // Calculate instantaneous temperature - ke = std::transform_reduce(mass.begin(), mass.end(), velocities.begin(), 0.0, std::plus<>(), - [](const auto m, const auto &v) { return 0.5 * m * v.dp(v); }); - tInstant = ke * 2.0 / (3.0 * atoms.size() * kb); - - // Rescale velocities for desired temperature - tScale = sqrt(temperature / tInstant); - std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto v) { return v * tScale; }); - - // Zero force arrays - std::fill(fInter.begin(), fInter.end(), Vec3()); - std::fill(fIntra.begin(), fIntra.end(), Vec3()); - - ForcesModule::totalForces(procPool, sp, potentialMap, ForcesModule::ForceCalculationType::Full, fInter, fIntra, rInit); - - // Must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) - std::transform(fInter.begin(), fInter.end(), fInter.begin(), [](auto f) { return f * 100.0; }); - std::transform(fIntra.begin(), fIntra.end(), fIntra.begin(), [](auto f) { return f * 100.0; }); - - // Check for suitable timestep - if (!determineTimeStep(TimestepType::Automatic, maxDeltaT, fInter, fIntra)) - { - Messenger::print("Forces are currently too high for species MD to proceed. Try decreasing the maximum timestep.\n"); - return rInit; - } - - // Copy coordinates ready for propagation - auto rNew = rInit; - - // Ready to do MD propagation of the species - for (auto step = 1; step <= nSteps; ++step) - { - // Get timestep - auto optDT = determineTimeStep(TimestepType::Automatic, maxDeltaT, fInter, fIntra); - if (!optDT) - { - Messenger::warn("A reasonable timestep could not be determined. Stopping evolution.\n"); - break; - } - auto dT = *optDT; - auto deltaTSq = dT * dT; - - // Velocity Verlet first stage (A) - // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 - // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt - // B: a(t+dt) = F(t+dt)/m - // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt - for (auto &&[r, v, a] : zip(rNew, velocities, accelerations)) - { - // Propagate positions (by whole step)... - r += v * dT + a * 0.5 * deltaTSq; - - // ...velocities (by half step)... - v += a * 0.5 * dT; - } - - // Zero force arrays - std::fill(fInter.begin(), fInter.end(), Vec3()); - std::fill(fIntra.begin(), fIntra.end(), Vec3()); - - // Calculate forces - must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) - ForcesModule::totalForces(procPool, sp, potentialMap, ForcesModule::ForceCalculationType::Full, fInter, fIntra, rNew); - std::transform(fInter.begin(), fInter.end(), fInter.begin(), [](auto f) { return f * 100.0; }); - std::transform(fIntra.begin(), fIntra.end(), fIntra.begin(), [](auto f) { return f * 100.0; }); - - // Velocity Verlet second stage (B) and velocity scaling - // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 - // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt - // B: a(t+dt) = F(t+dt)/m - // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt - ke = 0.0; - for (auto &&[f1, f2, v, a, m] : zip(fInter, fIntra, velocities, accelerations, mass)) - { - // Determine new accelerations - a = (f1 + f2) / m; - - // ..and finally velocities again (by second half-step) - v += a * 0.5 * dT; - - ke += 0.5 * m * v.dp(v); - } - - // Rescale velocities for desired temperature - tInstant = ke * 2.0 / (3.0 * sp->nAtoms() * kb); - tScale = sqrt(temperature / tInstant); - std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto &v) { return v * tScale; }); - } - - return rNew; -} +//#include "classes/configuration.h" +//#include "classes/species.h" +//#include "data/atomicMasses.h" +//#include "modules/forces/forces.h" +//#include "modules/dlPoly/dlPoly.h" +// +//// Cap forces in Configuration +//int DlPolyModule::capForces(double maxForce, std::vector> &fInter, std::vector> &fIntra) +//{ +// double fMag; +// const auto maxForceSq = maxForce * maxForce; +// auto nCapped = 0; +// for (auto &&[inter, intra] : zip(fInter, fIntra)) +// { +// fMag = (inter + intra).magnitudeSq(); +// if (fMag < maxForceSq) +// continue; +// +// fMag = maxForce / sqrt(fMag); +// inter *= fMag; +// intra *= fMag; +// +// ++nCapped; +// } +// +// return nCapped; +//} +// +//// Determine timestep to use +//std::optional DlPolyModule::determineTimeStep(TimestepType timestepType, double requestedTimeStep, +// const std::vector> &fInter, +// const std::vector> &fIntra) +//{ +// if (timestepType == TimestepType::Fixed) +// return requestedTimeStep; +// +// // Simple variable timestep +// if (timestepType == TimestepType::Variable) +// { +// auto absFMax = 0.0; +// for (auto &&[inter, intra] : zip(fInter, fIntra)) +// absFMax = std::max(absFMax, (inter + intra).absMax()); +// +// return 1.0 / absFMax; +// } +// +// // Automatic timestep determination, using maximal interatomic force to guide the timestep up to the current fixed timestep +// // value +// auto absFMaxInter = +// std::max_element(fInter.begin(), fInter.end(), [](auto &left, auto &right) { return left.absMax() < right.absMax(); }) +// ->absMax(); +// +// auto deltaT = 100.0 / absFMaxInter; +// if (deltaT < (requestedTimeStep / 100.0)) +// return {}; +// return deltaT > requestedTimeStep ? requestedTimeStep : deltaT; +//} +// +//// Evolve Species coordinates, returning new coordinates +//std::vector> DlPolyModule::evolve(const ProcessPool &procPool, const PotentialMap &potentialMap, const Species *sp, +// double temperature, int nSteps, double maxDeltaT, +// const std::vector> &rInit, std::vector> &velocities) +//{ +// assert(sp); +// assert(sp->nAtoms() == velocities.size()); +// +// // Create arrays +// std::vector mass(sp->nAtoms(), 0.0); +// std::vector> fInter(sp->nAtoms()), fIntra(sp->nAtoms()), accelerations(sp->nAtoms()); +// +// // Variables +// auto &atoms = sp->atoms(); +// double tInstant, ke, tScale; +// +// // Units +// // J = kg m2 s-2 --> 10 J = g Ang2 ps-2 +// // If ke is in units of [g mol-1 Angstroms2 ps-2] then must use kb in units of 10 J mol-1 K-1 (= 0.8314462) +// const auto kb = 0.8314462; +// +// // Store atomic masses for future use +// std::transform(atoms.begin(), atoms.end(), mass.begin(), [](const auto &atom) { return AtomicMass::mass(atom.Z()); }); +// +// // Calculate total velocity and mass over all atoms +// auto massSum = std::accumulate(mass.begin(), mass.end(), 0.0); +// auto vCom = std::transform_reduce(velocities.begin(), velocities.end(), mass.begin(), Vec3()); +// +// // Remove any velocity shift +// vCom /= massSum; +// std::transform(velocities.begin(), velocities.end(), velocities.begin(), [vCom](auto vel) { return vel - vCom; }); +// +// // Calculate instantaneous temperature +// ke = std::transform_reduce(mass.begin(), mass.end(), velocities.begin(), 0.0, std::plus<>(), +// [](const auto m, const auto &v) { return 0.5 * m * v.dp(v); }); +// tInstant = ke * 2.0 / (3.0 * atoms.size() * kb); +// +// // Rescale velocities for desired temperature +// tScale = sqrt(temperature / tInstant); +// std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto v) { return v * tScale; }); +// +// // Zero force arrays +// std::fill(fInter.begin(), fInter.end(), Vec3()); +// std::fill(fIntra.begin(), fIntra.end(), Vec3()); +// +// ForcesModule::totalForces(procPool, sp, potentialMap, ForcesModule::ForceCalculationType::Full, fInter, fIntra, rInit); +// +// // Must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) +// std::transform(fInter.begin(), fInter.end(), fInter.begin(), [](auto f) { return f * 100.0; }); +// std::transform(fIntra.begin(), fIntra.end(), fIntra.begin(), [](auto f) { return f * 100.0; }); +// +// // Check for suitable timestep +// if (!determineTimeStep(TimestepType::Automatic, maxDeltaT, fInter, fIntra)) +// { +// Messenger::print("Forces are currently too high for species MD to proceed. Try decreasing the maximum timestep.\n"); +// return rInit; +// } +// +// // Copy coordinates ready for propagation +// auto rNew = rInit; +// +// // Ready to do MD propagation of the species +// for (auto step = 1; step <= nSteps; ++step) +// { +// // Get timestep +// auto optDT = determineTimeStep(TimestepType::Automatic, maxDeltaT, fInter, fIntra); +// if (!optDT) +// { +// Messenger::warn("A reasonable timestep could not be determined. Stopping evolution.\n"); +// break; +// } +// auto dT = *optDT; +// auto deltaTSq = dT * dT; +// +// // Velocity Verlet first stage (A) +// // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 +// // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt +// // B: a(t+dt) = F(t+dt)/m +// // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt +// for (auto &&[r, v, a] : zip(rNew, velocities, accelerations)) +// { +// // Propagate positions (by whole step)... +// r += v * dT + a * 0.5 * deltaTSq; +// +// // ...velocities (by half step)... +// v += a * 0.5 * dT; +// } +// +// // Zero force arrays +// std::fill(fInter.begin(), fInter.end(), Vec3()); +// std::fill(fIntra.begin(), fIntra.end(), Vec3()); +// +// // Calculate forces - must multiply by 100.0 to convert from kJ/mol to 10J/mol (our internal MD units) +// ForcesModule::totalForces(procPool, sp, potentialMap, ForcesModule::ForceCalculationType::Full, fInter, fIntra, rNew); +// std::transform(fInter.begin(), fInter.end(), fInter.begin(), [](auto f) { return f * 100.0; }); +// std::transform(fIntra.begin(), fIntra.end(), fIntra.begin(), [](auto f) { return f * 100.0; }); +// +// // Velocity Verlet second stage (B) and velocity scaling +// // A: r(t+dt) = r(t) + v(t)*dt + 0.5*a(t)*dt**2 +// // A: v(t+dt/2) = v(t) + 0.5*a(t)*dt +// // B: a(t+dt) = F(t+dt)/m +// // B: v(t+dt) = v(t+dt/2) + 0.5*a(t+dt)*dt +// ke = 0.0; +// for (auto &&[f1, f2, v, a, m] : zip(fInter, fIntra, velocities, accelerations, mass)) +// { +// // Determine new accelerations +// a = (f1 + f2) / m; +// +// // ..and finally velocities again (by second half-step) +// v += a * 0.5 * dT; +// +// ke += 0.5 * m * v.dp(v); +// } +// +// // Rescale velocities for desired temperature +// tInstant = ke * 2.0 / (3.0 * sp->nAtoms() * kb); +// tScale = sqrt(temperature / tInstant); +// std::transform(velocities.begin(), velocities.end(), velocities.begin(), [tScale](auto &v) { return v * tScale; }); +// } +// +// return rNew; +//} diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 5d84d05bf9..815a43c80f 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -31,6 +31,7 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) capForces_, capForcesAt_, cutoffDistance_, + timestepVariable_, fixedTimestep_, energyFrequency_, nSteps_, From 644ac207640379cb493546af6fc396262d183cc8 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Fri, 24 Nov 2023 14:05:41 +0000 Subject: [PATCH 07/21] =?UTF-8?q?add=20coul=20methods,=20precision=20and?= =?UTF-8?q?=20trajectory=20key=C3=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/io/export/dlPolyControl.cpp | 34 +++++++++++++++------ src/io/export/dlPolyControl.h | 4 +-- src/io/export/dlPolyField.cpp | 54 +++++++++++++++++++++++++++++++++ src/io/export/dlPolyField.h | 50 ++++++++++++++++++++++++++++++ src/modules/dlPoly/dlPoly.cpp | 7 ++++- src/modules/dlPoly/dlPoly.h | 29 ++++++++++++++++++ src/modules/dlPoly/process.cpp | 21 ++++++++++++- 7 files changed, 186 insertions(+), 13 deletions(-) create mode 100644 src/io/export/dlPolyField.cpp create mode 100644 src/io/export/dlPolyField.h diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 694848e07a..8132a23e18 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -14,12 +14,12 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi : FileAndFormat(formats_, filename, (int)format) { formats_ = EnumOptions( - "CoordinateExportFileFormat", {{DlPolyControlExportFormat::DLPOLY, "dlpoly", "DL_POLY CONFIG File"}}); + "ControlExportFileFormat", {{DlPolyControlExportFormat::DLPOLY, "dlpoly", "DL_POLY CONFIG File"}}); } // Export DlPolyControl as CONTROL -bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) +bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) { // Export title if (!parser.writeLineF("title {} @ {}\n\n", cfg->name(), cfg->contentsVersion())) @@ -42,6 +42,12 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("cutoff {} ang\n", cutoffDistance.value())) return false; + if (!parser.writeLineF("ensemble nvt\n", cutoffDistance.value())) + return false; + if (!parser.writeLineF("ensemble_method hoover\n", cutoffDistance.value())) + return false; + if (!parser.writeLineF("ensemble_thermostat_coupling 0.01 ps\n", cutoffDistance.value())) + return false; if (capForces && !parser.writeLineF("equilibration_force_cap {}\n", capForcesAt)) return false; if (!parser.writeLineF("time_run {} steps\n", nSteps)) @@ -50,23 +56,30 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati { if (!parser.writeLineF("timestep_variable ON\n")) return false; - } else { - if (!parser.writeLineF("timestep_variable OFF\n")) + } + if (trajectoryFrequency.value_or(0) > 0) + { + if (!parser.writeLineF("traj_calculate ON\n")) + return false; + if (!parser.writeLineF("traj_interval {}\n", trajectoryFrequency.value())) + return false; + if (!parser.writeLineF("traj_key {}\n", trajectoryKey)) return false; } - + if (!parser.writeLineF("coul_method {}\n", coulMethod)) + return false; + if (!parser.writeLineF("coul_precision {}\n", coulPrecision)) + return false; //if (!parser.writeLineF("{}\n", fixedTimestep)) // return false; //if (!parser.writeLineF("{}\n", randomVelocities)) // return false; - //if (!parser.writeLineF("{}\n", trajectoryFrequency.value())) - // return false; return true; } // Export DlPolyControl using current filename and format -bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency) +bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) { // Open the file LineParser parser; @@ -92,7 +105,10 @@ bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForce nSteps, outputFrequency, randomVelocities, - trajectoryFrequency); + trajectoryFrequency, + trajectoryKey, + coulMethod, + coulPrecision); break; default: throw(std::runtime_error(fmt::format("DlPolyControl format '{}' export has not been implemented.\n", diff --git a/src/io/export/dlPolyControl.h b/src/io/export/dlPolyControl.h index fbb41d7da6..15c2fa5478 100644 --- a/src/io/export/dlPolyControl.h +++ b/src/io/export/dlPolyControl.h @@ -39,11 +39,11 @@ class DlPolyControlExportFileFormat : public FileAndFormat */ private: // Export Configuration as DL_POLY CONTROL - bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); + bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); public: // Export Configuration using current filename and format - bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency); + bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); }; diff --git a/src/io/export/dlPolyField.cpp b/src/io/export/dlPolyField.cpp new file mode 100644 index 0000000000..5023fa9fec --- /dev/null +++ b/src/io/export/dlPolyField.cpp @@ -0,0 +1,54 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#include "io/export/dlPolyField.h" +#include "base/lineParser.h" +#include "base/sysFunc.h" +#include "classes/atomType.h" +#include "classes/box.h" +#include "classes/configuration.h" +#include "classes/speciesAtom.h" +#include "data/atomicMasses.h" + +DlPolyControlExportFileFormat::DlPolyFieldExportFileFormat(std::string_view filename, DlPolyFieldExportFormat format) + : FileAndFormat(formats_, filename, (int)format) +{ + formats_ = EnumOptions( + "FieldExportFileFormat", {{DlPolyFieldExportFormat::DLPOLY, "dlpoly", "DL_POLY FIELD File"}}); +} + + +// Export DlPolyField as Field +bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg) +{ + // Export title + + + return true; +} + +// Export DlPolyField using current filename and format +bool DlPolyFieldExportFileFormat::exportData(Configuration *cfg) +{ + // Open the file + LineParser parser; + if (!parser.openOutput(filename_)) + { + parser.closeFiles(); + return false; + } + + // Write data + auto result = false; + switch (formats_.enumerationByIndex(*formatIndex_)) + { + case (DlPolyFieldExportFormat::DLPOLY): + result = exportDLPOLY(parser, cfg); + break; + default: + throw(std::runtime_error(fmt::format("DlPolyField format '{}' export has not been implemented.\n", + formats_.keywordByIndex(*formatIndex_)))); + } + + return true; +} diff --git a/src/io/export/dlPolyField.h b/src/io/export/dlPolyField.h new file mode 100644 index 0000000000..0c9adb2ddc --- /dev/null +++ b/src/io/export/dlPolyField.h @@ -0,0 +1,50 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2023 Team Dissolve and contributors + +#pragma once + +#include "io/fileAndFormat.h" + +// Forward Declarations +class Configuration; + +// Coordinate Export Formats +class DlPolyControlExportFileFormat : public FileAndFormat +{ + public: + // Available DlPolyControl formats + enum class DlPolyControlExportFormat + { + DLPOLY + }; + DlPolyFieldExportFileFormat(std::string_view filename = "", DlPolyFieldExportFormat format = DlPolyFieldExportFormat::DLPOLY); + ~DlPolyFieldExportFileFormat() override = default; + + /* + * Formats + */ + private: + // Format enum options + EnumOptions formats_; + + /* + * Filename / Basename + */ + public: + // Return whether the file must exist + bool fileMustExist() const override { return false; } + + /* + * Export Functions + */ + private: + // Export Configuration as DL_POLY FIELD + bool exportDLPOLY(LineParser &parser, Configuration *cfg); + + public: + // Export Configuration using current filename and format + bool exportData(Configuration *cfg); + +}; + + diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 7e5926026f..653be39b1f 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -8,6 +8,7 @@ #include "keywords/integer.h" #include "keywords/optionalDouble.h" #include "keywords/optionalInt.h" +#include "keywords/enumOptions.h" #include "keywords/speciesVector.h" #include "keywords/fileAndFormat.h" @@ -35,6 +36,11 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) 0, std::nullopt, 5, "Off"); keywords_.add("TrajectoryFrequency", "Write frequency for trajectory file", trajectoryFrequency_, 0, std::nullopt, 5, "Off"); + keywords_.add>( + "TrajectoryKey", "Set trajectory key", trajectoryKey_, DlPolyModule::trajectoryKey()); + keywords_.add>( + "CoulMethod", "Set Coul Method", coulMethod_, DlPolyModule::coulMethod()); + keywords_.add("Coul Precision", "Calculate electrostatics using Fennell damping (Ewald-like) with given precision", coulPrecision_, 0.0); keywords_.setOrganisation("Advanced"); keywords_.add("CapForces", "Control whether atomic forces are capped every step", capForces_); @@ -42,5 +48,4 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.add( "CutoffDistance", "Interatomic cutoff distance to use for energy calculation (0.0 to use pair potential range)", cutoffDistance_, 0.0, std::nullopt, 0.1, "Use PairPotential Range"); - } \ No newline at end of file diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 8e569ba756..136652c492 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -21,6 +21,29 @@ class DlPolyModule : public Module /* * Definition */ + public: + // TrajectoryKey Enum + enum TrajectoryKey + { + pos, + pos_vel, + pos_vel_force + }; + // Return EnumOptions for TrajectoryKey + static EnumOptions trajectoryKey(); + + // CoulMethod Enum + enum CoulMethod + { + off, + spme, + dddp, + pairwise, + reaction_field, + force_shifted + }; + // Return EnumOptions for CoulMethod + static EnumOptions coulMethod(); private: // Filename and format for CONTROL export @@ -47,6 +70,12 @@ class DlPolyModule : public Module bool randomVelocities_{false}; // Write frequency for trajectory file std::optional trajectoryFrequency_; + // Type of target TrajectoryKeySet + DlPolyModule::TrajectoryKey trajectoryKey_{TrajectoryKey::pos}; + // Set method for electrostatics method + DlPolyModule::CoulMethod coulMethod_{CoulMethod::off}; + // Calculate electrostatics with given precision + double coulPrecision_{0.0}; /* * Functions diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 815a43c80f..3068bdf50f 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -11,6 +11,22 @@ #include "modules/forces/forces.h" #include "modules/dlPoly/dlPoly.h" +// Return EnumOptions for TrajectoryKey +EnumOptions DlPolyModule::trajectoryKey() +{ + return EnumOptions( + "TrajectoryKey", + {{TrajectoryKey::pos, "pos"}, {TrajectoryKey::pos_vel, "pos-vel"}, {TrajectoryKey::pos_vel_force, "pos-vel-force"}}); +} + +// Return EnumOptions for CoulMethod +EnumOptions DlPolyModule::coulMethod() +{ + return EnumOptions( + "CoulMethod", + {{CoulMethod::off, "off"}, {CoulMethod::spme, "spme"}, {CoulMethod::dddp, "dddp"}, {CoulMethod::pairwise, "pairwise"}, {CoulMethod::reaction_field, "reaction_field"}, {CoulMethod::force_shifted, "force_shifted"}}); +} + // Run main processing Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) { @@ -37,7 +53,10 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) nSteps_, outputFrequency_, randomVelocities_, - trajectoryFrequency_)) + trajectoryFrequency_, + trajectoryKey().keyword(trajectoryKey_), + coulMethod().keyword(coulMethod_), + coulPrecision_)) { Messenger::print("Export: Failed to export coordinates file '{}'.\n", dlPolyControlFormat_.filename()); moduleContext.processPool().decideFalse(); From 3d8a5def896c989a7c3b19c6d42d8ec619907130 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Tue, 28 Nov 2023 09:55:54 +0000 Subject: [PATCH 08/21] all CONTROL file options implemented --- src/io/export/dlPolyControl.cpp | 16 ++++++++-------- src/io/export/dlPolyControl.h | 4 ++-- src/modules/dlPoly/dlPoly.cpp | 4 ++++ src/modules/dlPoly/dlPoly.h | 17 +++++++++++++++++ src/modules/dlPoly/process.cpp | 13 +++++++++++++ 5 files changed, 44 insertions(+), 10 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 8132a23e18..1a73891011 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -19,7 +19,7 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi // Export DlPolyControl as CONTROL -bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) +bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) { // Export title if (!parser.writeLineF("title {} @ {}\n\n", cfg->name(), cfg->contentsVersion())) @@ -42,11 +42,11 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("cutoff {} ang\n", cutoffDistance.value())) return false; - if (!parser.writeLineF("ensemble nvt\n", cutoffDistance.value())) + if (!parser.writeLineF("ensemble {}\n", ensemble)) return false; if (!parser.writeLineF("ensemble_method hoover\n", cutoffDistance.value())) return false; - if (!parser.writeLineF("ensemble_thermostat_coupling 0.01 ps\n", cutoffDistance.value())) + if (!parser.writeLineF("ensemble_thermostat_coupling {} ps\n", ensembleThermostatCoupling)) return false; if (capForces && !parser.writeLineF("equilibration_force_cap {}\n", capForcesAt)) return false; @@ -70,16 +70,14 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("coul_precision {}\n", coulPrecision)) return false; - //if (!parser.writeLineF("{}\n", fixedTimestep)) - // return false; - //if (!parser.writeLineF("{}\n", randomVelocities)) - // return false; + + Messenger::print("CONFIG {}\n", cfg); return true; } // Export DlPolyControl using current filename and format -bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) +bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) { // Open the file LineParser parser; @@ -99,6 +97,8 @@ bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForce capForces, capForcesAt, cutoffDistance, + ensemble, + ensembleThermostatCoupling, timestepVariable, fixedTimestep, energyFrequency, diff --git a/src/io/export/dlPolyControl.h b/src/io/export/dlPolyControl.h index 15c2fa5478..243b779494 100644 --- a/src/io/export/dlPolyControl.h +++ b/src/io/export/dlPolyControl.h @@ -39,11 +39,11 @@ class DlPolyControlExportFileFormat : public FileAndFormat */ private: // Export Configuration as DL_POLY CONTROL - bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); + bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); public: // Export Configuration using current filename and format - bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); + bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); }; diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 653be39b1f..0a32f33b81 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -48,4 +48,8 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.add( "CutoffDistance", "Interatomic cutoff distance to use for energy calculation (0.0 to use pair potential range)", cutoffDistance_, 0.0, std::nullopt, 0.1, "Use PairPotential Range"); + keywords_.add>( + "Ensemble", " Set ensemble constraints", ensemble_, DlPolyModule::ensemble()); + keywords_.add("EnsembleThermostatCoupling", "Set thermostat relaxation/decorrelation times (use ensemble_thermostat_friction for langevin)", ensembleThermostatCoupling_, 0.0); + } \ No newline at end of file diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 136652c492..05c7787c89 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -44,6 +44,18 @@ class DlPolyModule : public Module }; // Return EnumOptions for CoulMethod static EnumOptions coulMethod(); + + // Ensemble Enum + enum Ensemble + { + NVE, + PMF, + NVT, + NPT, + NST + }; + // Return EnumOptions for Ensemble + static EnumOptions ensemble(); private: // Filename and format for CONTROL export @@ -56,6 +68,10 @@ class DlPolyModule : public Module double capForcesAt_{1.0e7}; // Interatomic cutoff distance to employ std::optional cutoffDistance_; + // Set Ensemble + DlPolyModule::Ensemble ensemble_{Ensemble::NVE}; + // Set Ensemble Thermostat Coupling + double ensembleThermostatCoupling_{0.0}; // Timestep type to employ bool timestepVariable_{false}; // Fixed timestep (ps) to use in MD simulation @@ -76,6 +92,7 @@ class DlPolyModule : public Module DlPolyModule::CoulMethod coulMethod_{CoulMethod::off}; // Calculate electrostatics with given precision double coulPrecision_{0.0}; + /* * Functions diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 3068bdf50f..f3ee9cdc86 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -27,9 +27,20 @@ EnumOptions DlPolyModule::coulMethod() {{CoulMethod::off, "off"}, {CoulMethod::spme, "spme"}, {CoulMethod::dddp, "dddp"}, {CoulMethod::pairwise, "pairwise"}, {CoulMethod::reaction_field, "reaction_field"}, {CoulMethod::force_shifted, "force_shifted"}}); } +// Return EnumOptions for Ensemble +EnumOptions DlPolyModule::ensemble() +{ + return EnumOptions( + "ensemble", + {{Ensemble::NVE, "NVE"}, {Ensemble::PMF, "PMF"}, {Ensemble::NVT, "NVT"}, {Ensemble::NPT, "NPT"}, {Ensemble::NST, "NST"}}); +} + // Run main processing Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) { + + Messenger::print("CONFIG {}\n", targetConfiguration_); + // Check for zero Configuration targets if (!targetConfiguration_) { @@ -47,6 +58,8 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) capForces_, capForcesAt_, cutoffDistance_, + ensemble().keyword(ensemble_), + ensembleThermostatCoupling_, timestepVariable_, fixedTimestep_, energyFrequency_, From d5a29ea5c266ecd1da8822326792aed95fdcdb74 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 6 Dec 2023 16:51:50 +0000 Subject: [PATCH 09/21] field file export implemented, params not working... --- src/io/export/CMakeLists.txt | 2 + src/io/export/dlPolyControl.cpp | 2 - src/io/export/dlPolyField.cpp | 81 +++++++++++++++++++++++++++++++-- src/modules/dlPoly/dlPoly.cpp | 3 +- src/modules/dlPoly/dlPoly.h | 3 ++ src/modules/dlPoly/process.cpp | 16 ++++--- 6 files changed, 95 insertions(+), 12 deletions(-) diff --git a/src/io/export/CMakeLists.txt b/src/io/export/CMakeLists.txt index 36c40c1311..cdda418c81 100644 --- a/src/io/export/CMakeLists.txt +++ b/src/io/export/CMakeLists.txt @@ -10,6 +10,7 @@ add_library( trajectory.cpp values.cpp dlPolyControl.cpp + dlPolyField.cpp coordinates.h data1D.h data2D.h @@ -20,6 +21,7 @@ add_library( trajectory.h values.h dlPolyControl.h + dlPolyField.h ) target_include_directories(export PRIVATE ${PROJECT_SOURCE_DIR}/src ${PROJECT_BINARY_DIR}/src) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 1a73891011..22bb857b4f 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -70,8 +70,6 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("coul_precision {}\n", coulPrecision)) return false; - - Messenger::print("CONFIG {}\n", cfg); return true; } diff --git a/src/io/export/dlPolyField.cpp b/src/io/export/dlPolyField.cpp index 5023fa9fec..a230659af3 100644 --- a/src/io/export/dlPolyField.cpp +++ b/src/io/export/dlPolyField.cpp @@ -10,7 +10,46 @@ #include "classes/speciesAtom.h" #include "data/atomicMasses.h" -DlPolyControlExportFileFormat::DlPolyFieldExportFileFormat(std::string_view filename, DlPolyFieldExportFormat format) + +#include "classes/atom.h" +#include "classes/atomType.h" +#include "classes/atomTypeMix.h" +#include "classes/molecule.h" +#include +#include + +#include "base/serialiser.h" +#include "base/version.h" +#include "classes/atom.h" +#include "classes/atomTypeMix.h" +#include "classes/box.h" +#include "classes/cellArray.h" +#include "classes/molecule.h" +#include "classes/siteStack.h" +#include "io/import/coordinates.h" +#include "items/list.h" +#include "kernels/potentials/base.h" +#include "math/data1D.h" +#include "math/histogram1D.h" +#include "math/interpolator.h" +#include "module/layer.h" +#include "procedure/procedure.h" +#include "templates/vector3.h" +#include +#include +#include +#include + +#include "io/export/coordinates.h" +#include "base/lineParser.h" +#include "base/sysFunc.h" +#include "classes/atomType.h" +#include "classes/box.h" +#include "classes/configuration.h" +#include "classes/speciesAtom.h" +#include "data/atomicMasses.h" + +DlPolyFieldExportFileFormat::DlPolyFieldExportFileFormat(std::string_view filename, DlPolyFieldExportFormat format) : FileAndFormat(formats_, filename, (int)format) { formats_ = EnumOptions( @@ -22,8 +61,44 @@ DlPolyControlExportFileFormat::DlPolyFieldExportFileFormat(std::string_view file bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg) { // Export title - - + if (!parser.writeLineF("title {} @ {}\n", cfg->name(), cfg->contentsVersion())) + return false; + if (!parser.writeLineF("units Kj\n")) + return false; + + const std::vector> speciesPopulations = cfg->speciesPopulations(); + + if (!parser.writeLineF("moleculer types {}\n", speciesPopulations.size())) + return false; + + for (const auto &species : speciesPopulations){ + if (!parser.writeLineF("{}\nnummols {}\natoms {}\n", species.first->name(), species.second, species.first->nAtoms())){ + return false; + } + for (const auto &atom : species.first->atoms()){ + if (!parser.writeLineF("{} {} {} 1 0\n", atom.atomType()->name(), AtomicMass::mass(atom.Z()), atom.charge())) + return false; + } + if (!parser.writeLineF("bonds {}\n", species.first->nBonds())) + return false; + for (const auto &bond : species.first->bonds()){ + //const auto params = bond.interactionParameters(); + //Messenger::error("{}",params); + if (!parser.writeLineF("{} {} {} {} {}\n", BondFunctions::forms().keyword(bond.interactionForm()).substr(0, 4), bond.indexI() + 1, bond.indexJ() + 1), "bond.interactionParameters()[0]", "bond.interactionParameters()[1]") + return false; + } + if (!parser.writeLineF("angles {}\n", species.first->nAngles())) + return false; + for (const auto &angle : species.first->angles()){ + if (!parser.writeLineF("{} {} {} {} {}\n", AngleFunctions::forms().keyword(angle.interactionForm()).substr(0, 4), angle.indexI() + 1, angle.indexJ() + 1), "angle.interactionParameters()[0]", "angle.interactionParameters()[1]") + return false; + } + } + + if (!parser.writeLineF("finish")) + return false; + + return true; } diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 0a32f33b81..8a3a9da5d7 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -17,7 +17,8 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.addTarget("Configuration", "Set target configuration for the module", targetConfiguration_); keywords_.setOrganisation("Options", "File"); - keywords_.add("Format", "File / format for coordinates", dlPolyControlFormat_, "EndFormat"); + keywords_.add("CONTROL", "File / format for CONTROL", dlPolyControlFormat_, "EndFormat"); + keywords_.add("FIELD", "File / format for FIELD", dlPolyFieldFormat_, "EndFormat"); keywords_.setOrganisation("Options", "Simulation"); keywords_.add("NSteps", "Number of DlPoly steps to perform", nSteps_, 1); diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 05c7787c89..754e64ddc9 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -4,6 +4,7 @@ #pragma once #include "io/export/dlPolyControl.h" +#include "io/export/dlPolyField.h" #include "base/enumOptions.h" #include "module/module.h" @@ -60,6 +61,8 @@ class DlPolyModule : public Module private: // Filename and format for CONTROL export DlPolyControlExportFileFormat dlPolyControlFormat_; + // Filename and format for FIELD export + DlPolyFieldExportFileFormat dlPolyFieldFormat_; // Target configurations Configuration *targetConfiguration_{nullptr}; // Control whether atomic forces are capped every step diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index f3ee9cdc86..d6b8ed192c 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -37,10 +37,7 @@ EnumOptions DlPolyModule::ensemble() // Run main processing Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) -{ - - Messenger::print("CONFIG {}\n", targetConfiguration_); - +{ // Check for zero Configuration targets if (!targetConfiguration_) { @@ -48,7 +45,7 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } - // Save DlPoly CONTROL file with props + // Save DlPoly CONTROL/FIELD file with props if (moduleContext.processPool().isMaster()) { Messenger::print("Export: Writing TESTTEST file ({}) for Configuration '{}'...\n", @@ -71,7 +68,14 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) coulMethod().keyword(coulMethod_), coulPrecision_)) { - Messenger::print("Export: Failed to export coordinates file '{}'.\n", dlPolyControlFormat_.filename()); + Messenger::print("Export: Failed to export CONTROL file '{}'.\n", dlPolyControlFormat_.filename()); + moduleContext.processPool().decideFalse(); + return ExecutionResult::Failed; + } + + if (!dlPolyFieldFormat_.exportData(targetConfiguration_)) + { + Messenger::print("Export: Failed to export FIELD file '{}'.\n", dlPolyFieldFormat_.filename()); moduleContext.processPool().decideFalse(); return ExecutionResult::Failed; } From dff8c4ee1e90ae40d0b3d85e6b40fb9020f5f35d Mon Sep 17 00:00:00 2001 From: jounaidr Date: Thu, 14 Dec 2023 18:23:03 +0000 Subject: [PATCH 10/21] field file fully implemented --- src/io/export/dlPolyControl.cpp | 6 +++++- src/io/export/dlPolyField.cpp | 35 +++++++++++++++++++++++++++------ src/io/export/dlPolyField.h | 6 +++--- 3 files changed, 37 insertions(+), 10 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 22bb857b4f..eee2beeb87 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -48,7 +48,7 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("ensemble_thermostat_coupling {} ps\n", ensembleThermostatCoupling)) return false; - if (capForces && !parser.writeLineF("equilibration_force_cap {}\n", capForcesAt)) + if (capForces && !parser.writeLineF("equilibration_force_cap {} k_b.temp/ang\n", capForcesAt)) return false; if (!parser.writeLineF("time_run {} steps\n", nSteps)) return false; @@ -57,6 +57,8 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati if (!parser.writeLineF("timestep_variable ON\n")) return false; } + if (!parser.writeLineF("timestep {} internal_t\n", fixedTimestep)) + return false; if (trajectoryFrequency.value_or(0) > 0) { if (!parser.writeLineF("traj_calculate ON\n")) @@ -70,6 +72,8 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("coul_precision {}\n", coulPrecision)) return false; + if (!parser.writeLineF("vdw_mix_method Lorentz-Berthelot\n")) + return false; return true; } diff --git a/src/io/export/dlPolyField.cpp b/src/io/export/dlPolyField.cpp index a230659af3..3937f7cd0f 100644 --- a/src/io/export/dlPolyField.cpp +++ b/src/io/export/dlPolyField.cpp @@ -82,22 +82,45 @@ bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration if (!parser.writeLineF("bonds {}\n", species.first->nBonds())) return false; for (const auto &bond : species.first->bonds()){ - //const auto params = bond.interactionParameters(); - //Messenger::error("{}",params); - if (!parser.writeLineF("{} {} {} {} {}\n", BondFunctions::forms().keyword(bond.interactionForm()).substr(0, 4), bond.indexI() + 1, bond.indexJ() + 1), "bond.interactionParameters()[0]", "bond.interactionParameters()[1]") + const auto params = bond.interactionParameters(); + if (!parser.writeLineF("{} {} {} {} {}\n", BondFunctions::forms().keyword(bond.interactionForm()).substr(0, 4), bond.indexI() + 1, bond.indexJ() + 1, params[0], params[1])) return false; } if (!parser.writeLineF("angles {}\n", species.first->nAngles())) return false; for (const auto &angle : species.first->angles()){ - if (!parser.writeLineF("{} {} {} {} {}\n", AngleFunctions::forms().keyword(angle.interactionForm()).substr(0, 4), angle.indexI() + 1, angle.indexJ() + 1), "angle.interactionParameters()[0]", "angle.interactionParameters()[1]") + const auto params = angle.interactionParameters(); + if (!parser.writeLineF("{} {} {} {} {} {}\n", AngleFunctions::forms().keyword(angle.interactionForm()).substr(0, 4), angle.indexI() + 1, angle.indexJ() + 1, angle.indexK() + 1, params[0], params[1])) return false; } } - if (!parser.writeLineF("finish")) + if (!parser.writeLineF("finish\n")) return false; - + + int vdw = 0; + + for (const auto &species : speciesPopulations){ + for (const auto &atom : species.first->atoms()){ + const auto params = atom.atomType()->interactionPotential().parameters(); + if (params[0] != 0 || params[1] != 0) + vdw++; + } + } + + if (!parser.writeLineF("vdw {}\n", vdw)) + return false; + + for (const auto &species : speciesPopulations){ + for (const auto &atom : species.first->atoms()){ + const auto params = atom.atomType()->interactionPotential().parameters(); + if (params[0] != 0 || params[1] != 0) + if (!parser.writeLineF("{} {} {} {} {}\n", atom.atomType()->name(), atom.atomType()->name(), "LJ", params[0], params[1])) + return false; + } + } + if (!parser.writeLineF("close")) + return false; return true; } diff --git a/src/io/export/dlPolyField.h b/src/io/export/dlPolyField.h index 0c9adb2ddc..ac9a24b3e1 100644 --- a/src/io/export/dlPolyField.h +++ b/src/io/export/dlPolyField.h @@ -9,11 +9,11 @@ class Configuration; // Coordinate Export Formats -class DlPolyControlExportFileFormat : public FileAndFormat +class DlPolyFieldExportFileFormat : public FileAndFormat { public: - // Available DlPolyControl formats - enum class DlPolyControlExportFormat + // Available DlPolyField formats + enum class DlPolyFieldExportFormat { DLPOLY }; From 3288e7d6635060224f79f66372a8f8718e7510ef Mon Sep 17 00:00:00 2001 From: jounaidr Date: Thu, 14 Dec 2023 19:31:01 +0000 Subject: [PATCH 11/21] cleanup imports --- src/io/export/dlPolyControl.cpp | 4 ---- src/io/export/dlPolyField.cpp | 41 --------------------------------- 2 files changed, 45 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index eee2beeb87..c389d44dcc 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -4,11 +4,7 @@ #include "io/export/dlPolyControl.h" #include "base/lineParser.h" #include "base/sysFunc.h" -#include "classes/atomType.h" -#include "classes/box.h" #include "classes/configuration.h" -#include "classes/speciesAtom.h" -#include "data/atomicMasses.h" DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view filename, DlPolyControlExportFormat format) : FileAndFormat(formats_, filename, (int)format) diff --git a/src/io/export/dlPolyField.cpp b/src/io/export/dlPolyField.cpp index 3937f7cd0f..c06c754ade 100644 --- a/src/io/export/dlPolyField.cpp +++ b/src/io/export/dlPolyField.cpp @@ -4,51 +4,10 @@ #include "io/export/dlPolyField.h" #include "base/lineParser.h" #include "base/sysFunc.h" -#include "classes/atomType.h" -#include "classes/box.h" #include "classes/configuration.h" -#include "classes/speciesAtom.h" #include "data/atomicMasses.h" -#include "classes/atom.h" -#include "classes/atomType.h" -#include "classes/atomTypeMix.h" -#include "classes/molecule.h" -#include -#include - -#include "base/serialiser.h" -#include "base/version.h" -#include "classes/atom.h" -#include "classes/atomTypeMix.h" -#include "classes/box.h" -#include "classes/cellArray.h" -#include "classes/molecule.h" -#include "classes/siteStack.h" -#include "io/import/coordinates.h" -#include "items/list.h" -#include "kernels/potentials/base.h" -#include "math/data1D.h" -#include "math/histogram1D.h" -#include "math/interpolator.h" -#include "module/layer.h" -#include "procedure/procedure.h" -#include "templates/vector3.h" -#include -#include -#include -#include - -#include "io/export/coordinates.h" -#include "base/lineParser.h" -#include "base/sysFunc.h" -#include "classes/atomType.h" -#include "classes/box.h" -#include "classes/configuration.h" -#include "classes/speciesAtom.h" -#include "data/atomicMasses.h" - DlPolyFieldExportFileFormat::DlPolyFieldExportFileFormat(std::string_view filename, DlPolyFieldExportFormat format) : FileAndFormat(formats_, filename, (int)format) { From 895c36c4031d65e1fc0ec982ae771a9d8c83a309 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Tue, 19 Dec 2023 20:48:22 +0000 Subject: [PATCH 12/21] vdw_mix_method and ensemble_method implemented in CONTROL file --- src/io/export/dlPolyControl.cpp | 31 +++++++++++++++++++++++++++---- src/io/export/dlPolyControl.h | 4 ++-- src/modules/dlPoly/dlPoly.cpp | 6 +++++- src/modules/dlPoly/dlPoly.h | 21 +++++++++++++++++++++ src/modules/dlPoly/process.cpp | 10 ++++++++++ 5 files changed, 65 insertions(+), 7 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index c389d44dcc..7ec7c1ee73 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -15,7 +15,7 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi // Export DlPolyControl as CONTROL -bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) +bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) { // Export title if (!parser.writeLineF("title {} @ {}\n\n", cfg->name(), cfg->contentsVersion())) @@ -40,7 +40,7 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("ensemble {}\n", ensemble)) return false; - if (!parser.writeLineF("ensemble_method hoover\n", cutoffDistance.value())) + if (!parser.writeLineF("ensemble_method {}\n", ensembleMethod)) return false; if (!parser.writeLineF("ensemble_thermostat_coupling {} ps\n", ensembleThermostatCoupling)) return false; @@ -68,14 +68,35 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("coul_precision {}\n", coulPrecision)) return false; - if (!parser.writeLineF("vdw_mix_method Lorentz-Berthelot\n")) + + // Note interaction is taken from the first atom only, not sure how to do this ... + const auto sri = cfg->speciesPopulations()[0].first->atoms()[0].atomType()->interactionPotential().form(); + std::string vdw_mix_method; + + switch (sri) + { + case (ShortRangeFunctions::Form::None): + vdw_mix_method = "Off"; + break; + case (ShortRangeFunctions::Form::LennardJones): + vdw_mix_method = "Lorentz-Berthelot"; + break; + case (ShortRangeFunctions::Form::LennardJonesGeometric): + vdw_mix_method = "Hogervorst"; + break; + default: + throw(std::runtime_error(fmt::format("Short-range type {} is not accounted for in PairPotential::setUp().\n", + ShortRangeFunctions::forms().keyword(sri)))); + } + + if (!parser.writeLineF("vdw_mix_method {}\n", vdw_mix_method)) return false; return true; } // Export DlPolyControl using current filename and format -bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) +bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) { // Open the file LineParser parser; @@ -95,7 +116,9 @@ bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForce capForces, capForcesAt, cutoffDistance, + padding, ensemble, + ensembleMethod, ensembleThermostatCoupling, timestepVariable, fixedTimestep, diff --git a/src/io/export/dlPolyControl.h b/src/io/export/dlPolyControl.h index 243b779494..6e8b736f75 100644 --- a/src/io/export/dlPolyControl.h +++ b/src/io/export/dlPolyControl.h @@ -39,11 +39,11 @@ class DlPolyControlExportFileFormat : public FileAndFormat */ private: // Export Configuration as DL_POLY CONTROL - bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); + bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); public: // Export Configuration using current filename and format - bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, std::string ensemble, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); + bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); }; diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 8a3a9da5d7..1797772eb3 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -19,6 +19,7 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.setOrganisation("Options", "File"); keywords_.add("CONTROL", "File / format for CONTROL", dlPolyControlFormat_, "EndFormat"); keywords_.add("FIELD", "File / format for FIELD", dlPolyFieldFormat_, "EndFormat"); + //keywords_.add("CONFIG", "File / format for CONFIG", coordinatesFormat_, "EndFormat"); keywords_.setOrganisation("Options", "Simulation"); keywords_.add("NSteps", "Number of DlPoly steps to perform", nSteps_, 1); @@ -48,9 +49,12 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.add("CapForcesAt", "Set cap on allowable force (kJ/mol) per atom", capForcesAt_, 0.0); keywords_.add( "CutoffDistance", "Interatomic cutoff distance to use for energy calculation (0.0 to use pair potential range)", - cutoffDistance_, 0.0, std::nullopt, 0.1, "Use PairPotential Range"); + cutoffDistance_, 0.0, std::nullopt, 1.0, "Use PairPotential Range"); + keywords_.add("padding", "padding", padding_, 0.0); keywords_.add>( "Ensemble", " Set ensemble constraints", ensemble_, DlPolyModule::ensemble()); + keywords_.add>( + "EnsembleMethod", " Set ensemble method", ensembleMethod_, DlPolyModule::ensembleMethod()); keywords_.add("EnsembleThermostatCoupling", "Set thermostat relaxation/decorrelation times (use ensemble_thermostat_friction for langevin)", ensembleThermostatCoupling_, 0.0); } \ No newline at end of file diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 754e64ddc9..a9b7b55512 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -57,6 +57,23 @@ class DlPolyModule : public Module }; // Return EnumOptions for Ensemble static EnumOptions ensemble(); + + // Ensemble Enum + enum EnsembleMethod + { + Evans, + Langevin, + Andersen, + Berendsen, + Hoover, + gentle, + ttm, + dpds1, + dpds2, + MTK, + }; + // Return EnumOptions for Ensemble method + static EnumOptions ensembleMethod(); private: // Filename and format for CONTROL export @@ -71,8 +88,12 @@ class DlPolyModule : public Module double capForcesAt_{1.0e7}; // Interatomic cutoff distance to employ std::optional cutoffDistance_; + // Padding + double padding_{0.0}; // Set Ensemble DlPolyModule::Ensemble ensemble_{Ensemble::NVE}; + // Set Ensemble + DlPolyModule::EnsembleMethod ensembleMethod_{EnsembleMethod::Evans}; // Set Ensemble Thermostat Coupling double ensembleThermostatCoupling_{0.0}; // Timestep type to employ diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index d6b8ed192c..53a00eae12 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -35,6 +35,14 @@ EnumOptions DlPolyModule::ensemble() {{Ensemble::NVE, "NVE"}, {Ensemble::PMF, "PMF"}, {Ensemble::NVT, "NVT"}, {Ensemble::NPT, "NPT"}, {Ensemble::NST, "NST"}}); } +// Return EnumOptions for Ensemble +EnumOptions DlPolyModule::ensembleMethod() +{ + return EnumOptions( + "ensembleMethod", + {{EnsembleMethod::Evans, "Evans"}, {EnsembleMethod::Langevin, "Langevin"}, {EnsembleMethod::Andersen, "Andersen"}, {EnsembleMethod::Berendsen, "Berendsen"}, {EnsembleMethod::Hoover, "Hoover"}, {EnsembleMethod::gentle, "gentle"}, {EnsembleMethod::ttm, "ttm"}, {EnsembleMethod::dpds1, "dpds1"}, {EnsembleMethod::dpds2, "dpds2"}, {EnsembleMethod::MTK, "MTK"}}); +} + // Run main processing Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) { @@ -55,7 +63,9 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) capForces_, capForcesAt_, cutoffDistance_, + padding_, ensemble().keyword(ensemble_), + ensembleMethod().keyword(ensembleMethod_), ensembleThermostatCoupling_, timestepVariable_, fixedTimestep_, From c2b5f31be09384153999f34b8764ec9d4eb29155 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 20 Dec 2023 19:10:15 +0000 Subject: [PATCH 13/21] minimisation options implemented --- src/io/export/dlPolyControl.cpp | 21 ++++++++++++++++----- src/io/export/dlPolyControl.h | 4 ++-- src/modules/dlPoly/dlPoly.cpp | 5 ++++- src/modules/dlPoly/dlPoly.h | 16 +++++++++++++++- src/modules/dlPoly/process.cpp | 15 +++++++++++++-- 5 files changed, 50 insertions(+), 11 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 7ec7c1ee73..2484cfef23 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -15,7 +15,7 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi // Export DlPolyControl as CONTROL -bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) +bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision, std::string minimisationCriterion, double minimisationTolerance, double minimisationFrequency) { // Export title if (!parser.writeLineF("title {} @ {}\n\n", cfg->name(), cfg->contentsVersion())) @@ -36,10 +36,13 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("stats_frequency {} steps\n", outputFrequency.value())) return false; - if (!parser.writeLineF("cutoff {} ang\n", cutoffDistance.value())) + if (!parser.writeLineF("cutoff {} nm\n", cutoffDistance.value())) + return false; + if (!parser.writeLineF("padding {} ang\n", cutoffDistance.value())) return false; if (!parser.writeLineF("ensemble {}\n", ensemble)) return false; + //anything other than NVE method and thermsostat cooling if (!parser.writeLineF("ensemble_method {}\n", ensembleMethod)) return false; if (!parser.writeLineF("ensemble_thermostat_coupling {} ps\n", ensembleThermostatCoupling)) @@ -88,15 +91,20 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati throw(std::runtime_error(fmt::format("Short-range type {} is not accounted for in PairPotential::setUp().\n", ShortRangeFunctions::forms().keyword(sri)))); } - if (!parser.writeLineF("vdw_mix_method {}\n", vdw_mix_method)) return false; + if (!parser.writeLineF("minimisation_criterion {}\n", minimisationCriterion)) + return false; + if (!parser.writeLineF("minimisation_tolerance {}\n", minimisationTolerance)) + return false; + if (!parser.writeLineF("minimisation_frequency {}\n", minimisationFrequency)) + return false; return true; } // Export DlPolyControl using current filename and format -bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision) +bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision, std::string minimisationCriterion, double minimisationTolerance, double minimisationFrequency) { // Open the file LineParser parser; @@ -129,7 +137,10 @@ bool DlPolyControlExportFileFormat::exportData(Configuration *cfg, bool capForce trajectoryFrequency, trajectoryKey, coulMethod, - coulPrecision); + coulPrecision, + minimisationCriterion, + minimisationTolerance, + minimisationFrequency); break; default: throw(std::runtime_error(fmt::format("DlPolyControl format '{}' export has not been implemented.\n", diff --git a/src/io/export/dlPolyControl.h b/src/io/export/dlPolyControl.h index 6e8b736f75..c86bbdc010 100644 --- a/src/io/export/dlPolyControl.h +++ b/src/io/export/dlPolyControl.h @@ -39,11 +39,11 @@ class DlPolyControlExportFileFormat : public FileAndFormat */ private: // Export Configuration as DL_POLY CONTROL - bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); + bool exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision, std::string minimisationCriterion, double minimisationTolerance, double minimisationFrequency); public: // Export Configuration using current filename and format - bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision); + bool exportData(Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision, std::string minimisationCriterion, double minimisationTolerance, double minimisationFrequency); }; diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 1797772eb3..088ed463c8 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -56,5 +56,8 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.add>( "EnsembleMethod", " Set ensemble method", ensembleMethod_, DlPolyModule::ensembleMethod()); keywords_.add("EnsembleThermostatCoupling", "Set thermostat relaxation/decorrelation times (use ensemble_thermostat_friction for langevin)", ensembleThermostatCoupling_, 0.0); - + keywords_.add>( + "MinimisationCriterion", " Set minimisation criterion", minimisationCriterion_, DlPolyModule::minimisationCriterion()); + keywords_.add("minimisationTolerance", "Minimisation Tolerance", minimisationTolerance_, 0.0); + keywords_.add("minimisationFrequency", "Minimisation Frequency", minimisationFrequency_, 0.0); } \ No newline at end of file diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index a9b7b55512..f8b791d1e9 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -74,6 +74,15 @@ class DlPolyModule : public Module }; // Return EnumOptions for Ensemble method static EnumOptions ensembleMethod(); + + enum MinimisationCriterion + { + Off, + Force, + Energy, + Distance + }; + static EnumOptions minimisationCriterion(); private: // Filename and format for CONTROL export @@ -116,7 +125,12 @@ class DlPolyModule : public Module DlPolyModule::CoulMethod coulMethod_{CoulMethod::off}; // Calculate electrostatics with given precision double coulPrecision_{0.0}; - + // Set MinimisationCriterion + DlPolyModule::MinimisationCriterion minimisationCriterion_{MinimisationCriterion::Off}; + // Minimisation Tolerance + double minimisationTolerance_{0.0}; + // Minimsation Frequency + double minimisationFrequency_{0.0}; /* * Functions diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 53a00eae12..959496c4f7 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -35,7 +35,7 @@ EnumOptions DlPolyModule::ensemble() {{Ensemble::NVE, "NVE"}, {Ensemble::PMF, "PMF"}, {Ensemble::NVT, "NVT"}, {Ensemble::NPT, "NPT"}, {Ensemble::NST, "NST"}}); } -// Return EnumOptions for Ensemble +// Return EnumOptions for EnsembleMethod EnumOptions DlPolyModule::ensembleMethod() { return EnumOptions( @@ -43,6 +43,14 @@ EnumOptions DlPolyModule::ensembleMethod() {{EnsembleMethod::Evans, "Evans"}, {EnsembleMethod::Langevin, "Langevin"}, {EnsembleMethod::Andersen, "Andersen"}, {EnsembleMethod::Berendsen, "Berendsen"}, {EnsembleMethod::Hoover, "Hoover"}, {EnsembleMethod::gentle, "gentle"}, {EnsembleMethod::ttm, "ttm"}, {EnsembleMethod::dpds1, "dpds1"}, {EnsembleMethod::dpds2, "dpds2"}, {EnsembleMethod::MTK, "MTK"}}); } +// Return EnumOptions for MinimisationCriterion +EnumOptions DlPolyModule::minimisationCriterion() +{ + return EnumOptions( + "minimisationCriterion", + {{MinimisationCriterion::Off, "Off"}, {MinimisationCriterion::Force, "Force"}, {MinimisationCriterion::Energy, "Energy"}, {MinimisationCriterion::Distance, "Distance"}}); +} + // Run main processing Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) { @@ -76,7 +84,10 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) trajectoryFrequency_, trajectoryKey().keyword(trajectoryKey_), coulMethod().keyword(coulMethod_), - coulPrecision_)) + coulPrecision_, + minimisationCriterion().keyword(minimisationCriterion_), + minimisationTolerance_, + minimisationFrequency_)) { Messenger::print("Export: Failed to export CONTROL file '{}'.\n", dlPolyControlFormat_.filename()); moduleContext.processPool().decideFalse(); From 93e927041064793e1c71ebe0d41cbbd08b94689a Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 20 Dec 2023 19:18:07 +0000 Subject: [PATCH 14/21] conditional esmemble_method/emsemble_thermostat_coupling --- src/io/export/dlPolyControl.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 2484cfef23..91d21f442d 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -42,11 +42,13 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("ensemble {}\n", ensemble)) return false; - //anything other than NVE method and thermsostat cooling - if (!parser.writeLineF("ensemble_method {}\n", ensembleMethod)) - return false; - if (!parser.writeLineF("ensemble_thermostat_coupling {} ps\n", ensembleThermostatCoupling)) - return false; + //anything other than NVE method and thermsostat cooling + if (ensemble != "NVE"){ + if (!parser.writeLineF("ensemble_method {}\n", ensembleMethod)) + return false; + if (!parser.writeLineF("ensemble_thermostat_coupling {} ps\n", ensembleThermostatCoupling)) + return false; + } if (capForces && !parser.writeLineF("equilibration_force_cap {} k_b.temp/ang\n", capForcesAt)) return false; if (!parser.writeLineF("time_run {} steps\n", nSteps)) From 297754c3171c47b00d68f4612a6f65c7fcb9e56d Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 20 Dec 2023 19:31:18 +0000 Subject: [PATCH 15/21] add config file to dlpoly module --- src/modules/dlPoly/dlPoly.cpp | 2 +- src/modules/dlPoly/dlPoly.h | 3 +++ src/modules/dlPoly/process.cpp | 14 ++++++++++++++ 3 files changed, 18 insertions(+), 1 deletion(-) diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 088ed463c8..2f5dbe3b10 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -19,7 +19,7 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) keywords_.setOrganisation("Options", "File"); keywords_.add("CONTROL", "File / format for CONTROL", dlPolyControlFormat_, "EndFormat"); keywords_.add("FIELD", "File / format for FIELD", dlPolyFieldFormat_, "EndFormat"); - //keywords_.add("CONFIG", "File / format for CONFIG", coordinatesFormat_, "EndFormat"); + keywords_.add("CONFIG", "File / format for CONFIG", coordinatesFormat_, "EndFormat"); keywords_.setOrganisation("Options", "Simulation"); keywords_.add("NSteps", "Number of DlPoly steps to perform", nSteps_, 1); diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index f8b791d1e9..dfc24403e7 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -5,6 +5,7 @@ #include "io/export/dlPolyControl.h" #include "io/export/dlPolyField.h" +#include "io/export/coordinates.h" #include "base/enumOptions.h" #include "module/module.h" @@ -89,6 +90,8 @@ class DlPolyModule : public Module DlPolyControlExportFileFormat dlPolyControlFormat_; // Filename and format for FIELD export DlPolyFieldExportFileFormat dlPolyFieldFormat_; + // Filename and format for CONFIG export + CoordinateExportFileFormat coordinatesFormat_; // Target configurations Configuration *targetConfiguration_{nullptr}; // Control whether atomic forces are capped every step diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 959496c4f7..17bb5e7c33 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -100,6 +100,20 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) moduleContext.processPool().decideFalse(); return ExecutionResult::Failed; } + + if (!coordinatesFormat_.hasFilename()) + { + Messenger::error("No valid file/format set for CONFIG export.\n"); + return ExecutionResult::Failed; + } + Messenger::print("Export: Writing CONFIG file ({}) for Configuration '{}'...\n", coordinatesFormat_.formatDescription(), targetConfiguration_->name()); + + if (!coordinatesFormat_.exportData(targetConfiguration_)) + { + Messenger::print("Export: Failed to export CONFIG file '{}'.\n", coordinatesFormat_.filename()); + moduleContext.processPool().decideFalse(); + return ExecutionResult::Failed; + } moduleContext.processPool().decideTrue(); } From 585420b46d30f416d48a572e3ed55f3fc1945cfc Mon Sep 17 00:00:00 2001 From: jounaidr Date: Wed, 20 Dec 2023 19:45:08 +0000 Subject: [PATCH 16/21] add ./DLPOLY.Z execution to dlpoly module --- src/modules/dlPoly/process.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 17bb5e7c33..1799b6e192 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -51,6 +51,19 @@ EnumOptions DlPolyModule::minimisationCrite {{MinimisationCriterion::Off, "Off"}, {MinimisationCriterion::Force, "Force"}, {MinimisationCriterion::Energy, "Energy"}, {MinimisationCriterion::Distance, "Distance"}}); } +std::string exec(const char* cmd) { + std::array buffer; + std::string result; + std::unique_ptr pipe(popen(cmd, "r"), pclose); + if (!pipe) { + throw std::runtime_error("popen() failed!"); + } + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + result += buffer.data(); + } + return result; +} + // Run main processing Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) { @@ -114,6 +127,8 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) moduleContext.processPool().decideFalse(); return ExecutionResult::Failed; } + + auto run = exec("./DLPOLY.Z"); moduleContext.processPool().decideTrue(); } From 4d8de145c1172db13a35238f11e04032464de479 Mon Sep 17 00:00:00 2001 From: jounaidr Date: Tue, 2 Jan 2024 20:33:55 +0000 Subject: [PATCH 17/21] cleanup CONTROL option tool tips, attempt at hardcoding filenames (not working) --- src/io/export/dlPolyControl.cpp | 5 ++--- src/modules/dlPoly/dlPoly.cpp | 38 ++++++++++++++------------------- src/modules/dlPoly/dlPoly.h | 6 +++--- 3 files changed, 21 insertions(+), 28 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 91d21f442d..d66304d170 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -38,11 +38,10 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("cutoff {} nm\n", cutoffDistance.value())) return false; - if (!parser.writeLineF("padding {} ang\n", cutoffDistance.value())) + if (!parser.writeLineF("padding {} ang\n", padding.value())) return false; if (!parser.writeLineF("ensemble {}\n", ensemble)) return false; - //anything other than NVE method and thermsostat cooling if (ensemble != "NVE"){ if (!parser.writeLineF("ensemble_method {}\n", ensembleMethod)) return false; @@ -99,7 +98,7 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("minimisation_tolerance {}\n", minimisationTolerance)) return false; - if (!parser.writeLineF("minimisation_frequency {}\n", minimisationFrequency)) + if (!parser.writeLineF("minimisation_frequency {} steps\n", minimisationFrequency)) return false; return true; diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 2f5dbe3b10..13232e5b86 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -15,49 +15,43 @@ DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) { keywords_.addTarget("Configuration", "Set target configuration for the module", targetConfiguration_); - - keywords_.setOrganisation("Options", "File"); - keywords_.add("CONTROL", "File / format for CONTROL", dlPolyControlFormat_, "EndFormat"); - keywords_.add("FIELD", "File / format for FIELD", dlPolyFieldFormat_, "EndFormat"); - keywords_.add("CONFIG", "File / format for CONFIG", coordinatesFormat_, "EndFormat"); keywords_.setOrganisation("Options", "Simulation"); - keywords_.add("NSteps", "Number of DlPoly steps to perform", nSteps_, 1); + keywords_.add("NSteps", "Set duration of simulation (inc. equilibration)", nSteps_, 1); keywords_.add("VariableTimestep", - "Whether a variable timestep should be used, determined from the maximal force vector", + "Enable variable timestep", timestepVariable_); - keywords_.add("DeltaT", "Fixed timestep (ps) to use in DlPoly simulation", fixedTimestep_, 0.0); + keywords_.add("DeltaT", "Set calculation timestep or initial timestep for variable timestep calculations", fixedTimestep_, 0.0); keywords_.add("RandomVelocities", "Whether random velocities should always be assigned before beginning DlPoly simulation", randomVelocities_); keywords_.setOrganisation("Options", "Output"); - keywords_.add("EnergyFrequency", "Frequency at which to calculate total system energy", + keywords_.add("EnergyFrequency", "Set frequency of printing results to output", energyFrequency_, 0, std::nullopt, 5, "Off"); - keywords_.add("OutputFrequency", "Frequency at which to output step information", outputFrequency_, + keywords_.add("OutputFrequency", "Set frequency of stats sampling to statis file", outputFrequency_, 0, std::nullopt, 5, "Off"); - keywords_.add("TrajectoryFrequency", "Write frequency for trajectory file", trajectoryFrequency_, 0, - std::nullopt, 5, "Off"); + keywords_.add("TrajectoryFrequency", "Interval between dumping trajectory configurations", trajectoryFrequency_, 0, std::nullopt, 5, "Off"); keywords_.add>( - "TrajectoryKey", "Set trajectory key", trajectoryKey_, DlPolyModule::trajectoryKey()); + "TrajectoryKey", "Set trajectory output, options: pos, pos-vel, pos-vel-force, compressed", trajectoryKey_, DlPolyModule::trajectoryKey()); keywords_.add>( - "CoulMethod", "Set Coul Method", coulMethod_, DlPolyModule::coulMethod()); + "CoulMethod", "Set method for electrostatics method, options: off, spme, dddp, pairwise, reaction_field, force_shifted", coulMethod_, DlPolyModule::coulMethod()); keywords_.add("Coul Precision", "Calculate electrostatics using Fennell damping (Ewald-like) with given precision", coulPrecision_, 0.0); keywords_.setOrganisation("Advanced"); keywords_.add("CapForces", "Control whether atomic forces are capped every step", capForces_); - keywords_.add("CapForcesAt", "Set cap on allowable force (kJ/mol) per atom", capForcesAt_, 0.0); + keywords_.add("CapForcesAt", "Set force cap clamping maximum force during equilibration", capForcesAt_, 0.0); keywords_.add( - "CutoffDistance", "Interatomic cutoff distance to use for energy calculation (0.0 to use pair potential range)", + "CutoffDistance", "Set the global cutoff for real-speace potentials", cutoffDistance_, 0.0, std::nullopt, 1.0, "Use PairPotential Range"); - keywords_.add("padding", "padding", padding_, 0.0); + keywords_.add("padding", "Set padding for sizing of Verlet neighbour lists", padding_, 0.0); keywords_.add>( - "Ensemble", " Set ensemble constraints", ensemble_, DlPolyModule::ensemble()); + "Ensemble", "Set ensemble constraints, options: NVE, PMF, NVT, NPT, NST", ensemble_, DlPolyModule::ensemble()); keywords_.add>( - "EnsembleMethod", " Set ensemble method", ensembleMethod_, DlPolyModule::ensembleMethod()); + "EnsembleMethod", "Set ensemble method, options NVT: Evans, Langevin, Andersen, Berendsen, Hoover, gentle, ttm, dpds1, dpds2. NP|ST: Langevin, Berendsen, Hoover, MTK.", ensembleMethod_, DlPolyModule::ensembleMethod()); keywords_.add("EnsembleThermostatCoupling", "Set thermostat relaxation/decorrelation times (use ensemble_thermostat_friction for langevin)", ensembleThermostatCoupling_, 0.0); keywords_.add>( - "MinimisationCriterion", " Set minimisation criterion", minimisationCriterion_, DlPolyModule::minimisationCriterion()); - keywords_.add("minimisationTolerance", "Minimisation Tolerance", minimisationTolerance_, 0.0); - keywords_.add("minimisationFrequency", "Minimisation Frequency", minimisationFrequency_, 0.0); + "MinimisationCriterion", "Set minimisation criterion, options: off, force, energy, distance", minimisationCriterion_, DlPolyModule::minimisationCriterion()); + keywords_.add("minimisationTolerance", "Set minimisation tolerance, units: determined by criterion", minimisationTolerance_, 0.0); + keywords_.add("minimisationFrequency", "Set minimisation frequency", minimisationFrequency_, 0.0); } \ No newline at end of file diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index dfc24403e7..6878556b07 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -87,11 +87,11 @@ class DlPolyModule : public Module private: // Filename and format for CONTROL export - DlPolyControlExportFileFormat dlPolyControlFormat_; + DlPolyControlExportFileFormat(("CONTROL"), (DlPolyControlExportFormat::DLPOLY)) dlPolyControlFormat_; // Filename and format for FIELD export - DlPolyFieldExportFileFormat dlPolyFieldFormat_; + DlPolyFieldExportFileFormat(("FIELD"), (DlPolyFieldExportFormat::DLPOLY)) dlPolyFieldFormat_; // Filename and format for CONFIG export - CoordinateExportFileFormat coordinatesFormat_; + CoordinateExportFileFormat(("CONFIG"), (CoordinateExportFormat::DLPOLY)) coordinatesFormat_; // Target configurations Configuration *targetConfiguration_{nullptr}; // Control whether atomic forces are capped every step From 96117c1f90ba43a2ffce0f4d56fafb3ccb6a718e Mon Sep 17 00:00:00 2001 From: jounaidr Date: Fri, 5 Jan 2024 18:28:48 +0000 Subject: [PATCH 18/21] comments --- src/io/export/dlPolyControl.cpp | 17 +++++++++-------- src/io/export/dlPolyField.cpp | 11 +++++++---- src/modules/dlPoly/dlPoly.h | 22 ++++++++++++---------- src/modules/dlPoly/process.cpp | 16 +++++++--------- 4 files changed, 35 insertions(+), 31 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index d66304d170..3ef36a393e 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -14,10 +14,9 @@ DlPolyControlExportFileFormat::DlPolyControlExportFileFormat(std::string_view fi } -// Export DlPolyControl as CONTROL +// Export DlPolyControl as CONTROL file bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg, bool capForces, double capForcesAt, std::optional cutoffDistance, double padding, std::string ensemble, std::string ensembleMethod, double ensembleThermostatCoupling, bool timestepVariable, double fixedTimestep, std::optional energyFrequency, int nSteps, std::optional outputFrequency, bool randomVelocities, std::optional trajectoryFrequency, std::string trajectoryKey, std::string coulMethod, double coulPrecision, std::string minimisationCriterion, double minimisationTolerance, double minimisationFrequency) { - // Export title if (!parser.writeLineF("title {} @ {}\n\n", cfg->name(), cfg->contentsVersion())) return false; if (!parser.writeLineF("io_file_config CONFIG\n")) @@ -42,23 +41,24 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("ensemble {}\n", ensemble)) return false; + // Export ensemble_method and ensemble_thermostat_coupling for ensemble types other than NVE if (ensemble != "NVE"){ if (!parser.writeLineF("ensemble_method {}\n", ensembleMethod)) return false; if (!parser.writeLineF("ensemble_thermostat_coupling {} ps\n", ensembleThermostatCoupling)) return false; - } + } + // Export equilibration_force_cap if capForces is true if (capForces && !parser.writeLineF("equilibration_force_cap {} k_b.temp/ang\n", capForcesAt)) return false; if (!parser.writeLineF("time_run {} steps\n", nSteps)) return false; - if (timestepVariable) - { - if (!parser.writeLineF("timestep_variable ON\n")) - return false; - } + // Export timestep_variable as ON if timestepVariable is true + if (timestepVariable && !parser.writeLineF("timestep_variable ON\n")) + return false; if (!parser.writeLineF("timestep {} internal_t\n", fixedTimestep)) return false; + // Export traj_calculate, traj_interval, traj_key if trajectoryFrequency is greater than 0 if (trajectoryFrequency.value_or(0) > 0) { if (!parser.writeLineF("traj_calculate ON\n")) @@ -77,6 +77,7 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati const auto sri = cfg->speciesPopulations()[0].first->atoms()[0].atomType()->interactionPotential().form(); std::string vdw_mix_method; + // Export vdw_mix_method based on short range function type switch (sri) { case (ShortRangeFunctions::Form::None): diff --git a/src/io/export/dlPolyField.cpp b/src/io/export/dlPolyField.cpp index c06c754ade..acd6d22c7f 100644 --- a/src/io/export/dlPolyField.cpp +++ b/src/io/export/dlPolyField.cpp @@ -19,7 +19,6 @@ DlPolyFieldExportFileFormat::DlPolyFieldExportFileFormat(std::string_view filena // Export DlPolyField as Field bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration *cfg) { - // Export title if (!parser.writeLineF("title {} @ {}\n", cfg->name(), cfg->contentsVersion())) return false; if (!parser.writeLineF("units Kj\n")) @@ -27,9 +26,11 @@ bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration const std::vector> speciesPopulations = cfg->speciesPopulations(); + // Export number of moleculer types if (!parser.writeLineF("moleculer types {}\n", speciesPopulations.size())) return false; + // Export number and name of molecules, and number of atoms, bonds and angles with respective properties for (const auto &species : speciesPopulations){ if (!parser.writeLineF("{}\nnummols {}\natoms {}\n", species.first->name(), species.second, species.first->nAtoms())){ return false; @@ -58,7 +59,8 @@ bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration return false; int vdw = 0; - + + // Calculate vdw based on interaction potentials for each atom for (const auto &species : speciesPopulations){ for (const auto &atom : species.first->atoms()){ const auto params = atom.atomType()->interactionPotential().parameters(); @@ -69,7 +71,8 @@ bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration if (!parser.writeLineF("vdw {}\n", vdw)) return false; - + + // Export interaction potentials and respective properties for (const auto &species : speciesPopulations){ for (const auto &atom : species.first->atoms()){ const auto params = atom.atomType()->interactionPotential().parameters(); @@ -78,7 +81,7 @@ bool DlPolyFieldExportFileFormat::exportDLPOLY(LineParser &parser, Configuration return false; } } - if (!parser.writeLineF("close")) + if (!parser.writeLineF("close\n")) return false; return true; diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 6878556b07..9fd05fe3a0 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -59,7 +59,7 @@ class DlPolyModule : public Module // Return EnumOptions for Ensemble static EnumOptions ensemble(); - // Ensemble Enum + // EnsembleMethod Enum enum EnsembleMethod { Evans, @@ -73,9 +73,10 @@ class DlPolyModule : public Module dpds2, MTK, }; - // Return EnumOptions for Ensemble method + // Return EnumOptions for EnsembleMethod static EnumOptions ensembleMethod(); + // MinimisationCriterion Enum enum MinimisationCriterion { Off, @@ -83,6 +84,7 @@ class DlPolyModule : public Module Energy, Distance }; + // Return EnumOptions for MinimisationCriterion static EnumOptions minimisationCriterion(); private: @@ -100,23 +102,23 @@ class DlPolyModule : public Module double capForcesAt_{1.0e7}; // Interatomic cutoff distance to employ std::optional cutoffDistance_; - // Padding + // Set Padding double padding_{0.0}; // Set Ensemble DlPolyModule::Ensemble ensemble_{Ensemble::NVE}; - // Set Ensemble + // Set Ensemble Method DlPolyModule::EnsembleMethod ensembleMethod_{EnsembleMethod::Evans}; // Set Ensemble Thermostat Coupling double ensembleThermostatCoupling_{0.0}; // Timestep type to employ bool timestepVariable_{false}; - // Fixed timestep (ps) to use in MD simulation + // Set fixed timestep (ps) to use in MD simulation double fixedTimestep_{5.0e-4}; - // Frequency at which to calculate total system energy + // Set frequency at which to calculate total system energy std::optional energyFrequency_{10}; - // Number of steps to perform + // Set number of steps to perform int nSteps_{50}; - // Frequency at which to output step information + // Set frequency at which to output step information std::optional outputFrequency_{5}; // Whether random velocities should always be assigned before beginning MD simulation bool randomVelocities_{false}; @@ -130,9 +132,9 @@ class DlPolyModule : public Module double coulPrecision_{0.0}; // Set MinimisationCriterion DlPolyModule::MinimisationCriterion minimisationCriterion_{MinimisationCriterion::Off}; - // Minimisation Tolerance + // Set Minimisation Tolerance double minimisationTolerance_{0.0}; - // Minimsation Frequency + // Set Minimsation Frequency double minimisationFrequency_{0.0}; /* diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 1799b6e192..8c950f0156 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -51,6 +51,7 @@ EnumOptions DlPolyModule::minimisationCrite {{MinimisationCriterion::Off, "Off"}, {MinimisationCriterion::Force, "Force"}, {MinimisationCriterion::Energy, "Energy"}, {MinimisationCriterion::Distance, "Distance"}}); } +// Execute command on cmd (MOVE TO FUNCTIONS FILE?) std::string exec(const char* cmd) { std::array buffer; std::string result; @@ -79,7 +80,8 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) { Messenger::print("Export: Writing TESTTEST file ({}) for Configuration '{}'...\n", dlPolyControlFormat_.formatDescription(), targetConfiguration_->name()); - + + // Save CONTROL file with respective properties if (!dlPolyControlFormat_.exportData(targetConfiguration_, capForces_, capForcesAt_, @@ -107,6 +109,7 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } + // Save FIELD file if (!dlPolyFieldFormat_.exportData(targetConfiguration_)) { Messenger::print("Export: Failed to export FIELD file '{}'.\n", dlPolyFieldFormat_.filename()); @@ -114,13 +117,7 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } - if (!coordinatesFormat_.hasFilename()) - { - Messenger::error("No valid file/format set for CONFIG export.\n"); - return ExecutionResult::Failed; - } - Messenger::print("Export: Writing CONFIG file ({}) for Configuration '{}'...\n", coordinatesFormat_.formatDescription(), targetConfiguration_->name()); - + // Save CONFIG file in DLPOLY format if (!coordinatesFormat_.exportData(targetConfiguration_)) { Messenger::print("Export: Failed to export CONFIG file '{}'.\n", coordinatesFormat_.filename()); @@ -128,7 +125,8 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } - auto run = exec("./DLPOLY.Z"); + // Run DLPOLY + auto run = exec("./DLPOLY.Z"); //CHANGE WHEN INSTALLATION PATH HAS BEEN DECIDED moduleContext.processPool().decideTrue(); } From 20ca332928f6ec5a61709aa5f4ac175216727d4c Mon Sep 17 00:00:00 2001 From: jounaidr Date: Fri, 5 Jan 2024 19:13:44 +0000 Subject: [PATCH 19/21] hardcode filenames and formats, fix parsing error --- src/io/export/dlPolyControl.cpp | 2 +- src/modules/dlPoly/dlPoly.h | 6 +++--- src/modules/dlPoly/process.cpp | 14 +++++++++++++- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/io/export/dlPolyControl.cpp b/src/io/export/dlPolyControl.cpp index 3ef36a393e..d6d2c3f9f6 100644 --- a/src/io/export/dlPolyControl.cpp +++ b/src/io/export/dlPolyControl.cpp @@ -37,7 +37,7 @@ bool DlPolyControlExportFileFormat::exportDLPOLY(LineParser &parser, Configurati return false; if (!parser.writeLineF("cutoff {} nm\n", cutoffDistance.value())) return false; - if (!parser.writeLineF("padding {} ang\n", padding.value())) + if (!parser.writeLineF("padding {} ang\n", padding)) return false; if (!parser.writeLineF("ensemble {}\n", ensemble)) return false; diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 9fd05fe3a0..19114d8eab 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -89,11 +89,11 @@ class DlPolyModule : public Module private: // Filename and format for CONTROL export - DlPolyControlExportFileFormat(("CONTROL"), (DlPolyControlExportFormat::DLPOLY)) dlPolyControlFormat_; + DlPolyControlExportFileFormat dlPolyControlFormat_; // Filename and format for FIELD export - DlPolyFieldExportFileFormat(("FIELD"), (DlPolyFieldExportFormat::DLPOLY)) dlPolyFieldFormat_; + DlPolyFieldExportFileFormat dlPolyFieldFormat_; // Filename and format for CONFIG export - CoordinateExportFileFormat(("CONFIG"), (CoordinateExportFormat::DLPOLY)) coordinatesFormat_; + CoordinateExportFileFormat coordinatesFormat_; // Target configurations Configuration *targetConfiguration_{nullptr}; // Control whether atomic forces are capped every step diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index 8c950f0156..a45f865bd0 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -81,6 +81,10 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) Messenger::print("Export: Writing TESTTEST file ({}) for Configuration '{}'...\n", dlPolyControlFormat_.formatDescription(), targetConfiguration_->name()); + // Set CONTROL filename and format + dlPolyControlFormat_.setFilename("CONTROL"); + dlPolyControlFormat_.setFormatByIndex(0); + // Save CONTROL file with respective properties if (!dlPolyControlFormat_.exportData(targetConfiguration_, capForces_, @@ -109,6 +113,10 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } + // Set FIELD filename and format + dlPolyFieldFormat_.setFilename("FIELD"); + dlPolyFieldFormat_.setFormatByIndex(0); + // Save FIELD file if (!dlPolyFieldFormat_.exportData(targetConfiguration_)) { @@ -117,7 +125,11 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) return ExecutionResult::Failed; } - // Save CONFIG file in DLPOLY format + // Set CONFIG filename and format + coordinatesFormat_.setFilename("CONFIG"); + coordinatesFormat_.setFormatByIndex(1); + + // Save CONFIG file if (!coordinatesFormat_.exportData(targetConfiguration_)) { Messenger::print("Export: Failed to export CONFIG file '{}'.\n", coordinatesFormat_.filename()); From d23e7002c978b0d42e8df7af485987a01a95c5ff Mon Sep 17 00:00:00 2001 From: jounaidr Date: Tue, 9 Jan 2024 11:59:46 +0000 Subject: [PATCH 20/21] add filepath for DLPOLY execution to options (section needs to be fixed) --- src/modules/dlPoly/dlPoly.cpp | 8 ++++---- src/modules/dlPoly/dlPoly.h | 2 ++ src/modules/dlPoly/process.cpp | 3 ++- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/modules/dlPoly/dlPoly.cpp b/src/modules/dlPoly/dlPoly.cpp index 13232e5b86..61af6e660e 100644 --- a/src/modules/dlPoly/dlPoly.cpp +++ b/src/modules/dlPoly/dlPoly.cpp @@ -11,16 +11,16 @@ #include "keywords/enumOptions.h" #include "keywords/speciesVector.h" #include "keywords/fileAndFormat.h" +#include "keywords/stdString.h" DlPolyModule::DlPolyModule() : Module(ModuleTypes::DlPoly) { keywords_.addTarget("Configuration", "Set target configuration for the module", targetConfiguration_); - + keywords_.add("DlPolyPath", "Filepath of ./DLPOLY.Z installation", dlPolyPath_); + keywords_.setOrganisation("Options", "Simulation"); keywords_.add("NSteps", "Set duration of simulation (inc. equilibration)", nSteps_, 1); - keywords_.add("VariableTimestep", - "Enable variable timestep", - timestepVariable_); + keywords_.add("VariableTimestep", "Enable variable timestep", timestepVariable_); keywords_.add("DeltaT", "Set calculation timestep or initial timestep for variable timestep calculations", fixedTimestep_, 0.0); keywords_.add("RandomVelocities", "Whether random velocities should always be assigned before beginning DlPoly simulation", diff --git a/src/modules/dlPoly/dlPoly.h b/src/modules/dlPoly/dlPoly.h index 19114d8eab..22ac282907 100644 --- a/src/modules/dlPoly/dlPoly.h +++ b/src/modules/dlPoly/dlPoly.h @@ -94,6 +94,8 @@ class DlPolyModule : public Module DlPolyFieldExportFileFormat dlPolyFieldFormat_; // Filename and format for CONFIG export CoordinateExportFileFormat coordinatesFormat_; + // Set DLPOLY.Z file path + std::string dlPolyPath_; // Target configurations Configuration *targetConfiguration_{nullptr}; // Control whether atomic forces are capped every step diff --git a/src/modules/dlPoly/process.cpp b/src/modules/dlPoly/process.cpp index a45f865bd0..05b4834212 100644 --- a/src/modules/dlPoly/process.cpp +++ b/src/modules/dlPoly/process.cpp @@ -138,7 +138,8 @@ Module::ExecutionResult DlPolyModule::process(ModuleContext &moduleContext) } // Run DLPOLY - auto run = exec("./DLPOLY.Z"); //CHANGE WHEN INSTALLATION PATH HAS BEEN DECIDED + char charArr[dlPolyPath_.length() + 1]; + auto run = exec(strcpy(charArr, dlPolyPath_.c_str())); moduleContext.processPool().decideTrue(); } From c82771592981a6a008b58719fad758a047f9474f Mon Sep 17 00:00:00 2001 From: jounaidr Date: Tue, 9 Jan 2024 16:22:56 +0000 Subject: [PATCH 21/21] initial test file added --- tests/modules/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt index 05f445e6b2..bdebcc4717 100644 --- a/tests/modules/CMakeLists.txt +++ b/tests/modules/CMakeLists.txt @@ -18,3 +18,4 @@ dissolve_add_test(SRC sdf.cpp) dissolve_add_test(SRC siteRDF.cpp) dissolve_add_test(SRC sq.cpp) dissolve_add_test(SRC xRaySQ.cpp) +dissolve_add_test(SRC dlPoly.cpp)