From 3caa5b074b4de43790c4beb93c244df25b2a40af Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 14 Oct 2024 23:03:37 -0400 Subject: [PATCH 1/3] Initial support for reading charges and forces from extended XYZ Signed-off-by: Geoff Hutchison --- avogadro/io/xyzformat.cpp | 52 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/avogadro/io/xyzformat.cpp b/avogadro/io/xyzformat.cpp index 9c3e391c42..f063722141 100644 --- a/avogadro/io/xyzformat.cpp +++ b/avogadro/io/xyzformat.cpp @@ -87,6 +87,36 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) mol.setUnitCell(cell); } } + // check to see if there's an extended XYZ Properties= line + // e.g. Properties=species:S:1:pos:R:3 + // https://gitlab.com/ase/ase/-/merge_requests/62 + start = buffer.find("Properties="); + unsigned int chargeColumn = 0; + unsigned int forceColumn = 0; + vector charges; + if (start != std::string::npos) { + start = start + 11; // skip over "Properties=" + unsigned int stop = buffer.find(' ', start); + unsigned int length = stop - start; + // we want to track columns after the position + // (esp. charge, spin, force, velocity, etc.) + std::string properties = buffer.substr(start, length); + vector tokens(split(properties, ':')); + unsigned int column = 0; + for (size_t i = 0; i < tokens.size(); i += 3) { + // increment column based on the count of the property + if (i + 2 < tokens.size()) { + column += lexicalCast(tokens[i + 2]); + } + + // we can safely assume species and pos are present + if (tokens[i] == "charge") { + chargeColumn = column; + } else if (tokens[i] == "force" || tokens[i] == "forces") { + forceColumn = column; + } // TODO other properties (velocity, spin, selection, etc.) + } + } // Parse atoms for (size_t i = 0; i < numAtoms; ++i) { @@ -114,6 +144,18 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) Atom newAtom = mol.addAtom(atomicNum); newAtom.setPosition3d(pos); + + // check for charge and force columns + if (chargeColumn > 0 && chargeColumn < tokens.size()) { + charges.push_back(lexicalCast(tokens[chargeColumn])); + // we set the charges after all atoms are added + } + if (forceColumn > 0 && forceColumn < tokens.size()) { + Vector3 force(lexicalCast(tokens[forceColumn]), + lexicalCast(tokens[forceColumn + 1]), + lexicalCast(tokens[forceColumn + 2])); + newAtom.setForceVector(force); + } } // Check that all atoms were handled. @@ -126,6 +168,15 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) return false; } + // set the charges + if (!charges.empty()) { + MatrixX chargesMatrix = MatrixX::Zero(mol.atomCount(), 1); + for (size_t i = 0; i < charges.size(); ++i) { + chargesMatrix(i, 0) = charges[i]; + } + mol.setPartialCharges("From File", chargesMatrix); + } + // Do we have an animation? size_t numAtoms2; // check if the next frame has the same number of atoms @@ -251,6 +302,7 @@ std::vector XyzFormat::fileExtensions() const { std::vector ext; ext.emplace_back("xyz"); + ext.emplace_back("exyz"); ext.emplace_back("extxyz"); return ext; } From 67033219b9320870ad27ccf374b04ee4057c0dca Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Tue, 15 Oct 2024 17:34:35 -0400 Subject: [PATCH 2/3] Added some debugging info - charges are set but disappear Signed-off-by: Geoff Hutchison --- avogadro/io/xyzformat.cpp | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/avogadro/io/xyzformat.cpp b/avogadro/io/xyzformat.cpp index f063722141..21378c9715 100644 --- a/avogadro/io/xyzformat.cpp +++ b/avogadro/io/xyzformat.cpp @@ -104,17 +104,18 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) vector tokens(split(properties, ':')); unsigned int column = 0; for (size_t i = 0; i < tokens.size(); i += 3) { - // increment column based on the count of the property - if (i + 2 < tokens.size()) { - column += lexicalCast(tokens[i + 2]); - } - // we can safely assume species and pos are present if (tokens[i] == "charge") { chargeColumn = column; + std::cout << " found charges " << chargeColumn << std::endl; } else if (tokens[i] == "force" || tokens[i] == "forces") { forceColumn = column; } // TODO other properties (velocity, spin, selection, etc.) + + // increment column based on the count of the property + if (i + 2 < tokens.size()) { + column += lexicalCast(tokens[i + 2]); + } } } @@ -175,6 +176,9 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) chargesMatrix(i, 0) = charges[i]; } mol.setPartialCharges("From File", chargesMatrix); + + std::cout << "Charges set " << charges.size() << " " << mol.atomCount() + << std::endl; } // Do we have an animation? @@ -245,6 +249,12 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) mol.perceiveBondOrders(); } + // check partial charge types + auto types = mol.partialChargeTypes(); + for (const auto& type : types) { + std::cout << "Partial charge type: " << type << std::endl; + } + return true; } From 1afc41923eaad02607f758da08770fa4b24fe2fa Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Tue, 15 Oct 2024 21:50:13 -0400 Subject: [PATCH 3/3] Fix up partial charges -- need to add after setting bonds Signed-off-by: Geoff Hutchison --- avogadro/io/xyzformat.cpp | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/avogadro/io/xyzformat.cpp b/avogadro/io/xyzformat.cpp index 21378c9715..023fa5ceda 100644 --- a/avogadro/io/xyzformat.cpp +++ b/avogadro/io/xyzformat.cpp @@ -107,7 +107,6 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) // we can safely assume species and pos are present if (tokens[i] == "charge") { chargeColumn = column; - std::cout << " found charges " << chargeColumn << std::endl; } else if (tokens[i] == "force" || tokens[i] == "forces") { forceColumn = column; } // TODO other properties (velocity, spin, selection, etc.) @@ -169,18 +168,6 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) return false; } - // set the charges - if (!charges.empty()) { - MatrixX chargesMatrix = MatrixX::Zero(mol.atomCount(), 1); - for (size_t i = 0; i < charges.size(); ++i) { - chargesMatrix(i, 0) = charges[i]; - } - mol.setPartialCharges("From File", chargesMatrix); - - std::cout << "Charges set " << charges.size() << " " << mol.atomCount() - << std::endl; - } - // Do we have an animation? size_t numAtoms2; // check if the next frame has the same number of atoms @@ -249,10 +236,14 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) mol.perceiveBondOrders(); } - // check partial charge types - auto types = mol.partialChargeTypes(); - for (const auto& type : types) { - std::cout << "Partial charge type: " << type << std::endl; + // have to set the charges after creating bonds + // (since modifying bonds invalidates the partial charges) + if (!charges.empty()) { + MatrixX chargesMatrix = MatrixX::Zero(mol.atomCount(), 1); + for (size_t i = 0; i < charges.size(); ++i) { + chargesMatrix(i, 0) = charges[i]; + } + mol.setPartialCharges("From File", chargesMatrix); } return true;