diff --git a/src/base/parserLibrary.cpp b/src/base/parserLibrary.cpp index fe5202978f..df7009a88b 100644 --- a/src/base/parserLibrary.cpp +++ b/src/base/parserLibrary.cpp @@ -54,8 +54,25 @@ Parser real() // A parser that accepts a 3-vector of floating point numbers Parser vector3() { - return (real() & spaces() >> real() & spaces() >> real()) + // Annoyingly, the Moscito file format does not guarantee spaces + // between the numbers when the numbers are negative. Luckily, we + // can still parse this scenario. + return (real() & maybe(spaces()) >> real() & maybe(spaces()) >> real()) .apply([](double x, double y, double z) { return Vector3(x, y, z); }); } +// A parser that accepts a 3-vector of integers +Parser vector3i() +{ + return (integer() & maybe(spaces()) >> integer() & maybe(spaces()) >> integer()) + .apply([](auto x, auto y, auto z) { return Vector3i(x, y, z); }); +} + +// A parser that accepts a 3x3 matrix of floating point numbers +Parser matrix3() +{ + return (maybe(inlineSpaces()) >> vector3() << spaces() & maybe(inlineSpaces()) >> vector3() << spaces() & + maybe(inlineSpaces()) >> vector3() << spaces()) + .apply([](const auto m1, const auto m2, const auto m3) { return Matrix3(m1, m2, m3); }); +} } // namespace Parsers diff --git a/src/base/parserLibrary.h b/src/base/parserLibrary.h index 25772cc015..7f02f662d9 100644 --- a/src/base/parserLibrary.h +++ b/src/base/parserLibrary.h @@ -4,7 +4,9 @@ #pragma once #include "base/applicative.h" +#include "math/matrix3.h" #include "math/vector3.h" +#include "math/vector3i.h" namespace Parsers { @@ -20,5 +22,10 @@ Parser real(); // A parser that accepts a 3-vector of floating point numbers Parser vector3(); +// A parser that accepts a 3-vector of integers +Parser vector3i(); + +// A parser that accepts a 3x3 matrix of floating point numbers +Parser matrix3(); }; // namespace Parsers diff --git a/src/nodes/importDLPOLYStructure.cpp b/src/nodes/importDLPOLYStructure.cpp index 1eb19bfd36..c8377ea5fe 100644 --- a/src/nodes/importDLPOLYStructure.cpp +++ b/src/nodes/importDLPOLYStructure.cpp @@ -2,6 +2,9 @@ // Copyright (c) 2026 Team Dissolve and contributors #include "nodes/importDLPOLYStructure.h" +#include "base/applicative.h" +#include "base/parserLibrary.h" +#include "data/elements.h" ImportDLPOLYStructureNode::ImportDLPOLYStructureNode(Graph *parentGraph) : Node(parentGraph) { @@ -30,50 +33,42 @@ std::string_view ImportDLPOLYStructureNode::summary() const { return "Import a D // Perform processing NodeConstants::ProcessResult ImportDLPOLYStructureNode::process() { - // Open file and check that we're OK to proceed importing from it - LineParser parser; - if ((!parser.openInput(filePath_)) || (!parser.isFileGoodForReading())) + std::ifstream infile{filePath_}; + if (!infile) return error("Couldn't open file '{}' for loading coordinates data.\n", filePath_); - - // Skip title - if (parser.skipLines(1) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - - // Import in keytrj, imcon, and number of atoms, and initialise arrays - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - - auto keytrj = parser.argi(0); - auto imcon = parser.argi(1); - auto nAtoms = parser.hasArg(2) ? parser.argi(2) : 0; - if (nAtoms == 0) + auto head = header().parse(infile); + if (!head) + return error("Failed to parse file header"); + auto &[title, keytrj, imcon, natoms] = std::get<0>(*head); + if (!natoms) message(" --> Expecting coordinates for an unknown number of atoms (DLPOLY keytrj={}, imcon={}) - will read " "until end of file.\n", - nAtoms, keytrj, imcon); + keytrj, imcon); else - message(" --> Expecting coordinates for {} atoms (DLPOLY keytrj={}, imcon={}).\n", nAtoms, keytrj, imcon); - - return read(parser, keytrj, imcon, nAtoms, structure_, forces_); + message(" --> Expecting coordinates for {} atoms (DLPOLY keytrj={}, imcon={}).\n", *natoms, keytrj, imcon); + return read(infile, keytrj, imcon, natoms.value_or(0), structure_, forces_); } // Read structure from the specified file parser -NodeConstants::ProcessResult ImportDLPOLYStructureNode::read(LineParser &parser, int keytrj, int imcon, int nAtoms, +NodeConstants::ProcessResult ImportDLPOLYStructureNode::read(std::istream &input, int keytrj, int imcon, int nAtoms, Structure &structure, OptionalReferenceWrapper> optForces) { - /* - * Import DL_POLY coordinates information through the specified line parser. - * We assume HISf, CONFIG or REVCON format (only the first two lines differ) - * - * Line 1: Title - * Line 2: keytrj imcon natoms [] - * Line 3-5: cell matrix (if imcon > 0) - * Line 6: atomtype id - * Line 7: rx ry rz - * Line 8: vx vy vz if (keytrj > 0) - * Line 9: fx fy fz if (keytrj > 1) - * ... - */ + // /* + // * Import DL_POLY coordinates information through the specified line parser. + // * We assume HISf, CONFIG or REVCON format (only the first two lines differ) + // * + // * Line 1: Title + // * Line 2: keytrj imcon natoms [] + // * Line 3-5: cell matrix (if imcon > 0) + // * Line 6: atomtype id + // * Line 7: rx ry rz + // * Line 8: vx vy vz if (keytrj > 0) + // * Line 9: fx fy fz if (keytrj > 1) + // * ... + // */ + + using namespace Parsers; // Clear storage objects structure.clear(); @@ -82,47 +77,43 @@ NodeConstants::ProcessResult ImportDLPOLYStructureNode::read(LineParser &parser, // Read cell information if given if (imcon > 0) { - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) + auto mat = matrix3().parse(input); + if (!mat) return NodeConstants::ProcessResult::Failed; - auto m1 = parser.arg3d(0); - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - auto m2 = parser.arg3d(0); - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - auto m3 = parser.arg3d(0); - structure.createBox(Matrix3(m1, m2, m3)); + structure.createBox(std::get<0>(*mat)); } - // Loop over atoms (either a specified number, or until we reach the end of the file - while (!parser.eofOrBlank()) + auto atomType = graphs() << inlineSpaces() & natural() << inlineSpaces() & + maybe(real() << maybe(inlineSpaces() << real()) << newlines()); + + if (keytrj == 0) { - // Skip atomname line - if (parser.skipLines(1) != LineParser::Success) + auto terms = atomType >> spaces() >> vector3() << spaces(); + auto result = some(terms).parse(input); + if (!result) return NodeConstants::ProcessResult::Failed; - - // Read position - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) + for (auto position : std::get<0>(*result)) + structure.addAtom(Elements::Unknown, position); + } + else if (keytrj == 1) + { + auto terms = atomType >> vector3() << spaces() & vector3() << spaces(); + auto result = some(terms).parse(input); + if (!result) return NodeConstants::ProcessResult::Failed; - structure.addAtom(Elements::Unknown, parser.arg3d(0)); - - // Read velocity if present - if (keytrj > 0) - { - if (parser.skipLines(1) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - } - - // Read forces if present - if (keytrj > 1) + for (auto &[position, velocity] : std::get<0>(*result)) + structure.addAtom(Elements::Unknown, position); + } + else + { + auto result = some(atom()).parse(input); + if (!result) + return NodeConstants::ProcessResult::Failed; + for (auto &[position, velocity, force] : std::get<0>(*result)) { - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - forces.push_back(parser.arg3d(0)); + structure.addAtom(Elements::Unknown, position); + forces.push_back(*force); } - - if ((nAtoms > 0) && (structure.nAtoms() == nAtoms)) - break; } // Copy forces out? @@ -131,3 +122,18 @@ NodeConstants::ProcessResult ImportDLPOLYStructureNode::read(LineParser &parser, return NodeConstants::ProcessResult::Success; } + +Parsers::Parser>> ImportDLPOLYStructureNode::header() +{ + using namespace Parsers; + return inlines() << newlines() & maybe(inlineSpaces()) >> natural() << inlineSpaces() & natural() & + maybe(inlineSpaces() >> natural() << maybe(inlineSpaces() << real())) << spaces(); +} + +Parsers::Parser, std::optional>> ImportDLPOLYStructureNode::atom() +{ + using namespace Parsers; + auto atomType = + graphs() << inlineSpaces() & natural() & maybe(inlineSpaces() >> real()) & maybe(inlineSpaces() >> real()) << spaces(); + return atomType >> vector3() << spaces() & maybe(vector3() << spaces()) & maybe(vector3() << spaces()); +} diff --git a/src/nodes/importDLPOLYStructure.h b/src/nodes/importDLPOLYStructure.h index cefb1d0b21..d1dda9ab9e 100644 --- a/src/nodes/importDLPOLYStructure.h +++ b/src/nodes/importDLPOLYStructure.h @@ -3,6 +3,7 @@ #pragma once +#include "base/applicative.h" #include "classes/structure.h" #include "nodes/node.h" @@ -41,6 +42,8 @@ class ImportDLPOLYStructureNode : public Node public: // Read structure from the specified file parser - static NodeConstants::ProcessResult read(LineParser &parser, int keytrj, int imcon, int nAtoms, Structure &structure, + static NodeConstants::ProcessResult read(std::istream &input, int keytrj, int imcon, int nAtoms, Structure &structure, OptionalReferenceWrapper> forces = {}); -}; \ No newline at end of file + static Parsers::Parser>> header(); + static Parsers::Parser, std::optional>> atom(); +}; diff --git a/src/nodes/importDLPOLYTrajectory.cpp b/src/nodes/importDLPOLYTrajectory.cpp index db4cf34933..3fa49f3b57 100644 --- a/src/nodes/importDLPOLYTrajectory.cpp +++ b/src/nodes/importDLPOLYTrajectory.cpp @@ -2,6 +2,8 @@ // Copyright (c) 2026 Team Dissolve and contributors #include "nodes/importDLPOLYTrajectory.h" +#include "base/applicative.h" +#include "base/parserLibrary.h" #include "nodes/importDLPOLYStructure.h" ImportDLPOLYTrajectoryNode::ImportDLPOLYTrajectoryNode(Graph *parentGraph) : Node(parentGraph) @@ -38,35 +40,41 @@ std::string_view ImportDLPOLYTrajectoryNode::summary() const // Perform processing NodeConstants::ProcessResult ImportDLPOLYTrajectoryNode::process() { + using namespace Parsers; message("Reading DL_POLY trajectory file frame from '{}'...\n", filePath_); - // Open the file - LineParser parser; - if ((!parser.openInput(filePath_)) || (!parser.isFileGoodForReading())) + std::ifstream infile{filePath_}; + if (!infile) { error("Couldn't open trajectory file '{}'.\n", filePath_); return NodeConstants::ProcessResult::Failed; } // Seek to the next file position - parser.seekg(filePosition_); + infile.seekg(filePosition_); // Read first line: 'timestep - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) + auto head = header().parse(infile); + if (!head) return NodeConstants::ProcessResult::Failed; + auto &[stepno, natoms, keytrj, imcon, _] = std::get<0>(*head); - auto keytrj = parser.argi(3); - auto imcon = parser.argi(4); - auto nAtoms = parser.hasArg(2) ? parser.argi(2) : 0; - message(" --> Expecting coordinates for {} atoms (DLPOLY keytrj={}, imcon={}).\n", nAtoms, keytrj, imcon); + message(" --> Expecting coordinates for {} atoms (DLPOLY keytrj={}, imcon={}).\n", natoms, keytrj, imcon); // Get the frame read result - auto frameResult = ImportDLPOLYStructureNode::read(parser, parser.argi(3), parser.argi(4), parser.argi(2), structure_); + auto frameResult = ImportDLPOLYStructureNode::read(infile, keytrj, imcon, natoms, structure_); if (frameResult != NodeConstants::ProcessResult::Success) return frameResult; // Store the new trajectory file position - filePosition_ = parser.tellg(); + filePosition_ = infile.tellg(); return NodeConstants::ProcessResult::Success; } + +Parsers::Parser> ImportDLPOLYTrajectoryNode::header() +{ + using namespace Parsers; + return literal("timestep") >> inlineSpaces() >> natural() << inlineSpaces() & natural() << inlineSpaces() & + natural() << inlineSpaces() & natural() << inlineSpaces() & real() << spaces(); +} diff --git a/src/nodes/importDLPOLYTrajectory.h b/src/nodes/importDLPOLYTrajectory.h index 8d1bf7317a..7cc67c935a 100644 --- a/src/nodes/importDLPOLYTrajectory.h +++ b/src/nodes/importDLPOLYTrajectory.h @@ -3,6 +3,7 @@ #pragma once +#include "base/applicative.h" #include "classes/structure.h" #include "nodes/node.h" @@ -37,4 +38,8 @@ class ImportDLPOLYTrajectoryNode : public Node protected: // Perform processing NodeConstants::ProcessResult process() override; -}; \ No newline at end of file + + public: + // Parse file header + static Parsers::Parser> header(); +}; diff --git a/src/nodes/importDLPUtilsPDens.cpp b/src/nodes/importDLPUtilsPDens.cpp index 85ca44622c..5384aaa7f8 100644 --- a/src/nodes/importDLPUtilsPDens.cpp +++ b/src/nodes/importDLPUtilsPDens.cpp @@ -2,7 +2,8 @@ // Copyright (c) 2026 Team Dissolve and contributors #include "nodes/importDLPUtilsPDens.h" -#include "base/lineParser.h" +#include "base/applicative.h" +#include "base/parserLibrary.h" ImportDLPUtilsPDensNode::ImportDLPUtilsPDensNode(Graph *parentGraph) : Node(parentGraph) { @@ -51,49 +52,39 @@ bool ImportDLPUtilsPDensNode::read(Data3D &data, std::string filePath) * Line 4 : loop order (e.g. 'zyx') * Line 5+: data (N = gridx*gridy*gridz) */ + using namespace Parsers; data.clear(); + auto firstLine = inlineSpaces() >> vector3i() & inlineSpaces() >> vector3i() & inlineSpaces() >> vector3i() << spaces(); + auto secondLine = vector3() << inlineSpaces() & vector3() << inlineSpaces() & vector3() << spaces(); + auto thirdLine = vector3() << spaces(); + auto fourthLine = literal("zyx") << spaces(); + auto header = firstLine & secondLine & thirdLine << fourthLine; + // Open file and check that we're OK to proceed importing from it - LineParser parser; - if ((!parser.openInput(filePath)) || (!parser.isFileGoodForReading())) + std::ifstream infile(filePath); + if (!infile) return false; - // Get array dimensioos - if (parser.getArgsDelim() != LineParser::Success) + auto head = header.parse(infile); + if (!head) return false; - auto N = parser.argi(0); - // Get voxel sizes, assuming cubic grid - if (parser.getArgsDelim() != LineParser::Success) - return false; - auto delta = Vector3(parser.argd(0), parser.argd(4), parser.argd(8)); - - // Get grid origin coordinates - if (parser.getArgsDelim() != LineParser::Success) - return false; - auto axisOrigin = parser.arg3d(0); - - // Get loop order - we handle `zyx` and nothing else for now - if (parser.getArgsDelim() != LineParser::Success) - return false; - if (parser.args(0) != "zyx") - return Messenger::error("Only 'zyx' loop order is allowed.\n"); + auto &[n, imin, imax, a, b, c, axisOrigin] = std::get<0>(*head); // Set up our data - data.initialise(N, axisOrigin.x, delta.x, N, axisOrigin.y, delta.y, N, axisOrigin.z, delta.z); + data.initialise(n.x, axisOrigin.x, a.x, n.y, axisOrigin.y, b.y, n.z, axisOrigin.z, c.z); + auto points = some(maybe(spaces()) >> real() << spaces()).exact(infile); + if (!points) + return false; + auto idx = 0; // Loop over data values ('zyx' loop order, meaning fastest varying is z) - for (auto x = 0; x < N; ++x) - for (auto y = 0; y < N; ++y) - for (auto z = 0; z < N; ++z) - { - // Read line - if (parser.getArgsDelim() != LineParser::Success) - return false; - + for (auto x = 0; x < n.x; ++x) + for (auto y = 0; y < n.y; ++y) + for (auto z = 0; z < n.z; ++z) // Set the value - data.value(x, y, z) = parser.argd(0); - } + data.value(x, y, z) = (*points)[idx++]; return true; -} \ No newline at end of file +} diff --git a/src/nodes/importDLPUtilsSurface.cpp b/src/nodes/importDLPUtilsSurface.cpp index 90770f4616..4ca21e9ee1 100644 --- a/src/nodes/importDLPUtilsSurface.cpp +++ b/src/nodes/importDLPUtilsSurface.cpp @@ -2,7 +2,8 @@ // Copyright (c) 2026 Team Dissolve and contributors #include "nodes/importDLPUtilsSurface.h" -#include "base/lineParser.h" +#include "base/applicative.h" +#include "base/parserLibrary.h" ImportDLPUtilsSurfaceNode::ImportDLPUtilsSurfaceNode(Graph *parentGraph) : Node(parentGraph) { @@ -42,44 +43,34 @@ NodeConstants::ProcessResult ImportDLPUtilsSurfaceNode::process() // Read data specified bool ImportDLPUtilsSurfaceNode::read(Data2D &data, std::string filePath) { + using namespace Parsers; data.clear(); // Open file and check that we're OK to proceed importing from it - LineParser parser; - if ((!parser.openInput(filePath)) || (!parser.isFileGoodForReading())) + std::ifstream infile(filePath); + if (!infile) return false; // Data is in blocks of common Y value, three-columns: x y f(x,y) + auto block = some(maybe(inlineSpaces()) >> vector3() << newlines()); + auto blocks = some(block << spaces()); std::vector xAxis, yAxis, values; - auto firstLineOfBlock = true; - while (!parser.eofOrBlank()) - { - if (parser.getArgsDelim(LineParser::KeepBlanks) != LineParser::Success) - return false; - // Is this a blank line inbetween blocks? - if (parser.nArgs() == 0) - { - firstLineOfBlock = true; - continue; - } - - // If this is the first line of the block, re-start x-axis storage and push next y value - if (firstLineOfBlock) + auto result = blocks.exact(infile); + if (!result) + return false; + for (auto group : *result) + { + xAxis.clear(); + yAxis.push_back(group[0].y); + for (auto point : group) { - xAxis.clear(); - yAxis.push_back(parser.argd(1)); - firstLineOfBlock = false; + xAxis.push_back(point.x); + values.push_back(point.z); } - - // Store the x-axis value - xAxis.push_back(parser.argd(0)); - - // Store the value - values.push_back(parser.argd(2)); } - parser.closeFiles(); + infile.close(); // Validity check on number of points in loaded file if (xAxis.size() * yAxis.size() != values.size()) @@ -96,4 +87,4 @@ bool ImportDLPUtilsSurfaceNode::read(Data2D &data, std::string filePath) data.values()[{x, y}] = values[index++]; return true; -} \ No newline at end of file +} diff --git a/src/nodes/importMoscitoStructure.cpp b/src/nodes/importMoscitoStructure.cpp index fa01147cac..2ea9b58e29 100644 --- a/src/nodes/importMoscitoStructure.cpp +++ b/src/nodes/importMoscitoStructure.cpp @@ -2,6 +2,7 @@ // Copyright (c) 2026 Team Dissolve and contributors #include "nodes/importMoscitoStructure.h" +#include "base/parserLibrary.h" ImportMoscitoStructureNode::ImportMoscitoStructureNode(Graph *parentGraph) : Node(parentGraph) { @@ -47,58 +48,42 @@ NodeConstants::ProcessResult ImportMoscitoStructureNode::process() * * Units are: distance = nm, velocities = nm ps-1, forces = kJ mol-1 nm-1 */ + using namespace Parsers; + auto firstLine = inlineSpaces() >> vector3() << newlines(); + auto secondLine = inlineSpaces() >> natural() << newlines(); + auto thirdLine = maybe(inlines()) >> newlines(); + auto atom = alphanums() & inlineSpaces() >> natural() << spaces() & vector3() << spaces() & vector3() << spaces() & + vector3() << spaces(); + auto molecule = alphanums() << spaces() & natural() & inlineSpaces() >> natural() & + inlineSpaces() >> natural() << newlines() & some(atom); + auto fileStructure = firstLine & secondLine << thirdLine & some(molecule); + + std::ifstream infile(filePath_); + if (!infile) + return error("Couldn't open file '{}' for loading Moscito data.\n", filePath_); + + auto parsed = fileStructure.exact(infile); + if (!parsed) + return error("Couldn't parse file '{}' for loading Moscito data.\n", filePath_); + + auto &[box, nmolecules, molecules] = *parsed; // Clear storage objects structure_.clear(); forces_.clear(); - // Open file and check that we're OK to proceed importing from it - LineParser parser; - if ((!parser.openInput(filePath_)) || (!parser.isFileGoodForReading())) - return error("Couldn't open file '{}' for loading Moscito data.\n", filePath_); - - message(" --> Importing coordinates in Moscito (str) format...\n"); - // Read cell lengths - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - structure_.createBox(parser.arg3d(0), {90.0, 90.0, 90.0}); - - // Read nmolecules - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - auto nMolecules = parser.argi(0); - message(" --> Structure file contains {} molecules.\n", nMolecules); - - for (auto n = 0; n < nMolecules; ++n) + structure_.createBox(box, {90.0, 90.0, 90.0}); + assert(nmolecules == molecules.size()); + for (auto molecule : molecules) { - // Read and discard remark and molecule label lines - if (parser.skipLines(2) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - - // Get number of atoms in this molecule (second integer) - if (parser.getArgsDelim(LineParser::KeepBlanks) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - auto nAtoms = parser.argi(1); + auto &[molname, moltype, natoms, molindex, atoms] = molecule; + assert(natoms == atoms.size()); // Read in atom coordinates - for (auto i = 0; i < nAtoms; ++i) + for (auto atom : atoms) { - // Read atom label / index line - if (parser.getArgsDelim(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - auto name = parser.args(0); - - // Read coordinates (in nm) - // Coordinates are in fixed format (15.8e) with *no spacing between values* - if (parser.readNextLine(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - std::string line{parser.line()}; - structure_.addAtom(name, {std::stof(line.substr(0, 15)) * 10.0, std::stof(line.substr(15, 15)) * 10.0, - std::stof(line.substr(30)) * 10.0}); - - // Skip velocity line - if (parser.skipLines(1) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; + auto &[name, index, pos, vec, force] = atom; + structure_.addAtom(name, pos * 10.0); /* * Read forces (in kJ mol-1 nm-1) @@ -111,13 +96,35 @@ NodeConstants::ProcessResult ImportMoscitoStructureNode::process() * * Note: Forces are in fixed format (15.8e) with *no spacing between values* */ - if (parser.readNextLine(LineParser::Defaults) != LineParser::Success) - return NodeConstants::ProcessResult::Failed; - line = parser.line(); - forces_.emplace_back(std::stof(line.substr(0, 15)) * 10.0, std::stof(line.substr(15, 15)) * 10.0, - std::stof(line.substr(30)) * 10.0); + forces_.emplace_back(force * 10.0); } } return NodeConstants::ProcessResult::Success; } + +// parse the header of a moscito file +Parsers::Parser> ImportMoscitoStructureNode::header() +{ + using namespace Parsers; + auto firstLine = inlineSpaces() >> vector3() << newlines(); + auto secondLine = inlineSpaces() >> natural() << newlines(); + auto thirdLine = maybe(inlines()) >> newlines(); + return firstLine & secondLine << thirdLine; +} + +// parse an atom from a moscito file +Parsers::Parser> ImportMoscitoStructureNode::atom() +{ + using namespace Parsers; + return alphanums() & inlineSpaces() >> natural() << spaces() & vector3() << spaces() & vector3() << spaces() & + vector3() << spaces(); +} +// parse an molecule from a moscito file +Parsers::Parser>>> +ImportMoscitoStructureNode::molecule() +{ + using namespace Parsers; + return alphanums() << spaces() & natural() & inlineSpaces() >> natural() & inlineSpaces() >> natural() << newlines() & + some(ImportMoscitoStructureNode::atom()); +} diff --git a/src/nodes/importMoscitoStructure.h b/src/nodes/importMoscitoStructure.h index ba45cd12a5..4f7bd1210d 100644 --- a/src/nodes/importMoscitoStructure.h +++ b/src/nodes/importMoscitoStructure.h @@ -3,7 +3,9 @@ #pragma once +#include "base/applicative.h" #include "classes/structure.h" +#include "math/vector3.h" #include "nodes/node.h" class ImportMoscitoStructureNode : public Node @@ -38,4 +40,14 @@ class ImportMoscitoStructureNode : public Node protected: // Perform processing NodeConstants::ProcessResult process() override; -}; \ No newline at end of file + + public: + // Parse the header of a moscito file + static Parsers::Parser> header(); + // Parse an atom from a moscito file + static Parsers::Parser> atom(); + // Parse an molecule from a moscito file + static Parsers::Parser< + std::tuple>>> + molecule(); +}; diff --git a/src/nodes/md.cpp b/src/nodes/md.cpp index f279265a29..88ce066458 100644 --- a/src/nodes/md.cpp +++ b/src/nodes/md.cpp @@ -7,6 +7,7 @@ #include "kernels/force.h" #include "math/mathFunc.h" #include "nodes/dissolve.h" +#include MDNode::MDNode(Graph *parentGraph) : Node(parentGraph) { @@ -285,11 +286,13 @@ NodeConstants::ProcessResult MDNode::process() } // Open trajectory file (if requested) - LineParser trajParser; + std::ofstream outfile; + std::ostreambuf_iterator out(outfile); if (trajectoryFrequency && trajectoryFrequency > 0) { std::string trajectoryFile = std::format("{}.md.xyz", targetConfiguration_->name()); - if ((!trajParser.appendOutput(trajectoryFile)) || (!trajParser.isFileGoodForWriting())) + outfile.open(trajectoryFile, std::ios::app); + if (!outfile) { Messenger::error("Failed to open MD trajectory output file '{}'.\n", trajectoryFile); return NodeConstants::ProcessResult::Failed; @@ -421,30 +424,26 @@ NodeConstants::ProcessResult MDNode::process() if (trajectoryFrequency > 0 && (step % trajectoryFrequency == 0)) { // Write number of atoms - trajParser.writeLineF("{}\n", targetConfiguration_->nAtoms()); + std::format_to(out, "{}\n", targetConfiguration_->nAtoms()); // Construct and write header - std::string header = std::format("Step {} of {}, T = {:10.3e}, ke = {:10.3e}", step, nSteps, tInstant, ke); + std::format_to(out, "Step {} of {}, T = {:10.3e}, ke = {:10.3e}", step, nSteps, tInstant, ke); if (energyFrequency && (step % energyFrequency == 0)) - header += std::format(", inter = {:10.3e}, intra = {:10.3e}, tot = {:10.3e}", pe.pairPotential.total(), - pe.geometry.total(), ke + pe.pairPotential.total() + pe.geometry.total()); - if (!trajParser.writeLine(header)) - return NodeConstants::ProcessResult::Failed; + std::format_to(out, ", inter = {:10.3e}, intra = {:10.3e}, tot = {:10.3e}", pe.pairPotential.total(), + pe.geometry.total(), ke + pe.pairPotential.total() + pe.geometry.total()); + std::format_to(out, "\n"); // 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)) - return NodeConstants::ProcessResult::Failed; - } + std::format_to(out, "{:<3} {:10.3f} {:10.3f} {:10.3f}\n", Elements::symbol(i.speciesAtom()->Z()), i.r().x, + i.r().y, i.r().z); } } timer.stop(); // Close trajectory file if (trajectoryFrequency > 0) - trajParser.closeFiles(); + outfile.close(); if (capForces_) Messenger::print("A total of {} forces were capped over the course of the dynamics ({:9.3e} per step).\n", nCapped, diff --git a/tests/ff/cos-n.cpp b/tests/ff/cos-n.cpp index 60f39f7770..f9733d3898 100644 --- a/tests/ff/cos-n.cpp +++ b/tests/ff/cos-n.cpp @@ -1,9 +1,7 @@ // SPDX-License-Identifier: GPL-3.0-or-later // Copyright (c) 2026 Team Dissolve and contributors -#include "data/ff/library.h" #include "tests/graphData.h" -#include "tests/tempFile.h" #include "tests/testData.h" #include @@ -61,7 +59,7 @@ TEST(CosNTorsionForcesTest, POE) // Check agreement with external reference forces checkReferenceForceConsistency(pairPotentialForces, geometryForces, - importNode->getOutputValue>("Forces"), 6.0e-2); + importNode->getOutputValue>("Forces"), 6.2e-2); } TEST(CosNTorsionEnergyTest, Py4OHNTf2) diff --git a/tests/io/applicative.cpp b/tests/io/applicative.cpp index 7660ec5357..5bf819818a 100644 --- a/tests/io/applicative.cpp +++ b/tests/io/applicative.cpp @@ -3,6 +3,9 @@ #include "base/applicative.h" #include "base/parserLibrary.h" +#include "nodes/importDLPOLYStructure.h" +#include "nodes/importDLPOLYTrajectory.h" +#include "nodes/importMoscitoStructure.h" #include "nodes/importXYZStructure.h" #include #include @@ -158,4 +161,73 @@ TEST(ApplicativeTest, Helium) EXPECT_EQ(elem, "He"); } +TEST(ApplicativeTest, DLPOLYStructure) +{ + std::ifstream infile{"dlpoly/water1000/full.REVCON"}; + ASSERT_TRUE(infile); + auto header = ImportDLPOLYStructureNode::header().parse(infile); + ASSERT_TRUE(header); + auto &[value, rest] = *header; + auto [title, keytrj, imcon, natoms] = value; + EXPECT_TRUE(title.starts_with("Bulk")); + EXPECT_EQ(keytrj, 2); + EXPECT_EQ(imcon, 1); + EXPECT_EQ(natoms, 3000); + auto matrix = matrix3().parse(infile); + ASSERT_TRUE(matrix); + std::cout << "Next char:\t" << infile.peek() << std::endl; + auto example = ImportDLPOLYStructureNode::atom(); + auto more = some(example).parse(infile); + ASSERT_TRUE(more); + auto &[pos, vel, forces] = std::get<0>(*more)[0]; + EXPECT_EQ(pos.x, -10.97620028); + EXPECT_EQ(vel->y, -1.89534626301); + EXPECT_EQ(forces->z, -7400.12500402); +} + +TEST(ApplicativeTest, DLPOLYTrajectory) +{ + std::ifstream infile{"dlpoly/water267-npt/water-267-298K.HISf"}; + ASSERT_TRUE(infile); + auto header = ImportDLPOLYTrajectoryNode::header().parse(infile); + ASSERT_TRUE(header); + auto &[value, rest] = *header; + auto [stepno, natoms, keytrj, imcon, _] = value; + EXPECT_EQ(natoms, 801); + auto matrix = matrix3().parse(infile); + ASSERT_TRUE(matrix); + std::cout << "Next char:\t" << infile.peek() << std::endl; + auto example = ImportDLPOLYStructureNode::atom(); + auto more = some(example).parse(infile); + ASSERT_TRUE(more); + auto &[pos, vel, forces] = std::get<0>(*more)[0]; + EXPECT_EQ(pos.x, 8.278); +} + +TEST(ApplicativeTest, Moscito) +{ + std::ifstream infile{"moscito/py5_torsions/py5-ntf2-final.str"}; + ASSERT_TRUE(infile); + auto header = ImportMoscitoStructureNode::header().parse(infile); + ASSERT_TRUE(header); + auto &[c, nmolecules] = std::get<0>(*header); + EXPECT_EQ(c.x, 2.0); + EXPECT_EQ(nmolecules, 2); + + auto mols = some(ImportMoscitoStructureNode::molecule()).parse(infile); + ASSERT_TRUE(mols); + auto &[name, type, natoms, index, atoms] = std::get<0>(*mols)[0]; + EXPECT_EQ(name, "cat"); + EXPECT_EQ(type, 1); + EXPECT_EQ(natoms, 27); + EXPECT_EQ(index, 1); + + auto &[atomname, idx, pos, vel, force] = atoms[0]; + EXPECT_EQ(atomname, "nc1"); + EXPECT_EQ(idx, 1); + EXPECT_EQ(pos.x, 0.31811301); + EXPECT_EQ(vel.x, 0.43076653); + EXPECT_EQ(force.x, -0.12764812E+03); +} + } // namespace UnitTest