diff --git a/src/classes/structure.cpp b/src/classes/structure.cpp index 3a5d2db1f6..2c20f0d8a2 100644 --- a/src/classes/structure.cpp +++ b/src/classes/structure.cpp @@ -246,6 +246,82 @@ void Structure::createBox(const Matrix3 &axes) box_ = Box::generate(lengths, angles); } +/* + * Manipulations + */ + +// Recursive function for general manipulation +void Structure::recurseLocal(std::set &fragmentAtoms, StructureAtom *i, ManipulationFunction action) +{ + if (fragmentAtoms.contains(i)) + return; + + fragmentAtoms.insert(i); + + // Loop over attached atoms, performing minimum image repositioning w.r.t. i, and call the action + for (const auto *b : i->bonds()) + { + auto j = b->partner(i); + if (fragmentAtoms.contains(j)) + continue; + + action(j, box_.minimumImage(j->r(), i->r())); + + // Recurse into bound neighbours + recurseLocal(fragmentAtoms, j, action); + } +} +void Structure::recurseLocal(std::set &fragmentAtoms, StructureAtom *i, ConstManipulationFunction action) const +{ + if (fragmentAtoms.contains(i)) + return; + + fragmentAtoms.insert(i); + + // Loop over attached atoms, performing minimum image repositioning w.r.t. i, and call the action + for (const auto b : i->bonds()) + { + auto j = b->partner(i); + if (fragmentAtoms.contains(j)) + continue; + + action(j, box_.minimumImage(j->r(), i->r())); + + // Recurse into bound neighbours + recurseLocal(fragmentAtoms, j, action); + } +} + +// Return atoms in the same fragment as the specified atom, unfolding the fragment at the same time +std::set Structure::getUnfoldedFragment(StructureAtom *containing, ManipulationFunction action) +{ + std::set fragmentAtoms; + action(containing, containing->r()); + recurseLocal(fragmentAtoms, containing, action); + return fragmentAtoms; +} +std::set Structure::getUnfoldedFragment(StructureAtom *containing, ConstManipulationFunction action) const +{ + std::set fragmentAtoms; + action(containing, containing->r()); + recurseLocal(fragmentAtoms, containing, action); + return fragmentAtoms; +} + +// Un-fold bound fragments in the structure +void Structure::unFold() +{ + std::set fragmentAtoms; + + for (auto &atom : atoms()) + { + if (fragmentAtoms.contains(atom.get())) + break; + + fragmentAtoms.merge(getUnfoldedFragment(atom.get(), [](StructureAtom *j, Vector3 rJ) { j->setR(rJ); })); + } +} + /* * Serialisation */ diff --git a/src/classes/structure.h b/src/classes/structure.h index 5c9c359ca3..79f084ec61 100644 --- a/src/classes/structure.h +++ b/src/classes/structure.h @@ -7,6 +7,7 @@ #include "classes/atom.h" #include "classes/bond.h" #include "classes/box.h" +#include #include // StructureAtom @@ -116,6 +117,24 @@ class Structure : public Serialisable<> // Create Box definition from axes matrix void createBox(const Matrix3 &axes); + /* + * Manipulations + */ + private: + // Typedef for manipulation functions + using ManipulationFunction = std::function; + using ConstManipulationFunction = std::function; + // Recursive function for general manipulation + void recurseLocal(std::set &fragmentAtoms, StructureAtom *i, ManipulationFunction action); + void recurseLocal(std::set &fragmentAtoms, StructureAtom *i, ConstManipulationFunction action) const; + // Return atoms in the same fragment as the specified atom, unfolding the fragment at the same time + std::set getUnfoldedFragment(StructureAtom *containing, ManipulationFunction action); + std::set getUnfoldedFragment(StructureAtom *containing, ConstManipulationFunction action) const; + + public: + // Un-fold bound fragments in the structure + void unFold(); + /* * Serialisation */ diff --git a/src/nodes/registry.cpp b/src/nodes/registry.cpp index d2a9ded0e2..695042afa4 100644 --- a/src/nodes/registry.cpp +++ b/src/nodes/registry.cpp @@ -56,6 +56,7 @@ #include "nodes/species.h" #include "nodes/sq.h" #include "nodes/subtract.h" +#include "nodes/supercellConfiguration.h" #include "nodes/vector3Assemble.h" #include "nodes/vector3Decompose.h" #include "nodes/voxelDensity.h" @@ -133,6 +134,7 @@ void NodeRegistry::instantiateNodeProducers() {"SQ", makeDerivedNode()}, {"Species", makeDerivedNode()}, {"Subtract", makeDerivedNode()}, + {"SupercellConfiguration", makeDerivedNode()}, {"Vector3Assemble", makeDerivedNode()}, {"Vector3Decompose", makeDerivedNode()}, {"XRaySQ", makeDerivedNode()}, diff --git a/src/nodes/supercellConfiguration.cpp b/src/nodes/supercellConfiguration.cpp new file mode 100644 index 0000000000..213f7d4b0e --- /dev/null +++ b/src/nodes/supercellConfiguration.cpp @@ -0,0 +1,69 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2026 Team Dissolve and contributors + +#include "nodes/supercellConfiguration.h" +#include "math/vector3.h" + +SupercellConfigurationNode::SupercellConfigurationNode(Graph *parentGraph) : Node(parentGraph) +{ + // Inputs + addInput("Configuration", "Target configuration", targetConfiguration_); + + // Outputs + addPointerOutput("SupercellConfiguration", "Supercell configuration", supercellConfiguration_); + + // Options + addOption("SupercellRepeat", "Integer coefficients by which unit cell will be repeated along its axes", supercellRepeat_); +} + +/* + * Definition + */ + +std::string_view SupercellConfigurationNode::type() const { return "SupercellConfiguration"; } + +std::string_view SupercellConfigurationNode::summary() const +{ + return "Create a repeated instance (supercell) of a configuration"; +} + +/* + * Processing + */ + +// Perform processing +NodeConstants::ProcessResult SupercellConfigurationNode::process() +{ + supercellConfiguration_.empty(); + + const auto &box = targetConfiguration_->box(); + + // Set up configuration + auto supercellLengths = box.axisLengths(); + supercellLengths.multiply(supercellRepeat_.x, supercellRepeat_.y, supercellRepeat_.z); + supercellConfiguration_.createBoxAndCells(supercellLengths, box.axisAngles(), false); + + // Create images of all molecular unit cell species + for (auto &mol : targetConfiguration_->molecules()) + { + const auto *sp = mol->species(); + + // Loop over cell images + for (auto ix = 0; ix < supercellRepeat_.x; ++ix) + { + for (auto iy = 0; iy < supercellRepeat_.y; ++iy) + { + // Create and translate molecule + for (auto iz = 0; iz < supercellRepeat_.z; ++iz) + supercellConfiguration_.addMolecule(sp)->translate(box.axes() * Vector3(ix, iy, iz)); + } + } + } + + supercellConfiguration_.updateObjectRelationships(); + + message("Created ({}, {}, {}) supercell - {} atoms total.\n", supercellRepeat_.x, supercellRepeat_.y, supercellRepeat_.z, + supercellConfiguration_.nAtoms()); + + return NodeConstants::ProcessResult::Success; +} \ No newline at end of file diff --git a/src/nodes/supercellConfiguration.h b/src/nodes/supercellConfiguration.h new file mode 100644 index 0000000000..02b66a4113 --- /dev/null +++ b/src/nodes/supercellConfiguration.h @@ -0,0 +1,39 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2026 Team Dissolve and contributors + +#pragma once + +#include "classes/configuration.h" +#include "nodes/node.h" + +class SupercellConfigurationNode : public Node +{ + public: + SupercellConfigurationNode(Graph *parentGraph); + ~SupercellConfigurationNode() override = default; + + /* + * Definition + */ + public: + std::string_view type() const override; + std::string_view summary() const override; + + /* + * Data + */ + private: + // Input configuration + Configuration *targetConfiguration_{nullptr}; + // Supercell repeat + Vector3i supercellRepeat_; + // Supercell configuration + Configuration supercellConfiguration_; + + /* + * Processing + */ + protected: + // Perform processing + NodeConstants::ProcessResult process() override; +}; diff --git a/tests/classes/structure.cpp b/tests/classes/structure.cpp index f45717d007..e07c5032b5 100644 --- a/tests/classes/structure.cpp +++ b/tests/classes/structure.cpp @@ -3,6 +3,8 @@ #include "classes/structure.h" #include "nodes/calculateBonding.h" +#include "tests/graphData.h" +#include #include namespace UnitTest @@ -68,4 +70,37 @@ TEST(StructureTest, Molecule2) EXPECT_EQ(s.nAtoms(), 0); } +TEST(StructureTest, Unfold) +{ + TestGraph graph; + + // Import XYZ node + auto importXYZNode = graph.createNode("ImportXYZStructure"); + ASSERT_TRUE(importXYZNode); + ASSERT_TRUE(importXYZNode->setOption("FilePath", std::string("xyz/ch4_folded.xyz"))); + + // Calculate bonding node + auto calculateBondingNode = graph.createNode("CalculateBonding"); + ASSERT_TRUE(calculateBondingNode); + + // Connect graph + ASSERT_TRUE(graph.addEdge({"ImportXYZStructure", "Structure", "CalculateBonding", "Structure"})); + + // Run + ASSERT_EQ(calculateBondingNode->run(), NodeConstants::ProcessResult::Success); + + auto structure = calculateBondingNode->getOutputValue("Structure"); + ASSERT_TRUE(structure.bonds().size() != 0); + structure.unFold(); + + // After unfolding, the distances between C and H should all be unity + for (const auto &bond : structure.bonds()) + { + auto iR = bond->i()->r(); + auto jR = bond->j()->r(); + auto distance = abs((iR - jR).magnitude()); + ASSERT_NEAR(distance, 1, 10e-9); + } +} + } // namespace UnitTest diff --git a/tests/data/xyz/ch4_folded.xyz b/tests/data/xyz/ch4_folded.xyz new file mode 100644 index 0000000000..e59f7b7701 --- /dev/null +++ b/tests/data/xyz/ch4_folded.xyz @@ -0,0 +1,17 @@ +15 +Unnamed001 +C 5.000000 0.100000 5.000000 0.000000 +H 5.000000 1.100000 5.000000 0.000000 +H 5.000000 9.766193 4.057359 0.000000 +H 5.816351 9.766193 5.471321 0.000000 +H 4.183649 9.766193 5.471321 0.000000 +C 0.100000 5.000000 5.000000 0.000000 +H 1.100000 5.000000 5.000000 0.000000 +H 9.766193 5.000000 5.942641 0.000000 +H 9.766193 5.816351 4.528679 0.000000 +H 9.766193 4.183649 4.528679 0.000000 +C 5.000000 5.000000 0.100000 0.000000 +H 5.000000 5.000000 1.100000 0.000000 +H 5.000000 5.942641 9.766193 0.000000 +H 5.816351 4.528679 9.766193 0.000000 +H 4.183649 4.528679 9.766193 0.000000