Merge branch 'rework_executable' of github.com:peterspackman/occ into…

… rework_executable
peterspackman · Nov 20, 2023 · b350a10 · b350a10
2 parents d0d4016 + fb9e5f5
commit b350a10
Show file tree

Hide file tree

Showing 21 changed files with 150 additions and 127 deletions.
diff --git a/include/occ/core/dimer.h b/include/occ/core/dimer.h
@@ -287,6 +287,7 @@ class Dimer {
 
     inline const auto &name() const { return m_name; }
     inline void set_name(const std::string &name) { m_name = name; }
+    std::string xyz_string() const;
 
   private:
     Molecule m_a, m_b;

diff --git a/include/occ/dft/dft.h b/include/occ/dft/dft.h
@@ -103,10 +103,15 @@ class DFT {
         return m_hf.frozen_electrons();
     }
 
-    void set_density_fitting_basis(const std::string &density_fitting_basis) {
+    inline void
+    set_density_fitting_basis(const std::string &density_fitting_basis) {
         m_hf.set_density_fitting_basis(density_fitting_basis);
     }
 
+    inline void set_precision(double precision) {
+        m_hf.set_precision(precision);
+    }
+
     double exchange_correlation_energy() const { return m_exc_dft; }
     double exchange_energy_total() const { return m_exchange_energy; }
 
@@ -203,8 +208,7 @@ class DFT {
 
     template <int derivative_order,
               SpinorbitalKind spinorbital_kind = SpinorbitalKind::Restricted>
-    Mat compute_K_dft(const MolecularOrbitals &mo, double precision,
-                      const Mat &Schwarz) {
+    Mat compute_K_dft(const MolecularOrbitals &mo, const Mat &Schwarz) {
         using occ::parallel::nthreads;
         const auto &basis = m_hf.aobasis();
         const auto &atoms = m_hf.atoms();
@@ -322,31 +326,27 @@ class DFT {
         return K;
     }
 
-    inline Mat
-    compute_J(const MolecularOrbitals &mo,
-              double precision = std::numeric_limits<double>::epsilon(),
-              const Mat &Schwarz = Mat()) const {
-        return m_hf.compute_J(mo, precision, Schwarz);
+    inline Mat compute_J(const MolecularOrbitals &mo,
+                         const Mat &Schwarz = Mat()) const {
+        return m_hf.compute_J(mo, Schwarz);
     }
 
-    inline Mat
-    compute_vxc(const MolecularOrbitals &mo,
-                double precision = std::numeric_limits<double>::epsilon(),
-                const Mat &Schwarz = Mat()) {
+    inline Mat compute_vxc(const MolecularOrbitals &mo,
+                           const Mat &Schwarz = Mat()) {
         int deriv = density_derivative();
         switch (mo.kind) {
         case SpinorbitalKind::Unrestricted: {
             occ::log::debug("Unrestricted vxc evaluation");
             switch (deriv) {
             case 0:
-                return compute_K_dft<0, SpinorbitalKind::Unrestricted>(
-                    mo, precision, Schwarz);
+                return compute_K_dft<0, SpinorbitalKind::Unrestricted>(mo,
+                                                                       Schwarz);
             case 1:
-                return compute_K_dft<1, SpinorbitalKind::Unrestricted>(
-                    mo, precision, Schwarz);
+                return compute_K_dft<1, SpinorbitalKind::Unrestricted>(mo,
+                                                                       Schwarz);
             case 2:
-                return compute_K_dft<2, SpinorbitalKind::Unrestricted>(
-                    mo, precision, Schwarz);
+                return compute_K_dft<2, SpinorbitalKind::Unrestricted>(mo,
+                                                                       Schwarz);
             default:
                 throw std::runtime_error(
                     "Not implemented: DFT for derivative order > 2");
@@ -355,14 +355,14 @@ class DFT {
         case SpinorbitalKind::Restricted: {
             switch (deriv) {
             case 0:
-                return compute_K_dft<0, SpinorbitalKind::Restricted>(
-                    mo, precision, Schwarz);
+                return compute_K_dft<0, SpinorbitalKind::Restricted>(mo,
+                                                                     Schwarz);
             case 1:
-                return compute_K_dft<1, SpinorbitalKind::Restricted>(
-                    mo, precision, Schwarz);
+                return compute_K_dft<1, SpinorbitalKind::Restricted>(mo,
+                                                                     Schwarz);
             case 2:
-                return compute_K_dft<2, SpinorbitalKind::Restricted>(
-                    mo, precision, Schwarz);
+                return compute_K_dft<2, SpinorbitalKind::Restricted>(mo,
+                                                                     Schwarz);
             default:
                 throw std::runtime_error(
                     "Not implemented: DFT for derivative order > 2");
@@ -374,23 +374,21 @@ class DFT {
         }
     }
 
-    inline std::pair<Mat, Mat>
-    compute_JK(const MolecularOrbitals &mo,
-               double precision = std::numeric_limits<double>::epsilon(),
-               const Mat &Schwarz = Mat()) {
+    inline std::pair<Mat, Mat> compute_JK(const MolecularOrbitals &mo,
+                                          const Mat &Schwarz = Mat()) {
 
         Mat J;
-        Mat K = -compute_vxc(mo, precision, Schwarz);
+        Mat K = -compute_vxc(mo, Schwarz);
         double ecoul{0.0}, exc{0.0};
         double exchange_factor = exact_exchange_factor();
         RangeSeparatedParameters rs = range_separated_parameters();
         if (rs.omega != 0.0) {
             // range separated hybrid
             Mat Ksr;
-            std::tie(J, Ksr) = m_hf.compute_JK(mo, precision, Schwarz);
+            std::tie(J, Ksr) = m_hf.compute_JK(mo, Schwarz);
             m_hf.set_range_separated_omega(rs.omega);
             Mat Jlr, Klr;
-            std::tie(Jlr, Klr) = m_hf.compute_JK(mo, precision, Schwarz);
+            std::tie(Jlr, Klr) = m_hf.compute_JK(mo, Schwarz);
             Mat Khf = Ksr * (rs.alpha + rs.beta) + Klr * (-rs.beta);
             ecoul = expectation(mo.kind, mo.D, J);
             exc = -expectation(mo.kind, mo.D, Khf);
@@ -399,12 +397,12 @@ class DFT {
         } else if (exchange_factor != 0.0) {
             // global hybrid
             Mat Khf;
-            std::tie(J, Khf) = m_hf.compute_JK(mo, precision, Schwarz);
+            std::tie(J, Khf) = m_hf.compute_JK(mo, Schwarz);
             ecoul = expectation(mo.kind, mo.D, J);
             exc = -expectation(mo.kind, mo.D, Khf) * exchange_factor;
             K.noalias() += Khf * exchange_factor;
         } else {
-            J = m_hf.compute_J(mo, precision, Schwarz);
+            J = m_hf.compute_J(mo, Schwarz);
             ecoul = expectation(mo.kind, mo.D, J);
         }
         occ::log::debug(
@@ -433,10 +431,8 @@ class DFT {
         return m_nlc_energy;
     }
 
-    Mat compute_fock(const MolecularOrbitals &mo,
-                     double precision = std::numeric_limits<double>::epsilon(),
-                     const Mat &Schwarz = Mat()) {
-        auto [J, K] = compute_JK(mo, precision, Schwarz);
+    Mat compute_fock(const MolecularOrbitals &mo, const Mat &Schwarz = Mat()) {
+        auto [J, K] = compute_JK(mo, Schwarz);
         return J - K;
     }
 

diff --git a/include/occ/interaction/pairinteraction.h b/include/occ/interaction/pairinteraction.h
@@ -60,7 +60,6 @@ inline CEParameterizedModel CE2_XDM = CE2_XDM_WB97MV;
 inline CEParameterizedModel CE5_XDM = CE5_XDM_WB97MV;
 
 struct CEMonomerCalculationParameters {
-    double precision{1e-12};
     Mat Schwarz;
     bool xdm{false};
     bool neglect_exchange{false};

diff --git a/include/occ/io/occ_input.h b/include/occ/io/occ_input.h
@@ -38,6 +38,7 @@ struct DriverInput {
 struct MethodInput {
     std::string name{"rhf"};
     BeckeGridSettings dft_grid;
+    double integral_precision{1e-12};
 };
 
 struct BasisSetInput {

diff --git a/include/occ/qm/hf.h b/include/occ/qm/hf.h
@@ -37,23 +37,30 @@ class HartreeFock {
     inline bool have_effective_core_potentials() const {
         return m_engine.have_effective_core_potentials();
     }
+
     void set_system_charge(int charge);
     void set_density_fitting_basis(const std::string &);
+
+    inline void set_precision(double precision) {
+        m_engine.set_precision(precision);
+
+        if (m_df_engine != nullptr) {
+            m_df_engine->set_precision(precision);
+        }
+    }
+
     double nuclear_repulsion_energy() const;
-    double nuclear_point_charge_interaction_energy(const PointChargeList &) const;
+    double
+    nuclear_point_charge_interaction_energy(const PointChargeList &) const;
 
     Mat compute_fock(const MolecularOrbitals &mo,
-                     double precision = std::numeric_limits<double>::epsilon(),
                      const Mat &Schwarz = Mat()) const;
 
     Mat compute_fock_mixed_basis(const MolecularOrbitals &mo_minbs,
                                  const qm::AOBasis &bs, bool is_shell_diagonal);
-    std::pair<Mat, Mat>
-    compute_JK(const MolecularOrbitals &mo,
-               double precision = std::numeric_limits<double>::epsilon(),
-               const Mat &Schwarz = Mat()) const;
+    std::pair<Mat, Mat> compute_JK(const MolecularOrbitals &mo,
+                                   const Mat &Schwarz = Mat()) const;
     Mat compute_J(const MolecularOrbitals &mo,
-                  double precision = std::numeric_limits<double>::epsilon(),
                   const Mat &Schwarz = Mat()) const;
 
     Mat compute_kinetic_matrix() const;

diff --git a/include/occ/qm/integral_engine.h b/include/occ/qm/integral_engine.h
@@ -197,6 +197,8 @@ class IntegralEngine {
         m_env.set_range_separated_omega(omega);
     }
 
+    inline void set_precision(double precision) { m_precision = precision; }
+
   private:
     double m_precision{1e-12};
     AOBasis m_aobasis, m_auxbasis;

diff --git a/include/occ/qm/integral_engine_df.h b/include/occ/qm/integral_engine_df.h
@@ -26,6 +26,7 @@ class IntegralEngineDF {
     inline void set_integral_policy(Policy p) { m_policy = p; }
 
     void set_range_separated_omega(double omega);
+    void set_precision(double precision);
 
   private:
     inline size_t num_rows() const {
@@ -56,6 +57,8 @@ class IntegralEngineDF {
         return (m_policy == Policy::Stored);
     }
 
+    double m_precision{1e-12};
+
     mutable IntegralEngine m_ao_engine;  // engine with ao basis & aux basis
     mutable IntegralEngine m_aux_engine; // engine with just aux basis
     Eigen::LLT<Mat> V_LLt;

diff --git a/include/occ/qm/scf.h b/include/occ/qm/scf.h
@@ -490,15 +490,8 @@ template <typename Procedure> struct SCF {
             }
 
             // build a new Fock matrix
-            // totally empirical precision variation, involves the condition
-            // number
-            const auto precision_F =
-                std::min(std::min(1e-3 / XtX_condition_number, 1e-7),
-                         std::max(diis_error / 1e4,
-                                  std::numeric_limits<double>::epsilon()));
-
             std::swap(mo.D, D_diff);
-            F += m_procedure.compute_fock(mo, precision_F, K);
+            F += m_procedure.compute_fock(mo, K);
             std::swap(mo.D, D_diff);
 
             // compute HF energy with the non-extrapolated Fock matrix

diff --git a/include/occ/solvent/solvation_correction.h b/include/occ/solvent/solvation_correction.h
@@ -147,6 +147,11 @@ template <typename Proc> class SolvationCorrectedProcedure {
     inline Vec3 center_of_mass() const { return m_proc.center_of_mass(); }
 
     void set_system_charge(int charge) { m_proc.set_system_charge(charge); }
+
+    inline void set_precision(double precision) {
+        m_proc.set_precision(precision);
+    }
+
     int system_charge() const { return m_proc.system_charge(); }
     int total_electrons() const { return m_proc.total_electrons(); }
     int active_electrons() const { return m_proc.active_electrons(); }
@@ -160,7 +165,8 @@ template <typename Proc> class SolvationCorrectedProcedure {
         return m_proc.nuclear_repulsion_energy();
     }
 
-    inline double nuclear_point_charge_interaction_energy(const PointChargeList &pc) const {
+    inline double
+    nuclear_point_charge_interaction_energy(const PointChargeList &pc) const {
         return m_proc.nuclear_point_charge_interaction_energy(pc);
     }
 
@@ -201,7 +207,8 @@ template <typename Proc> class SolvationCorrectedProcedure {
 
     auto compute_schwarz_ints() const { return m_proc.compute_schwarz_ints(); }
 
-    inline auto compute_point_charge_interaction_matrix(const PointChargeList &pc) {
+    inline auto
+    compute_point_charge_interaction_matrix(const PointChargeList &pc) {
         return m_proc.compute_point_charge_interaction_matrix(pc);
     }
 
@@ -262,9 +269,8 @@ template <typename Proc> class SolvationCorrectedProcedure {
     }
 
     Mat compute_fock(const MolecularOrbitals &mo,
-                     double precision = std::numeric_limits<double>::epsilon(),
                      const Mat &Schwarz = Mat()) const {
-        return m_proc.compute_fock(mo, precision, Schwarz);
+        return m_proc.compute_fock(mo, Schwarz);
     }
 
     Mat compute_fock_mixed_basis(const MolecularOrbitals &mo_bs,

diff --git a/src/core/dimer.cpp b/src/core/dimer.cpp
@@ -202,4 +202,19 @@ bool Dimer::equivalent_under_rotation(const occ::core::Dimer &rhs,
     return occ::util::all_close(posd1_rot, posd2, 1e-5, 1e-5);
 }
 
+std::string Dimer::xyz_string() const {
+    using occ::core::Element;
+    std::string result;
+
+    const auto &pos = positions();
+    const auto &nums = atomic_numbers();
+    result += fmt::format("{}\n\n", nums.rows());
+    for (size_t i = 0; i < nums.rows(); i++) {
+        result += fmt::format("{:5s} {:12.5f} {:12.5f} {:12.5f}\n",
+                              Element(nums(i)).symbol(), pos(0, i), pos(1, i),
+                              pos(2, i));
+    }
+    return result;
+}
+
 } // namespace occ::core
diff --git a/src/crystal/crystal.cpp b/src/crystal/crystal.cpp
@@ -485,19 +485,19 @@ CrystalDimers Crystal::symmetry_unique_dimers(double radius) const {
 
     for (size_t i = 0; i < m_asymmetric_unit.size(); i++) {
         const auto &pos = m_asymmetric_unit.positions.col(i);
-        upper.h =
-            std::max(upper.h, static_cast<int>(ceil(pos(0) + frac_radius(0))));
-        upper.k =
-            std::max(upper.k, static_cast<int>(ceil(pos(1) + frac_radius(1))));
-        upper.l =
-            std::max(upper.l, static_cast<int>(ceil(pos(2) + frac_radius(2))));
-
-        lower.h =
-            std::min(lower.h, static_cast<int>(floor(pos(0) - frac_radius(0))));
-        lower.k =
-            std::min(lower.k, static_cast<int>(floor(pos(1) - frac_radius(1))));
-        lower.l =
-            std::min(lower.l, static_cast<int>(floor(pos(2) - frac_radius(2))));
+        upper.h = std::max(upper.h,
+                           static_cast<int>(ceil(pos(0) + frac_radius(0))) + 1);
+        upper.k = std::max(upper.k,
+                           static_cast<int>(ceil(pos(1) + frac_radius(1))) + 1);
+        upper.l = std::max(upper.l,
+                           static_cast<int>(ceil(pos(2) + frac_radius(2))) + 1);
+
+        lower.h = std::min(
+            lower.h, static_cast<int>(floor(pos(0) - frac_radius(0))) - 1);
+        lower.k = std::min(
+            lower.k, static_cast<int>(floor(pos(1) - frac_radius(1))) - 1);
+        lower.l = std::min(
+            lower.l, static_cast<int>(floor(pos(2) - frac_radius(2))) - 1);
     }
 
     const auto &uc_mols = unit_cell_molecules();
@@ -573,19 +573,19 @@ CrystalDimers Crystal::unit_cell_dimers(double radius) const {
         Mat3N pos_frac = to_fractional(mol.positions());
         for (size_t i = 0; i < pos_frac.cols(); i++) {
             const auto &pos = pos_frac.col(i);
-            upper.h = std::max(upper.h,
-                               static_cast<int>(ceil(pos(0) + frac_radius(0))));
-            upper.k = std::max(upper.k,
-                               static_cast<int>(ceil(pos(1) + frac_radius(1))));
-            upper.l = std::max(upper.l,
-                               static_cast<int>(ceil(pos(2) + frac_radius(2))));
+            upper.h = std::max(
+                upper.h, static_cast<int>(ceil(pos(0) + frac_radius(0))) + 1);
+            upper.k = std::max(
+                upper.k, static_cast<int>(ceil(pos(1) + frac_radius(1))) + 1);
+            upper.l = std::max(
+                upper.l, static_cast<int>(ceil(pos(2) + frac_radius(2))) + 1);
 
             lower.h = std::min(
-                lower.h, static_cast<int>(floor(pos(0) - frac_radius(0))));
+                lower.h, static_cast<int>(floor(pos(0) - frac_radius(0))) - 1);
             lower.k = std::min(
-                lower.k, static_cast<int>(floor(pos(1) - frac_radius(1))));
+                lower.k, static_cast<int>(floor(pos(1) - frac_radius(1))) - 1);
             lower.l = std::min(
-                lower.l, static_cast<int>(floor(pos(2) - frac_radius(2))));
+                lower.l, static_cast<int>(floor(pos(2) - frac_radius(2))) - 1);
         }
     }