Add write support for molden

avogadro/quantumio/molden.cpp

@@ -5,11 +5,14 @@
 
 #include "molden.h"
 
+#include <avogadro/core/elements.h>
 #include <avogadro/core/gaussianset.h>
 #include <avogadro/core/molecule.h>
 #include <avogadro/core/utilities.h>
 
+#include <iomanip>
 #include <iostream>
+#include <sstream>
 
 using std::cout;
 using std::endl;
@@ -455,4 +458,301 @@ void MoldenFile::outputAll()
     cout << m_MOcoeff << "\t";
   cout << endl;
 }
+
+bool MoldenFile::write(std::ostream& out, const Core::Molecule& molecule)
+{
+  // Set output precision for floating point numbers
+  out << std::setprecision(10) << std::scientific;
+
+  // Write the Molden format header
+  out << "[Molden Format]\n";
+
+  // Get the basis set if available
+  const auto* basis = dynamic_cast<const GaussianSet*>(molecule.basisSet());
+
+  // Check for spherical basis functions and write appropriate flags
+  if (basis != nullptr) {
+    // We need to check what types of shells are present
+    std::vector<int> symmetry = basis->symmetry();
+    bool hasD5 = false, hasF7 = false, hasG9 = false;
+
+    for (int sym : symmetry) {
+      if (sym == GaussianSet::D5)
+        hasD5 = true;
+      else if (sym == GaussianSet::F7)
+        hasF7 = true;
+      else if (sym == GaussianSet::G9)
+        hasG9 = true;
+    }
+
+    if (hasD5)
+      out << "[5D]\n";
+    if (hasF7)
+      out << "[7F]\n";
+    if (hasG9)
+      out << "[9G]\n";
+  }
+
+  // Write atoms section
+  writeAtoms(out, molecule);
+
+  // Write GTO basis set if available
+  if (basis != nullptr) {
+    writeGTO(out, basis);
+    writeMO(out, basis);
+  }
+
+  // Write multiple coordinate sets if available
+  if (molecule.coordinate3dCount() > 1) {
+    writeGeometries(out, molecule);
+  }
+
+  // Write vibrational data if available
+  if (molecule.vibrationFrequencies().size() > 0) {
+    writeFrequencies(out, molecule);
+  }
+
+  return true;
+}
+
+void MoldenFile::writeAtoms(std::ostream& out, const Core::Molecule& molecule)
+{
+  out << "[Atoms] Angs\n";
+
+  for (Index i = 0; i < molecule.atomCount(); ++i) {
+    unsigned char atomicNum = molecule.atomicNumber(i);
+    Vector3 pos = molecule.atomPosition3d(i);
+
+    // Format: element_name number atomic_number x y z
+    out << std::left << std::setw(4) << Core::Elements::symbol(atomicNum)
+        << std::right << std::setw(6) << (i + 1) << std::setw(6)
+        << static_cast<int>(atomicNum) << "  " << std::setw(18) << pos.x()
+        << std::setw(18) << pos.y() << std::setw(18) << pos.z() << "\n";
+  }
+}
+
+void MoldenFile::writeGTO(std::ostream& out, const Core::GaussianSet* basis)
+{
+  out << "[GTO]\n";
+
+  std::vector<int> symmetry = basis->symmetry();
+  std::vector<unsigned int> atomIndices = basis->atomIndices();
+  std::vector<unsigned int> gtoIndices = basis->gtoIndices();
+  std::vector<double> gtoA = basis->gtoA();
+  std::vector<double> gtoC = basis->gtoC();
+
+  // Group shells by atom
+  unsigned int currentAtom = UINT_MAX;
+  for (unsigned int i = 0; i < symmetry.size(); ++i) {
+    unsigned int atomIdx = atomIndices[i];
+
+    // Start a new atom block if needed
+    if (atomIdx != currentAtom) {
+      if (currentAtom != UINT_MAX) {
+        out << "\n"; // Blank line between atoms
+      }
+      currentAtom = atomIdx;
+      out << (atomIdx + 1)
+          << " 0\n"; // atom number (1-indexed) and unused field
+    }
+
+    // Get shell type string
+    string shellType;
+    int shellSym = symmetry[i];
+    switch (shellSym) {
+      case GaussianSet::S:
+        shellType = "s";
+        break;
+      case GaussianSet::P:
+        shellType = "p";
+        break;
+      case GaussianSet::D:
+      case GaussianSet::D5:
+        shellType = "d";
+        break;
+      case GaussianSet::F:
+      case GaussianSet::F7:
+        shellType = "f";
+        break;
+      case GaussianSet::G:
+      case GaussianSet::G9:
+        shellType = "g";
+        break;
+      case GaussianSet::H:
+      case GaussianSet::H11:
+        shellType = "h";
+        break;
+      case GaussianSet::I:
+      case GaussianSet::I13:
+        shellType = "i";
+        break;
+      default:
+        shellType = "s"; // fallback
+        break;
+    }
+
+    // Get the number of primitives in this shell
+    unsigned int startGTO = gtoIndices[i];
+    unsigned int endGTO =
+      (i + 1 < gtoIndices.size()) ? gtoIndices[i + 1] : gtoA.size();
+    unsigned int numPrimitives = endGTO - startGTO;
+
+    out << shellType << " " << numPrimitives << " 1.00\n";
+
+    // Write exponents and contraction coefficients
+    for (unsigned int j = startGTO; j < endGTO; ++j) {
+      out << "  " << std::setw(18) << gtoA[j] << std::setw(18) << gtoC[j]
+          << "\n";
+    }
+  }
+  out << "\n"; // Final blank line to end GTO section
+}
+
+void MoldenFile::writeMO(std::ostream& out, const Core::GaussianSet* basis)
+{
+  out << "[MO]\n";
+
+  Core::ScfType scfType = basis->scfType();
+  bool isOpenShell = (scfType == Core::Uhf || scfType == Core::Rohf);
+
+  // Helper lambda to write orbitals for a given electron type
+  auto writeOrbitals = [&](BasisSet::ElectronType type,
+                           const string& spinLabel) {
+    MatrixX moMatrix = basis->moMatrix(type);
+    std::vector<double> energies = basis->moEnergy(type);
+    std::vector<std::string> symLabels = basis->symmetryLabels(type);
+    std::vector<unsigned char> occupancies = basis->moOccupancy(type);
+
+    if (moMatrix.cols() == 0)
+      return;
+
+    unsigned int numMOs = moMatrix.cols();
+    unsigned int numBasis = moMatrix.rows();
+
+    for (unsigned int mo = 0; mo < numMOs; ++mo) {
+      // Write symmetry label
+      string symLabel = (mo < symLabels.size() && !symLabels[mo].empty())
+                          ? symLabels[mo]
+                          : "a1";
+      out << " Sym= " << symLabel << "\n";
+
+      // Write energy (convert from eV back to Hartree for Molden format)
+      double energy =
+        (mo < energies.size()) ? energies[mo] / HARTREE_TO_EV_D : 0.0;
+      out << " Ene= " << std::setw(18) << energy << "\n";
+
+      // Write spin
+      out << " Spin= " << spinLabel << "\n";
+
+      // Write occupation
+      int occup = 0;
+      if (mo < occupancies.size()) {
+        occup = occupancies[mo];
+      } else if (!isOpenShell) {
+        // For closed shell, assume doubly occupied up to HOMO
+        occup = (mo < basis->lumo(type)) ? 2 : 0;
+      } else {
+        // For open shell, assume singly occupied up to HOMO
+        occup = (mo < basis->lumo(type)) ? 1 : 0;
+      }
+      out << " Occup= " << occup << "\n";
+
+      // Write MO coefficients
+      for (unsigned int bf = 0; bf < numBasis; ++bf) {
+        out << std::setw(6) << (bf + 1) << std::setw(18) << moMatrix(bf, mo)
+            << "\n";
+      }
+    }
+  };
+
+  if (isOpenShell) {
+    // Write alpha and beta orbitals separately
+    writeOrbitals(BasisSet::Alpha, "Alpha");
+    writeOrbitals(BasisSet::Beta, "Beta");
+  } else {
+    // Write paired orbitals as Alpha (standard convention)
+    writeOrbitals(BasisSet::Paired, "Alpha");
+  }
+}
+
+void MoldenFile::writeFrequencies(std::ostream& out,
+                                  const Core::Molecule& molecule)
+{
+  Core::Array<double> frequencies = molecule.vibrationFrequencies();
+  Core::Array<double> irIntensities = molecule.vibrationIRIntensities();
+  Core::Array<double> ramanIntensities = molecule.vibrationRamanIntensities();
+
+  if (frequencies.size() == 0)
+    return;
+
+  // Write frequencies
+  out << "[FREQ]\n";
+  for (size_t i = 0; i < frequencies.size(); ++i) {
+    out << std::setw(18) << frequencies[i] << "\n";
+  }
+
+  // Write coordinates for vibrations (in Bohr)
+  out << "[FR-COORD]\n";
+  for (Index i = 0; i < molecule.atomCount(); ++i) {
+    unsigned char atomicNum = molecule.atomicNumber(i);
+    Vector3 pos = molecule.atomPosition3d(i) * ANGSTROM_TO_BOHR_D;
+
+    out << std::left << std::setw(4) << Core::Elements::symbol(atomicNum)
+        << std::right << std::setw(18) << pos.x() << std::setw(18) << pos.y()
+        << std::setw(18) << pos.z() << "\n";
+  }
+
+  // Write normal mode displacements (in Bohr)
+  out << "[FR-NORM-COORD]\n";
+  for (size_t mode = 0; mode < frequencies.size(); ++mode) {
+    out << "vibration " << (mode + 1) << "\n";
+    Core::Array<Vector3> lx = molecule.vibrationLx(static_cast<int>(mode));
+    for (size_t atom = 0; atom < lx.size(); ++atom) {
+      // Convert from Angstrom to Bohr for Molden format
+      Vector3 disp = lx[atom] * ANGSTROM_TO_BOHR_D;
+      out << std::setw(18) << disp.x() << std::setw(18) << disp.y()
+          << std::setw(18) << disp.z() << "\n";
+    }
+  }
+
+  // Write intensities if available
+  if (irIntensities.size() > 0) {
+    out << "[INT]\n";
+    for (size_t i = 0; i < irIntensities.size(); ++i) {
+      out << std::setw(18) << irIntensities[i];
+      if (i < ramanIntensities.size()) {
+        out << std::setw(18) << ramanIntensities[i];
+      }
+      out << "\n";
+    }
+  }
+}
+
+void MoldenFile::writeGeometries(std::ostream& out,
+                                 const Core::Molecule& molecule)
+{
+  size_t numCoordSets = molecule.coordinate3dCount();
+  if (numCoordSets <= 1)
+    return;
+
+  out << "[GEOMETRIES] XYZ\n";
+
+  for (size_t coordIdx = 0; coordIdx < numCoordSets; ++coordIdx) {
+    Core::Array<Vector3> coords = molecule.coordinate3d(coordIdx);
+
+    // Write number of atoms and a comment line (required for XYZ format)
+    out << molecule.atomCount() << "\n";
+    out << "Frame " << (coordIdx + 1) << "\n";
+
+    for (Index i = 0; i < molecule.atomCount(); ++i) {
+      unsigned char atomicNum = molecule.atomicNumber(i);
+      Vector3 pos = (i < coords.size()) ? coords[i] : Vector3::Zero();
+
+      out << std::left << std::setw(4) << Core::Elements::symbol(atomicNum)
+          << std::right << std::setw(18) << pos.x() << std::setw(18) << pos.y()
+          << std::setw(18) << pos.z() << "\n";
+    }
+  }
+}
+
 } // namespace Avogadro::QuantumIO

