add function to generate tikz mo occupation images

Signed-off-by: Conrad Hübler <[email protected]>

src/capabilities/confscan.cpp

@@ -1113,8 +1113,9 @@ void ConfScan::Reorder(double dLE, double dLI, double dLH, bool reuse_only, bool
                     threads[i]->addReorderRule(rules[j]);
             }
 
-            if (m_RMSDmethod.compare("molalign") != 0 || m_threads == 1) {
-                p->StaticPool();
+            if (m_RMSDmethod.compare("molalign") != 0 || m_threads != 1) {
+                if (m_threads > 2)
+                    p->StaticPool();
                 p->StartAndWait();
             } else {
                 for (auto* t : threads) {

src/capabilities/curcumaopt.cpp

@@ -151,6 +151,10 @@ void CurcumaOpt::LoadControlJson()
     m_lambda = Json2KeyWord<double>(m_defaults, "lambda");
     m_diis_hist = Json2KeyWord<int>(m_defaults, "diis_hist");
     m_diis_start = Json2KeyWord<int>(m_defaults, "diis_start");
+    m_mo_scheme = Json2KeyWord<bool>(m_defaults, "mo_scheme");
+    m_mo_scale = Json2KeyWord<double>(m_defaults, "mo_scale");
+    m_mo_homo = Json2KeyWord<int>(m_defaults, "mo_homo");
+    m_mo_lumo = Json2KeyWord<int>(m_defaults, "mo_lumo");
 
     if (m_optimethod == 0) {
         std::cout << "Using external lBFGS module" << std::endl;
@@ -319,7 +323,7 @@ void CurcumaOpt::clear()
 
 double CurcumaOpt::SinglePoint(const Molecule* initial, std::string& output, std::vector<double>& charges)
 {
-    std::string method = Json2KeyWord<std::string>(m_controller, "method");
+    std::string method = m_method; // Json2KeyWord<std::string>(m_controller, "method");
 
     Geometry geometry = initial->getGeometry();
     Molecule tmp(initial);
@@ -348,9 +352,67 @@ double CurcumaOpt::SinglePoint(const Molecule* initial, std::string& output, std
     }
 #endif
     charges = interface.Charges();
+    if (m_mo_scheme) {
+        m_orbital_energies = interface.Energies();
+        m_num_electrons = interface.NumElectrons();
+        WriteMO(m_mo_homo, m_mo_lumo);
+        WriteMOAscii();
+    }
     return energy;
 }
 
+void CurcumaOpt::WriteMO(int n, int m)
+{
+    std::ofstream file(Basename() + ".inc");
+    std::ostream& out = std::cout;
+    auto write_tikz = [&](std::ostream& os) {
+        os << "\\begin{tikzpicture}\n";
+        int highest_occupied_orbital = m_num_electrons / 2 - 1;
+        int start_orbital = std::max(0, highest_occupied_orbital - m);
+        int end_orbital = std::min(static_cast<int>(m_orbital_energies.size()), highest_occupied_orbital + n + 1);
+
+        for (int i = start_orbital; i < end_orbital; ++i) {
+            double energy = m_orbital_energies[i] * m_mo_scale;
+            os << std::fixed << std::setprecision(4);
+            os << "\t\\draw[thick] (1," << energy << ") -- (0.0," << energy << ");\n";
+            os << "\t\\node at (1.1," << energy << ") {\\tiny " << m_orbital_energies[i] << "};\n";
+
+            if (i <= highest_occupied_orbital) {
+                os << "\t\\draw[thick,<-, color=orange] (0.25," << energy << ") -- (0.25," << energy << ");\n";
+                os << "\t\\draw[thick,->, color=orange] (0.75," << energy << ") -- (0.75," << energy << ");\n";
+            }
+        }
+        os << "\t\\node at (0.5,-1) {{" + Basename() + "}};\n";
+        os << "\t\\node at (0.5,-1.25) {\\hphantom{\\tiny (" + Basename() + ")}};\n";
+        os << "\\end{tikzpicture}\n";
+    };
+
+    write_tikz(out);
+    write_tikz(file);
+
+    file.close();
+}
+
+void CurcumaOpt::WriteMOAscii()
+{
+    std::vector<std::string> levels;
+    for (size_t i = 0; i < m_orbital_energies.size(); ++i) {
+        double energy = m_orbital_energies[i] * m_mo_scale;
+        std::string level = std::to_string(m_orbital_energies[i]);
+        level.resize(10, ' '); // Adjust the width for alignment
+        if (i < m_num_electrons / 2) {
+            levels.push_back("|" + std::string(8, '-') + "|" + level + "|<--->|");
+        } else {
+            levels.push_back("|" + std::string(8, '-') + "|" + level + "|     |");
+        }
+    }
+
+    std::reverse(levels.begin(), levels.end()); // Reverse to print from top to bottom
+    for (const auto& level : levels) {
+        std::cout << level << std::endl;
+    }
+}
+
 Molecule CurcumaOpt::LBFGSOptimise(Molecule* initial, std::string& output, std::vector<Molecule>* intermediate, std::vector<double>& charges, int thread, const std::string& basename)
 {
     std::vector<int> constrain;

src/capabilities/curcumaopt.h

@@ -66,7 +66,12 @@ static json CurcumaOptJson{
     { "inithess", false },
     { "lambda", 0.1 },
     { "diis_hist", 5 },
-    { "diis_start", 5 }
+    { "diis_start", 5 },
+    { "mo_scheme", false },
+    { "mo_scale", 1.0 },
+    { "mo_homo", -1 },
+    { "mo_lumo", -1 }
+
 };
 
 const json OptJsonPrivate{
@@ -230,6 +235,8 @@ class CurcumaOpt : public CurcumaMethod {
     double SinglePoint(const Molecule* initial, std::string& output, std::vector<double>& charges);
 
     void clear();
+    void WriteMO(int n, int m);
+    void WriteMOAscii();
 
 private:
     /* Lets have this for all modules */
@@ -254,11 +261,14 @@ class CurcumaOpt : public CurcumaMethod {
     Molecule m_molecule;
     json m_parameters;
     std::vector<Molecule> m_molecules;
-    bool m_file_set = false, m_mol_set = false, m_mols_set = false, m_writeXYZ = true, m_printoutput = true, m_singlepoint = false, m_fusion = false, m_optH = false, m_inithess = false;
-    int m_hessian = 0;
+    Vector m_orbital_energies;
+    Matrix m_molecular_orbitals;
+
+    bool m_file_set = false, m_mol_set = false, m_mols_set = false, m_writeXYZ = true, m_printoutput = true, m_singlepoint = false, m_fusion = false, m_optH = false, m_inithess = false, m_mo_scheme = false;
+    int m_hessian = 0, m_num_electrons = 0, m_mo_homo = -1, m_mo_lumo = -1;
     int m_threads = 1;
     int m_ConvCount = 1;
-    double m_dE = 0.1, m_dRMSD = 0.01, m_maxenergy = 100, m_GradNorm = 1e-5, m_lambda = 0.1;
+    double m_dE = 0.1, m_dRMSD = 0.01, m_maxenergy = 100, m_GradNorm = 1e-5, m_lambda = 0.1, m_mo_scale = 1.0;
     int m_charge = 0, m_spin = 0;
     int m_serial = false;
     int m_maxiter = 100, m_maxrise = 10, m_optimethod = 1, m_diis_hist = 10, m_diis_start = 10;

src/core/eht.cpp

@@ -31,14 +31,29 @@
 
 #include "eht.h"
 
+#include <iomanip>
+
 EHT::EHT()
 {
 }
 
-void EHT::start()
+void EHT::CalculateEHT(bool gradient, bool verbose)
 {
-    m_molecule.print_geom();
+    m_verbose = verbose;
+    if (gradient) {
+        std::cout << "EHT does not support gradients." << std::endl;
+    }
 
+    if (m_molecule.AtomCount() == 0) {
+        std::cout << "No molecule set." << std::endl;
+        return;
+    }
+
+    start();
+}
+
+std::vector<STO_6G> EHT::MakeBasis()
+{
     std::vector<STO_6G> basisset;
     for (int i = 0; i < m_molecule.Atoms().size(); ++i) {
         if (m_molecule.Atom(i).first == 1) {
@@ -148,14 +163,23 @@ void EHT::start()
             std::cout << "O" << std::endl;
         }
     }
-    std::cout << basisset.size() << std::endl;
+    return basisset;
+}
+
+void EHT::start()
+{
+    m_molecule.print_geom();
+    auto basisset = MakeBasis();
+    if (m_verbose)
+        std::cout << basisset.size() << std::endl;
     Matrix S = MakeOverlap(basisset);
     // std::cout << S << std::endl;
     Matrix H = MakeH(S, basisset);
-    std::cout << std::endl
-              << std::endl;
-    // std::cout << H << std::endl;
-
+    if (m_verbose) {
+        std::cout << std::endl
+                  << std::endl;
+        // std::cout << H << std::endl;
+    }
     Matrix S_1_2 = Matrix::Zero(basisset.size(), basisset.size());
     Eigen::JacobiSVD<Matrix> svd(S, Eigen::ComputeThinU | Eigen::ComputeThinV);
 
@@ -171,14 +195,17 @@ void EHT::start()
     Matrix F = S_1_2.transpose() * H * S_1_2;
     diag_F.compute(F);
     std::cout << std::endl;
-    // std::cout << diag_F.eigenvalues() <<std::endl;
+    m_energies = diag_F.eigenvalues();
+    m_mo = diag_F.eigenvectors();
 
     double energy = 0;
     for (int i = 0; i < m_num_electrons / 2; ++i) {
         energy += diag_F.eigenvalues()(i) * 2;
-        std::cout << diag_F.eigenvalues()(i) << " 2" << std::endl;
+        if (m_verbose)
+            std::cout << diag_F.eigenvalues()(i) << " 2" << std::endl;
     }
     std::cout << "Total electronic energy = " << energy << " Eh." << std::endl;
+
     // std::cout << diag_F.eigenvalues().sum();
 }
 

src/core/eht.h

@@ -80,7 +80,14 @@ class EHT {
     void setMolecule(const Molecule& molecule) { m_molecule = molecule; }
     void start();
 
+    Matrix MolecularOrbitals() const { return m_mo; }
+    Vector Energies() const { return m_energies; }
+    void CalculateEHT(bool gradient, bool verbose = false);
+    int NumElectrons() const { return m_num_electrons; }
+
 private:
+    std::vector<STO_6G> MakeBasis();
+
     Molecule m_molecule;
     /* Some integrals */
     /* s - s - sigma bond */
@@ -97,5 +104,11 @@ class EHT {
 
     Matrix MakeOverlap(const std::vector<STO_6G>& basisset);
     Matrix MakeH(const Matrix& S, const std::vector<STO_6G>& basisset);
+
     int m_num_electrons = 0;
+
+    Matrix m_mo;
+    Vector m_energies;
+
+    bool m_verbose = false;
 };

src/core/energycalculator.cpp

@@ -149,6 +149,13 @@ EnergyCalculator::EnergyCalculator(const std::string& method, const json& contro
             this->CalculateFF(gradient, verbose);
         };
 
+    } else if (m_method.compare("eht") == 0) {
+        m_eht = new EHT();
+        m_ecengine = [this](bool gradient, bool verbose) {
+            this->m_eht->CalculateEHT(gradient, verbose);
+            m_orbital_energies = this->m_eht->Energies();
+            m_num_electrons = this->m_eht->NumElectrons();
+        };
     } else { // Fall back to UFF?
         m_uff = new eigenUFF(controller);
     }
@@ -251,7 +258,11 @@ void EnergyCalculator::setMolecule(const Molecule& molecule)
 
         m_forcefield->setParameter(m_parameter);
 
-    } else { // Fall back to UFF?
+    } else if (m_method.compare("eht") == 0) {
+        m_eht->setMolecule(molecule);
+        // m_eht->Initialise();
+    } else {
+        // Fall back to UFF?
     }
     m_initialised = true;
 }
@@ -352,6 +363,10 @@ void EnergyCalculator::CalculateTBlite(bool gradient, bool verbose)
         }
     } else
         m_energy = m_tblite->GFNCalculation(m_gfn);
+
+    m_orbital_energies = m_tblite->OrbitalEnergies();
+    m_orbital_occupation = m_tblite->OrbitalOccupations();
+    m_num_electrons = m_orbital_occupation.sum();
 #endif
 }
 

src/core/energycalculator.h

@@ -37,6 +37,7 @@
 #include "src/core/dftd4interface.h"
 #endif
 
+#include "src/core/eht.h"
 #include "src/core/eigen_uff.h"
 #include "src/core/forcefield.h"
 #include "src/core/qmdff.h"
@@ -115,10 +116,14 @@ class EnergyCalculator {
     }
 
     inline json Parameter() const { return m_parameter; }
+    Vector Energies() const { return m_orbital_energies; }
+    Vector OrbitalOccuptations() const { return m_orbital_occupation; }
+
     std::vector<double> Charges() const;
     Position Dipole() const;
 
     std::vector<std::vector<double>> BondOrders() const;
+    int NumElectrons() const { return m_num_electrons; }
 
 private:
     void InitialiseUFF();
@@ -161,24 +166,28 @@ class EnergyCalculator {
     eigenUFF* m_uff = NULL;
     QMDFF* m_qmdff = NULL;
     ForceField* m_forcefield = NULL;
+    EHT* m_eht = NULL;
     StringList m_uff_methods = { "fuff" };
     StringList m_ff_methods = { "uff", "uff-d3", "qmdff" };
     StringList m_qmdff_method = { "fqmdff" };
     StringList m_tblite_methods = { "ipea1", "gfn1", "gfn2" };
     StringList m_xtb_methods = { "gfnff", "xtb-gfn1", "xtb-gfn2" };
     StringList m_d3_methods = { "d3" };
     StringList m_d4_methods = { "d4" };
+
     std::function<void(bool, bool)> m_ecengine;
     std::function<std::vector<double>()> m_charges;
     std::function<Position()> m_dipole;
     std::function<std::vector<std::vector<double>>()> m_bonds;
     json m_parameter;
     std::string m_method, m_param_file;
-    Matrix m_geometry, m_gradient;
+    Matrix m_geometry, m_gradient, m_molecular_orbitals;
+    Vector m_orbital_energies, m_orbital_occupation;
+
     double m_energy;
     double *m_coord, *m_grad;
 
-    int m_atoms;
+    int m_atoms, m_num_electrons = 0;
     int m_gfn = 2;
     int* m_atom_type;
 

src/core/tbliteinterface.cpp

@@ -352,6 +352,38 @@ std::vector<std::vector<double>> TBLiteInterface::BondOrders() const
     return bond_orders;
 }
 
+Vector TBLiteInterface::OrbitalEnergies() const
+{
+    int num_orbitals;
+    tblite_get_result_number_of_orbitals(m_error, m_tblite_res, &num_orbitals);
+
+    Vector orbital_energies(num_orbitals);
+    std::vector<double> energies(num_orbitals);
+    tblite_get_result_orbital_energies(m_error, m_tblite_res, energies.data());
+
+    for (int i = 0; i < num_orbitals; ++i) {
+        orbital_energies[i] = energies[i];
+    }
+
+    return orbital_energies;
+}
+
+Vector TBLiteInterface::OrbitalOccupations() const
+{
+    int num_orbitals;
+    tblite_get_result_number_of_orbitals(m_error, m_tblite_res, &num_orbitals);
+
+    Vector orbital_occupations(num_orbitals);
+    std::vector<double> occupations(num_orbitals);
+    tblite_get_result_orbital_occupations(m_error, m_tblite_res, occupations.data());
+
+    for (int i = 0; i < num_orbitals; ++i) {
+        orbital_occupations[i] = occupations[i];
+    }
+
+    return orbital_occupations;
+}
+
 void TBLiteInterface::tbliteError()
 {
     char message[512];