load mol2 files, fix opt before md

Signed-off-by: Conrad Hübler <Conrad.Huebler@gmx.net>

src/capabilities/curcumaopt.cpp

@@ -90,8 +90,8 @@ void CurcumaOpt::LoadControlJson()
 
 void CurcumaOpt::start()
 {
-    getBasename(m_filename);
     if (m_file_set) {
+        getBasename(m_filename);
         FileIterator file(m_filename);
         std::multimap<double, Molecule> results;
         while (!file.AtEnd()) {
@@ -246,19 +246,21 @@ double CurcumaOpt::SinglePoint(const Molecule* initial, const json& controller,
 Molecule CurcumaOpt::LBFGSOptimise(Molecule* initial, const json& controller, std::string& output, std::vector<Molecule>* intermediate, int thread, const std::string& basename)
 {
     bool printOutput = Json2KeyWord<bool>(controller, "printOutput");
+    bool fusion = Json2KeyWord<bool>(controller, "fusion");
     double dE = Json2KeyWord<double>(controller, "dE");
     double dRMSD = Json2KeyWord<double>(controller, "dRMSD");
     double GradNorm = Json2KeyWord<double>(controller, "GradNorm");
-
+    double maxrise = Json2KeyWord<double>(controller, "maxrise");
     std::string method = Json2KeyWord<std::string>(controller, "method");
     int MaxIter = Json2KeyWord<int>(controller, "MaxIter");
     int ConvCount = Json2KeyWord<int>(controller, "ConvCount");
     int SingleStep = Json2KeyWord<int>(controller, "SingleStep");
-
-    if(Json2KeyWord<int>(controller, "threads") > 1 )
-    {
-        printOutput = false;
-    }
+    /*
+        if(Json2KeyWord<int>(controller, "threads") > 1 )
+        {
+            printOutput = false;
+        }
+        */
     bool optH = Json2KeyWord<bool>(controller, "optH");
     std::vector<int> constrain;
     Geometry geometry = initial->getGeometry();
@@ -276,6 +278,7 @@ Molecule CurcumaOpt::LBFGSOptimise(Molecule* initial, const json& controller, st
     }
 
     EnergyCalculator interface(method, controller);
+
     interface.setMolecule(*initial);
 
     double final_energy = interface.CalculateEnergy(true);
@@ -284,7 +287,7 @@ Molecule CurcumaOpt::LBFGSOptimise(Molecule* initial, const json& controller, st
     std::cout << "Initial energy " << final_energy << "Eh" << std::endl;
     LBFGSParam<double> param;
     param.m = Json2KeyWord<int>(controller, "LBFGS_m");
-    ;
+
     param.epsilon = Json2KeyWord<double>(controller, "LBFGS_eps_abs");
     param.epsilon_rel = Json2KeyWord<double>(controller, "LBFGS_eps_rel");
     param.past = Json2KeyWord<int>(controller, "LBFGS_past");
@@ -388,7 +391,7 @@ Molecule CurcumaOpt::LBFGSOptimise(Molecule* initial, const json& controller, st
             }
             next.setGeometry(geometry);
         }
-        if ((fun.m_energy - final_energy) * 2625.5 > 10 && iteration > 10) {
+        if ((fun.m_energy - final_energy) * 2625.5 > maxrise && iteration > 10) {
             if (printOutput) {
                 output += fmt::format("Energy rises too much!\n");
                 output += fmt::format("{0: ^75}\n\n", "*** Geometry Optimisation sufficiantly converged ***");
@@ -464,13 +467,15 @@ Molecule CurcumaOpt::LBFGSOptimise(Molecule* initial, const json& controller, st
         }
         */
         final_energy = fun.m_energy;
-        if (next.Check() == 0) {
+        if (next.Check() == 0 || (next.Check() == 1 && fusion)) {
             previous = next;
             next.setEnergy(final_energy);
             intermediate->push_back(next);
             next.appendXYZFile(basename + ".t" + std::to_string(thread) + ".xyz");
 
         } else {
+            output += fmt::format("{0: ^75}\n\n", "*** Check next failed! ***");
+
             perform_optimisation = false;
             error = true;
         }

src/capabilities/curcumaopt.h

@@ -52,7 +52,9 @@ static json CurcumaOptJson{
     { "SinglePoint", false },
     { "optH", false },
     { "serial", false },
-    { "hessian", 0 }
+    { "hessian", 0 },
+    { "fusion", false },
+    { "maxrise", 100 }
 };
 
 const json OptJsonPrivate{
@@ -174,10 +176,10 @@ class CurcumaOpt : public CurcumaMethod {
     std::string m_method = "UFF";
     Molecule m_molecule;
     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;
+    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;
     int m_hessian = 0;
     int m_threads = 1;
-    double m_dE = 0.1, m_dRMSD = 0.01;
+    double m_dE = 0.1, m_dRMSD = 0.01, m_maxenergy = 100;
     int m_charge = 0, m_spin = 0;
     int m_serial = false;
 };

src/capabilities/simplemd.cpp

@@ -239,9 +239,10 @@ bool SimpleMD::Initialise()
         optimise.addMolecule(&m_molecule);
         optimise.start();
         auto mol = optimise.Molecules();
-        std::cout << mol->size() << std::endl;
+
         auto molecule = ((*mol)[0]);
         m_molecule.setGeometry(molecule.getGeometry());
+        m_molecule.appendXYZFile(Basename() + ".opt.xyz");
     }
 
     for (int i = 0; i < m_natoms; ++i) {
@@ -337,7 +338,6 @@ void SimpleMD::InitVelocities(double scaling)
     double Px = 0.0, Py = 0.0, Pz = 0.0;
     for (int i = 0; i < m_natoms; ++i) {
         double v0 = sqrt(kb_Eh * m_T0 * amu2au / (m_mass[i])) * scaling / fs2amu;
-        // double v0 = sqrt(kb_SI * m_T0 / (m_mass[i] * atomic_mass)) * scaling;
         m_velocities[3 * i + 0] = v0 * d(gen);
         m_velocities[3 * i + 1] = v0 * d(gen);
         m_velocities[3 * i + 2] = v0 * d(gen);

src/core/eigen_uff.cpp

@@ -84,7 +84,7 @@ double UFFThread::Distance(double x1, double x2, double y1, double y2, double z1
 
 double UFFThread::BondEnergy(double distance, double r, double kij, double D_ij)
 {
-    double energy = (0.5 * kij * (distance - r) * (distance - r)) * m_final_factor * m_bond_scaling;
+    double energy = (0.25 * kij * (distance - r) * (distance - r)) * m_final_factor * m_bond_scaling;
     if (isnan(energy))
         return 0;
     else
@@ -106,7 +106,7 @@ double UFFThread::CalculateBondStretching()
         Matrix derivate;
         double r0 = BondStretching(x, y, derivate, m_CalculateGradient);
 
-        energy += (0.5 * bond.kij * (r0 - bond.r0) * (r0 - bond.r0)) * m_final_factor * m_bond_scaling;
+        energy += (0.25 * bond.kij * (r0 - bond.r0) * (r0 - bond.r0)) * m_final_factor * m_bond_scaling;
         if (m_CalculateGradient) {
             if (m_calc_gradient == 0) {
                 double diff = (bond.kij) * (r0 - bond.r0) * m_final_factor * m_bond_scaling;
@@ -560,6 +560,7 @@ double UFFThread::CalculateElectrostatic()
 eigenUFF::eigenUFF(const json& controller)
 {
     json parameter = MergeJson(UFFParameterJson, controller);
+    std::cout << parameter["threads"] << std::endl;
     m_threadpool = new CxxThreadPool();
     m_threadpool->setProgressBar(CxxThreadPool::ProgressBarType::None);
 #ifdef USE_D3
@@ -605,7 +606,7 @@ eigenUFF::~eigenUFF()
         delete m_stored_threads[i];
 }
 
-void eigenUFF::Initialise()
+void eigenUFF::Initialise(const std::vector<std::pair<int, int>>& formed_bonds)
 {
     if (m_initialised)
         return;
@@ -617,33 +618,58 @@ void eigenUFF::Initialise()
     TContainer bonds, nonbonds, angles, dihedrals, inversions;
     m_scaling = 1.4;
     m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3);
-    for (int i = 0; i < m_atom_types.size(); ++i) {
-        m_stored_bonds.push_back(std::vector<int>());
-        ignored_vdw.push_back(std::set<int>({ i }));
-        for (int j = 0; j < m_atom_types.size() && m_stored_bonds[i].size() < CoordinationNumber[m_atom_types[i]]; ++j) {
-            if (i == j)
-                continue;
-            double x_i = m_geometry(i, 0) * m_au;
-            double x_j = m_geometry(j, 0) * m_au;
+    if (formed_bonds.size() == 0) {
+        for (int i = 0; i < m_atom_types.size(); ++i) {
+            m_stored_bonds.push_back(std::vector<int>());
+            ignored_vdw.push_back(std::set<int>({ i }));
+            for (int j = 0; j < m_atom_types.size() && m_stored_bonds[i].size() < CoordinationNumber[m_atom_types[i]]; ++j) {
+                if (i == j)
+                    continue;
+                double x_i = m_geometry(i, 0) * m_au;
+                double x_j = m_geometry(j, 0) * m_au;
 
-            double y_i = m_geometry(i, 1) * m_au;
-            double y_j = m_geometry(j, 1) * m_au;
+                double y_i = m_geometry(i, 1) * m_au;
+                double y_j = m_geometry(j, 1) * m_au;
 
-            double z_i = m_geometry(i, 2) * m_au;
-            double z_j = m_geometry(j, 2) * m_au;
+                double z_i = m_geometry(i, 2) * m_au;
+                double z_j = m_geometry(j, 2) * m_au;
 
-            double r_ij = sqrt((((x_i - x_j) * (x_i - x_j)) + ((y_i - y_j) * (y_i - y_j)) + ((z_i - z_j) * (z_i - z_j))));
+                double r_ij = sqrt((((x_i - x_j) * (x_i - x_j)) + ((y_i - y_j) * (y_i - y_j)) + ((z_i - z_j) * (z_i - z_j))));
 
-            if (r_ij <= (Elements::CovalentRadius[m_atom_types[i]] + Elements::CovalentRadius[m_atom_types[j]]) * m_scaling * m_au) {
-                if (bonds.insert({ std::min(i, j), std::max(i, j) })) {
-                    m_coordination[i]++;
-                    m_stored_bonds[i].push_back(j);
-                    ignored_vdw[i].insert(j);
+                if (r_ij <= (Elements::CovalentRadius[m_atom_types[i]] + Elements::CovalentRadius[m_atom_types[j]]) * m_scaling * m_au) {
+                    if (bonds.insert({ std::min(i, j), std::max(i, j) })) {
+                        m_coordination[i]++;
+                        m_stored_bonds[i].push_back(j);
+                        ignored_vdw[i].insert(j);
+                    }
+                    m_topo(i, j) = 1;
+                    m_topo(j, i) = 1;
                 }
-                m_topo(i, j) = 1;
-                m_topo(j, i) = 1;
             }
         }
+    } else {
+        for (int i = 0; i < m_atom_types.size(); ++i) {
+            m_stored_bonds.push_back(std::vector<int>());
+            ignored_vdw.push_back(std::set<int>({ i }));
+        }
+        for (const std::pair<int, int>& bond : formed_bonds) {
+
+            int i = bond.first - 1;
+            int j = bond.second - 1;
+
+            if (bonds.insert({ std::min(i, j), std::max(i, j) })) {
+                m_coordination[i]++;
+                m_coordination[j]++;
+
+                m_stored_bonds[i].push_back(j);
+                m_stored_bonds[j].push_back(i);
+
+                ignored_vdw[i].insert(j);
+                ignored_vdw[j].insert(i);
+            }
+            m_topo(i, j) = 1;
+            m_topo(j, i) = 1;
+        }
     }
     AssignUffAtomTypes();
     if (m_rings)
@@ -1723,13 +1749,14 @@ double eigenUFF::Calculate(bool grd, bool verbose)
     double dihedral_energy = 0.0;
     double inversion_energy = 0.0;
     double vdw_energy = 0.0;
-    m_threadpool->setActiveThreadCount(m_threads);
 
     for (int i = 0; i < m_stored_threads.size(); ++i) {
         m_stored_threads[i]->UpdateGeometry(&m_geometry);
     }
 
     m_threadpool->Reset();
+    m_threadpool->setActiveThreadCount(m_threads);
+
     m_threadpool->StartAndWait();
     m_threadpool->setWakeUp(m_threadpool->WakeUp() / 2.0);
     for (int i = 0; i < m_stored_threads.size(); ++i) {

src/core/eigen_uff.h

@@ -226,7 +226,7 @@ class eigenUFF {
         }
     }
 
-    void Initialise();
+    void Initialise(const std::vector<std::pair<int, int>>& formed_bonds = std::vector<std::pair<int, int>>());
 
     double Calculate(bool gradient = true, bool verbose = false);
 
@@ -250,6 +250,7 @@ class eigenUFF {
 
     void setBonds(const json& bonds);
     void setBonds(const TContainer& bonds, std::vector<std::set<int>>& ignored_vdw, TContainer& angels, TContainer& dihedrals, TContainer& inversions);
+    void setBonds(const std::vector<std::pair<int, int>>& bonds);
 
     void setAngles(const json& angles);
     void setAngles(const TContainer& angles, const std::vector<std::set<int>>& ignored_vdw);

src/core/energycalculator.cpp

@@ -165,7 +165,7 @@ void EnergyCalculator::setMolecule(const Molecule& molecule)
     m_geometry = geom;
     if (std::find(m_uff_methods.begin(), m_uff_methods.end(), m_method) != m_uff_methods.end()) { // UFF energy calculator requested
         m_uff->setMolecule(atoms, geom);
-        m_uff->Initialise();
+        m_uff->Initialise(molecule.Bonds());
     } else if (std::find(m_tblite_methods.begin(), m_tblite_methods.end(), m_method) != m_tblite_methods.end()) { // TBLite energy calculator requested
 #ifdef USE_TBLITE
         m_tblite->InitialiseMolecule(molecule);

src/core/molecule.cpp

@@ -57,6 +57,7 @@ Molecule::Molecule(const Molecule& other)
     m_name = other.m_name;
     m_energy = other.m_energy;
     m_spin = other.m_spin;
+    m_bonds = other.m_bonds;
 }
 /*
 Molecule& Molecule::operator=(const Molecule& other)
@@ -80,6 +81,7 @@ Molecule::Molecule(const Molecule* other)
     m_name = other->m_name;
     m_energy = other->m_energy;
     m_spin = other->m_spin;
+    m_bonds = other->m_bonds;
 }
 /*
 Molecule& Molecule::operator=(const Molecule* other)
@@ -102,6 +104,7 @@ Molecule::Molecule(const Mol& other)
     m_atoms = other.m_atoms;
     m_energy = other.m_energy;
     m_spin = other.m_spin;
+    m_bonds = other.m_bonds;
 }
 
 Molecule::Molecule(const Mol* other)
@@ -112,6 +115,7 @@ Molecule::Molecule(const Mol* other)
     m_atoms = other->m_atoms;
     m_energy = other->m_energy;
     m_spin = other->m_spin;
+    m_bonds = other->m_bonds;
 }
 
 Molecule::Molecule(const std::string& file)
@@ -239,9 +243,11 @@ void Molecule::print_geom(bool moreinfo) const
 
 int Molecule::Check() const
 {
-    for (int i = 1; i < AtomCount(); ++i) {
+    return 0;
+
+    for (int i = 0; i < AtomCount(); ++i) {
         if (std::isnan(m_geometry[i][0]) || std::isnan(m_geometry[i][1]) || std::isnan(m_geometry[i][2]))
-            return 1;
+            return 2;
         for (int j = 0; j < i; ++j) {
             if (CalculateDistance(i, j) < 1e-1)
                 return 1;
@@ -619,6 +625,7 @@ void Molecule::LoadMolecule(const Molecule& molecule)
     m_charge = molecule.Charge();
     m_atoms = molecule.Atoms();
     m_energy = molecule.Energy();
+    m_bonds = molecule.m_bonds;
     InitialiseEmptyGeometry(molecule.AtomCount());
     setGeometry(molecule.getGeometry());
 }
@@ -628,6 +635,8 @@ void Molecule::LoadMolecule(const Molecule* molecule)
     clear();
     m_charge = molecule->Charge();
     m_atoms = molecule->Atoms();
+    m_bonds = molecule->m_bonds;
+
     InitialiseEmptyGeometry(molecule->AtomCount());
     setGeometry(molecule->getGeometry());
 }
@@ -642,6 +651,7 @@ void Molecule::LoadMolecule(const Mol& molecule)
     m_energy = molecule.m_energy;
     InitialiseEmptyGeometry(AtomCount());
     m_geometry = (molecule.m_geometry);
+    m_bonds = molecule.m_bonds;
 }
 
 void Molecule::LoadMolecule(const Mol* molecule)
@@ -654,6 +664,7 @@ void Molecule::LoadMolecule(const Mol* molecule)
     m_energy = molecule->m_energy;
     InitialiseEmptyGeometry(AtomCount());
     m_geometry = (molecule->m_geometry);
+    m_bonds = molecule->m_bonds;
 }
 
 Molecule Molecule::getFragmentMolecule(int fragment) const

src/core/molecule.h

@@ -44,6 +44,7 @@ struct Mol {
     std::string m_commentline;
 
     std::vector<std::array<double, 3>> m_geometry;
+    std::vector<std::pair<int, int>> m_bonds;
     std::vector<int> m_atoms;
 };
 
@@ -213,6 +214,8 @@ class Molecule
 
     void setPartialCharges(const std::vector<double>& charges) { m_charges = charges; }
 
+    inline std::vector<std::pair<int, int>> Bonds() const { return m_bonds; }
+
 private:
     void ParseString(const std::string& internal, std::vector<std::string>& elements);
 
@@ -243,6 +246,8 @@ class Molecule
     mutable std::map<int, int> m_fragment_assignment;
 
     mutable std::vector<double> m_mass_fragments;
+    std::vector<std::pair<int, int>> m_bonds;
+
     mutable bool m_dirty = true;
     std::string m_name;
     double m_energy = 0, m_Ia = 0, m_Ib = 0, m_Ic = 0, m_mass = 0, m_hbond_cutoff = 3;