avogadro/quantumio/molden.h

@@ -24,7 +24,7 @@ class AVOGADROQUANTUMIO_EXPORT MoldenFile : public Io::FileFormat
 
   Operations supportedOperations() const override
   {
-    return Read | File | Stream | String;
+    return Read | Write | File | Stream | String;
   }
 
   FileFormat* newInstance() const override { return new MoldenFile; }
@@ -41,11 +41,8 @@ class AVOGADROQUANTUMIO_EXPORT MoldenFile : public Io::FileFormat
   std::vector<std::string> mimeTypes() const override;
 
   [[nodiscard]] bool read(std::istream& in, Core::Molecule& molecule) override;
-  [[nodiscard]] bool write(std::ostream&, const Core::Molecule&) override
-  {
-    // Empty, as we do not write out Molden files.
-    return false;
-  }
+  [[nodiscard]] bool write(std::ostream& out,
+                           const Core::Molecule& molecule) override;
 
 private:
   void outputAll();
@@ -54,6 +51,13 @@ class AVOGADROQUANTUMIO_EXPORT MoldenFile : public Io::FileFormat
   void readAtom(const std::vector<std::string>& list);
   void load(Core::GaussianSet* basis);
 
+  // Write helper methods
+  void writeAtoms(std::ostream& out, const Core::Molecule& molecule);
+  void writeGTO(std::ostream& out, const Core::GaussianSet* basis);
+  void writeMO(std::ostream& out, const Core::GaussianSet* basis);
+  void writeFrequencies(std::ostream& out, const Core::Molecule& molecule);
+  void writeGeometries(std::ostream& out, const Core::Molecule& molecule);
+
   bool m_cartesianD = true;
   bool m_cartesianF = true;
   bool m_cartesianG = true;