diff --git a/include/occ/core/dimer.h b/include/occ/core/dimer.h index 378928de6..0d34eb2d1 100644 --- 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 index bef05beb1..05e181f68 100644 --- 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 - 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::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::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 - compute_JK(const MolecularOrbitals &mo, - double precision = std::numeric_limits::epsilon(), - const Mat &Schwarz = Mat()) { + inline std::pair 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::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 index 439b59486..5c8a1f57e 100644 --- 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 index fbf1ebfcd..39ceeaff4 100644 --- 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 index 3b7c88e0b..93aa8a5ee 100644 --- 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::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 - compute_JK(const MolecularOrbitals &mo, - double precision = std::numeric_limits::epsilon(), - const Mat &Schwarz = Mat()) const; + std::pair compute_JK(const MolecularOrbitals &mo, + const Mat &Schwarz = Mat()) const; Mat compute_J(const MolecularOrbitals &mo, - double precision = std::numeric_limits::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 index 958bb73b6..e63131ffe 100644 --- 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 index 17bd79f23..77bb2ba1c 100644 --- 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 V_LLt; diff --git a/include/occ/qm/scf.h b/include/occ/qm/scf.h index b69b63176..1eca0d93e 100644 --- a/include/occ/qm/scf.h +++ b/include/occ/qm/scf.h @@ -490,15 +490,8 @@ template 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::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 index 89e51413b..13a7d805a 100644 --- a/include/occ/solvent/solvation_correction.h +++ b/include/occ/solvent/solvation_correction.h @@ -147,6 +147,11 @@ template 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 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 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 class SolvationCorrectedProcedure { } Mat compute_fock(const MolecularOrbitals &mo, - double precision = std::numeric_limits::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 index 45bb79d2d..d83cf9e00 100644 --- 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 index 12486bbeb..86d6838a6 100644 --- 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(ceil(pos(0) + frac_radius(0)))); - upper.k = - std::max(upper.k, static_cast(ceil(pos(1) + frac_radius(1)))); - upper.l = - std::max(upper.l, static_cast(ceil(pos(2) + frac_radius(2)))); - - lower.h = - std::min(lower.h, static_cast(floor(pos(0) - frac_radius(0)))); - lower.k = - std::min(lower.k, static_cast(floor(pos(1) - frac_radius(1)))); - lower.l = - std::min(lower.l, static_cast(floor(pos(2) - frac_radius(2)))); + upper.h = std::max(upper.h, + static_cast(ceil(pos(0) + frac_radius(0))) + 1); + upper.k = std::max(upper.k, + static_cast(ceil(pos(1) + frac_radius(1))) + 1); + upper.l = std::max(upper.l, + static_cast(ceil(pos(2) + frac_radius(2))) + 1); + + lower.h = std::min( + lower.h, static_cast(floor(pos(0) - frac_radius(0))) - 1); + lower.k = std::min( + lower.k, static_cast(floor(pos(1) - frac_radius(1))) - 1); + lower.l = std::min( + lower.l, static_cast(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(ceil(pos(0) + frac_radius(0)))); - upper.k = std::max(upper.k, - static_cast(ceil(pos(1) + frac_radius(1)))); - upper.l = std::max(upper.l, - static_cast(ceil(pos(2) + frac_radius(2)))); + upper.h = std::max( + upper.h, static_cast(ceil(pos(0) + frac_radius(0))) + 1); + upper.k = std::max( + upper.k, static_cast(ceil(pos(1) + frac_radius(1))) + 1); + upper.l = std::max( + upper.l, static_cast(ceil(pos(2) + frac_radius(2))) + 1); lower.h = std::min( - lower.h, static_cast(floor(pos(0) - frac_radius(0)))); + lower.h, static_cast(floor(pos(0) - frac_radius(0))) - 1); lower.k = std::min( - lower.k, static_cast(floor(pos(1) - frac_radius(1)))); + lower.k, static_cast(floor(pos(1) - frac_radius(1))) - 1); lower.l = std::min( - lower.l, static_cast(floor(pos(2) - frac_radius(2)))); + lower.l, static_cast(floor(pos(2) - frac_radius(2))) - 1); } } diff --git a/src/interaction/pairinteraction.cpp b/src/interaction/pairinteraction.cpp index 02a9b3d2d..a76c58f2b 100644 --- a/src/interaction/pairinteraction.cpp +++ b/src/interaction/pairinteraction.cpp @@ -35,12 +35,11 @@ void compute_ce_model_energies(Wavefunction &wfn, Proc &proc, wfn.H = wfn.V + wfn.T; if (params.neglect_exchange) { occ::log::debug("neglecting K, only computing J"); - wfn.J = proc.compute_J(wfn.mo, params.precision, params.Schwarz); + wfn.J = proc.compute_J(wfn.mo, params.Schwarz); wfn.K = Mat::Zero(wfn.J.rows(), wfn.J.cols()); } else { occ::log::debug("computing J with K"); - std::tie(wfn.J, wfn.K) = - proc.compute_JK(wfn.mo, params.precision, params.Schwarz); + std::tie(wfn.J, wfn.K) = proc.compute_JK(wfn.mo, params.Schwarz); } wfn.energy.coulomb = expectation(wfn.mo.D, wfn.J); if constexpr (std::is_same::value) { @@ -72,11 +71,10 @@ void compute_ce_model_energies(Wavefunction &wfn, Proc &proc, wfn.energy.nuclear_attraction = 2 * expectation(wfn.mo.D, wfn.V); wfn.energy.kinetic = 2 * expectation(wfn.mo.D, wfn.T); if (params.neglect_exchange) { - wfn.J = proc.compute_J(wfn.mo, params.precision, params.Schwarz); + wfn.J = proc.compute_J(wfn.mo, params.Schwarz); wfn.K = Mat::Zero(wfn.J.rows(), wfn.J.cols()); } else { - std::tie(wfn.J, wfn.K) = - proc.compute_JK(wfn.mo, params.precision, params.Schwarz); + std::tie(wfn.J, wfn.K) = proc.compute_JK(wfn.mo, params.Schwarz); } wfn.energy.coulomb = expectation(wfn.mo.D, wfn.J); if constexpr (std::is_same::value) { @@ -159,14 +157,12 @@ void CEModelInteraction::set_use_xdm_dimer_parameters(bool value) { } void CEModelInteraction::compute_monomer_energies(Wavefunction &wfn) const { - constexpr double precision = std::numeric_limits::epsilon(); HartreeFock hf(wfn.basis); if (m_use_density_fitting) { occ::log::debug("Setting DF basis: def2-universal-jkfit"); hf.set_density_fitting_basis("def2-universal-jkfit"); } CEMonomerCalculationParameters params; - params.precision = precision; params.Schwarz = hf.compute_schwarz_ints(); params.xdm = m_scale_factors.xdm; @@ -199,7 +195,6 @@ CEEnergyComponents CEModelInteraction::operator()(Wavefunction &A, Wavefunction &B) const { using occ::disp::ce_model_dispersion_energy; using occ::qm::Energy; - constexpr double precision = std::numeric_limits::epsilon(); HartreeFock hf_a(A.basis); HartreeFock hf_b(B.basis); @@ -210,12 +205,10 @@ CEEnergyComponents CEModelInteraction::operator()(Wavefunction &A, } CEMonomerCalculationParameters params_a; - params_a.precision = precision; params_a.Schwarz = hf_a.compute_schwarz_ints(); params_a.xdm = m_scale_factors.xdm; CEMonomerCalculationParameters params_b; - params_b.precision = precision; params_b.Schwarz = hf_b.compute_schwarz_ints(); params_b.xdm = m_scale_factors.xdm; @@ -236,7 +229,6 @@ CEEnergyComponents CEModelInteraction::operator()(Wavefunction &A, // basis and atoms auto hf_AB = HartreeFock(ABn.basis); CEMonomerCalculationParameters params_ab; - params_ab.precision = precision; params_ab.Schwarz = hf_AB.compute_schwarz_ints(); params_ab.xdm = m_scale_factors.xdm && m_use_xdm_dimer_parameters; @@ -378,7 +370,6 @@ CEEnergyComponents CEModelInteraction::dft_pair(const std::string &functional, using occ::dft::DFT; using occ::disp::ce_model_dispersion_energy; using occ::qm::Energy; - constexpr double precision = std::numeric_limits::epsilon(); occ::io::BeckeGridSettings grid_settings{110, 30, 30, 1e-6}; @@ -390,7 +381,6 @@ CEEnergyComponents CEModelInteraction::dft_pair(const std::string &functional, } CEMonomerCalculationParameters params_a; - params_a.precision = precision; params_a.Schwarz = dft_a.compute_schwarz_ints(); params_a.xdm = m_scale_factors.xdm; @@ -401,7 +391,6 @@ CEEnergyComponents CEModelInteraction::dft_pair(const std::string &functional, } CEMonomerCalculationParameters params_b; - params_b.precision = precision; params_b.Schwarz = dft_b.compute_schwarz_ints(); params_b.xdm = m_scale_factors.xdm; @@ -422,7 +411,6 @@ CEEnergyComponents CEModelInteraction::dft_pair(const std::string &functional, // basis and atoms auto dft_AB = DFT(functional, ABn.basis, grid_settings); CEMonomerCalculationParameters params_ab; - params_ab.precision = precision; params_ab.Schwarz = dft_AB.compute_schwarz_ints(); params_ab.xdm = m_scale_factors.xdm && m_use_xdm_dimer_parameters; diff --git a/src/main/occ_cg.cpp b/src/main/occ_cg.cpp index d8c675d92..daa2b7e2a 100644 --- a/src/main/occ_cg.cpp +++ b/src/main/occ_cg.cpp @@ -165,7 +165,6 @@ void compute_monomer_energies(const std::string &basename, std::cout << std::flush; HartreeFock hf(wfn.basis); occ::interaction::CEMonomerCalculationParameters params; - params.precision = 1e-8; params.Schwarz = hf.compute_schwarz_ints(); occ::interaction::compute_ce_model_energies(wfn, hf, params); occ::log::debug("Writing monomer energies to {}", diff --git a/src/main/occ_scf.cpp b/src/main/occ_scf.cpp index fa9845368..83ce521b8 100644 --- a/src/main/occ_scf.cpp +++ b/src/main/occ_scf.cpp @@ -111,6 +111,9 @@ CLI::App *add_scf_subcommand(CLI::App &app) { }, "use unrestricted SCF"); + scf->add_option("--integral-precision,--integral_precision", + config->method.integral_precision, + "cutoff for integral screening"); // dft grid scf->add_option("--dft-grid-max-angular,--dft_grid_max_angular", config->method.dft_grid.max_angular_points, diff --git a/src/main/single_point.cpp b/src/main/single_point.cpp index 8f0e9e334..325d9263f 100644 --- a/src/main/single_point.cpp +++ b/src/main/single_point.cpp @@ -87,6 +87,11 @@ Wavefunction run_method(Molecule &m, const occ::qm::AOBasis &basis, if (!config.basis.df_name.empty()) proc.set_density_fitting_basis(config.basis.df_name); + + occ::log::trace("Setting integral precision: {}", + config.method.integral_precision); + proc.set_precision(config.method.integral_precision); + SCF scf(proc, SK); occ::log::trace("Setting system charge: {}", config.electronic.charge); occ::log::trace("Setting system multiplicity: {}", diff --git a/src/occpy.cpp b/src/occpy.cpp index 95c0fc69c..f7d3962c9 100644 --- a/src/occpy.cpp +++ b/src/occpy.cpp @@ -176,6 +176,7 @@ NB_MODULE(_occpy, m) { .def("overlap_matrix", &HartreeFock::compute_overlap_matrix) .def("nuclear_repulsion", &HartreeFock::nuclear_repulsion_energy) .def("scf", [](HartreeFock &hf) { return HF(hf); }) + .def("set_precision", &HartreeFock::set_precision) .def("coulomb_matrix", [](const HartreeFock &hf, const MolecularOrbitals &mo) { return hf.compute_J(mo); @@ -210,6 +211,7 @@ NB_MODULE(_occpy, m) { .def("kinetic_matrix", &DFT::compute_kinetic_matrix) .def("overlap_matrix", &DFT::compute_overlap_matrix) .def("nuclear_repulsion", &DFT::nuclear_repulsion_energy) + .def("set_precision", &HartreeFock::set_precision) .def("set_method", &DFT::set_method, nb::arg("method_string"), nb::arg("unrestricted") = false) .def("set_unrestricted", &DFT::set_unrestricted) diff --git a/src/qm/hf.cpp b/src/qm/hf.cpp index 1a9fb0728..2e7a173f8 100644 --- a/src/qm/hf.cpp +++ b/src/qm/hf.cpp @@ -72,7 +72,7 @@ double HartreeFock::nuclear_point_charge_interaction_energy( return etot; } -Mat HartreeFock::compute_fock(const MolecularOrbitals &mo, double precision, +Mat HartreeFock::compute_fock(const MolecularOrbitals &mo, const Mat &Schwarz) const { if (m_df_engine) { return m_df_engine->fock_operator(mo); @@ -114,7 +114,7 @@ Mat HartreeFock::compute_fock_mixed_basis(const MolecularOrbitals &mo_minbs, } std::pair HartreeFock::compute_JK(const MolecularOrbitals &mo, - double precision, + const Mat &Schwarz) const { if (m_df_engine) { return m_df_engine->coulomb_and_exchange(mo); @@ -123,7 +123,7 @@ std::pair HartreeFock::compute_JK(const MolecularOrbitals &mo, } } -Mat HartreeFock::compute_J(const MolecularOrbitals &mo, double precision, +Mat HartreeFock::compute_J(const MolecularOrbitals &mo, const Mat &Schwarz) const { if (m_df_engine) { return m_df_engine->coulomb(mo); diff --git a/src/qm/integral_engine.cpp b/src/qm/integral_engine.cpp index da1ed75f3..408c21f2d 100644 --- a/src/qm/integral_engine.cpp +++ b/src/qm/integral_engine.cpp @@ -993,7 +993,7 @@ void evaluate_four_center(Lambda &f, cint::IntegralEnvironment &env, template Mat fock_operator_kernel(cint::IntegralEnvironment &env, const AOBasis &basis, const ShellPairList &shellpairs, - const MolecularOrbitals &mo, + const MolecularOrbitals &mo, double precision = 1e-12, const Mat &Schwarz = Mat()) { using Result = IntegralEngine::IntegralResult<4>; auto nthreads = occ::parallel::get_num_threads(); @@ -1029,7 +1029,7 @@ Mat fock_operator_kernel(cint::IntegralEnvironment &env, const AOBasis &basis, }; auto lambda = [&](int thread_id) { evaluate_four_center(f, env, basis, shellpairs, Dnorm, - Schwarz, 1e-12, thread_id); + Schwarz, precision, thread_id); }; occ::timing::start(occ::timing::category::fock); occ::parallel::parallel_do(lambda); @@ -1064,7 +1064,7 @@ Mat fock_operator_kernel(cint::IntegralEnvironment &env, const AOBasis &basis, template Mat coulomb_kernel(cint::IntegralEnvironment &env, const AOBasis &basis, const ShellPairList &shellpairs, const MolecularOrbitals &mo, - const Mat &Schwarz = Mat()) { + double precision = 1e-12, const Mat &Schwarz = Mat()) { using Result = IntegralEngine::IntegralResult<4>; auto nthreads = occ::parallel::get_num_threads(); constexpr Op op = Op::coulomb; @@ -1098,7 +1098,7 @@ Mat coulomb_kernel(cint::IntegralEnvironment &env, const AOBasis &basis, }; auto lambda = [&](int thread_id) { evaluate_four_center(f, env, basis, shellpairs, Dnorm, - Schwarz, 1e-12, thread_id); + Schwarz, precision, thread_id); }; occ::timing::start(occ::timing::category::fock); occ::parallel::parallel_do(lambda); @@ -1138,6 +1138,7 @@ std::pair coulomb_and_exchange_kernel(cint::IntegralEnvironment &env, const AOBasis &basis, const ShellPairList &shellpairs, const MolecularOrbitals &mo, + double precision = 1e-12, const Mat &Schwarz = Mat()) { using Result = IntegralEngine::IntegralResult<4>; auto nthreads = occ::parallel::get_num_threads(); @@ -1175,7 +1176,7 @@ std::pair coulomb_and_exchange_kernel(cint::IntegralEnvironment &env, }; auto lambda = [&](int thread_id) { evaluate_four_center(f, env, basis, shellpairs, Dnorm, - Schwarz, 1e-12, thread_id); + Schwarz, precision, thread_id); }; occ::timing::start(occ::timing::category::fock); occ::parallel::parallel_do(lambda); @@ -1246,28 +1247,28 @@ Mat IntegralEngine::fock_operator(SpinorbitalKind sk, case R: if (spherical) { return fock_operator_kernel(m_env, m_aobasis, m_shellpairs, - mo, Schwarz); + mo, m_precision, Schwarz); } else { return fock_operator_kernel(m_env, m_aobasis, m_shellpairs, - mo, Schwarz); + mo, m_precision, Schwarz); } break; case U: if (spherical) { return fock_operator_kernel(m_env, m_aobasis, m_shellpairs, - mo, Schwarz); + mo, m_precision, Schwarz); } else { return fock_operator_kernel(m_env, m_aobasis, m_shellpairs, - mo, Schwarz); + mo, m_precision, Schwarz); } case G: if (spherical) { return fock_operator_kernel(m_env, m_aobasis, m_shellpairs, - mo, Schwarz); + mo, m_precision, Schwarz); } else { return fock_operator_kernel(m_env, m_aobasis, m_shellpairs, - mo, Schwarz); + mo, m_precision, Schwarz); } } } @@ -1285,28 +1286,28 @@ Mat IntegralEngine::coulomb(SpinorbitalKind sk, const MolecularOrbitals &mo, case R: if (spherical) { return coulomb_kernel(m_env, m_aobasis, m_shellpairs, mo, - Schwarz); + m_precision, Schwarz); } else { return coulomb_kernel(m_env, m_aobasis, m_shellpairs, mo, - Schwarz); + m_precision, Schwarz); } break; case U: if (spherical) { return coulomb_kernel(m_env, m_aobasis, m_shellpairs, mo, - Schwarz); + m_precision, Schwarz); } else { return coulomb_kernel(m_env, m_aobasis, m_shellpairs, mo, - Schwarz); + m_precision, Schwarz); } case G: if (spherical) { return coulomb_kernel(m_env, m_aobasis, m_shellpairs, mo, - Schwarz); + m_precision, Schwarz); } else { return coulomb_kernel(m_env, m_aobasis, m_shellpairs, mo, - Schwarz); + m_precision, Schwarz); } } } @@ -1324,28 +1325,28 @@ std::pair IntegralEngine::coulomb_and_exchange( case R: if (spherical) { return coulomb_and_exchange_kernel( - m_env, m_aobasis, m_shellpairs, mo, Schwarz); + m_env, m_aobasis, m_shellpairs, mo, m_precision, Schwarz); } else { return coulomb_and_exchange_kernel( - m_env, m_aobasis, m_shellpairs, mo, Schwarz); + m_env, m_aobasis, m_shellpairs, mo, m_precision, Schwarz); } break; case U: if (spherical) { return coulomb_and_exchange_kernel( - m_env, m_aobasis, m_shellpairs, mo, Schwarz); + m_env, m_aobasis, m_shellpairs, mo, m_precision, Schwarz); } else { return coulomb_and_exchange_kernel( - m_env, m_aobasis, m_shellpairs, mo, Schwarz); + m_env, m_aobasis, m_shellpairs, mo, m_precision, Schwarz); } case G: if (spherical) { return coulomb_and_exchange_kernel( - m_env, m_aobasis, m_shellpairs, mo, Schwarz); + m_env, m_aobasis, m_shellpairs, mo, m_precision, Schwarz); } else { return coulomb_and_exchange_kernel( - m_env, m_aobasis, m_shellpairs, mo, Schwarz); + m_env, m_aobasis, m_shellpairs, mo, m_precision, Schwarz); } } } diff --git a/src/qm/integral_engine_df.cpp b/src/qm/integral_engine_df.cpp index bdf5e2969..6e5210097 100644 --- a/src/qm/integral_engine_df.cpp +++ b/src/qm/integral_engine_df.cpp @@ -1464,4 +1464,9 @@ void IntegralEngineDF::set_range_separated_omega(double omega) { set_integral_policy(Policy::Direct); } +void IntegralEngineDF::set_precision(double precision) { + m_precision = precision; + m_ao_engine.set_precision(precision); + m_aux_engine.set_precision(precision); +} } // namespace occ::qm diff --git a/src/qm/wavefunction.cpp b/src/qm/wavefunction.cpp index e8780458f..8ce016da9 100644 --- a/src/qm/wavefunction.cpp +++ b/src/qm/wavefunction.cpp @@ -365,10 +365,7 @@ void Wavefunction::save(FchkWriter &fchk) { } // TODO fix this is wrong - fchk.set_scalar("Virial ratio", - -(energy.nuclear_repulsion + energy.nuclear_attraction + - energy.coulomb + energy.exchange) / - energy.kinetic); + fchk.set_scalar("Virial ratio", -2 * energy.total / energy.kinetic); fchk.set_scalar("SCF ratio", energy.total); fchk.set_scalar("Total ratio", energy.total); diff --git a/tests/dft_tests.cpp b/tests/dft_tests.cpp index d43088e83..cde0bf2aa 100644 --- a/tests/dft_tests.cpp +++ b/tests/dft_tests.cpp @@ -334,7 +334,7 @@ TEST_CASE("Water seminumerical exchange approximation", "[scf]") { fmt::print("Compute K SGX done\n"); occ::Mat Jexact, Kexact; sw.start(1); - std::tie(Jexact, Kexact) = hf.compute_JK(scf.mo, 1e-12, occ::Mat()); + std::tie(Jexact, Kexact) = hf.compute_JK(scf.mo, occ::Mat()); sw.stop(1); int i, j; fmt::print("K - Kexact: {:12.8f}\n",