From 06121a6e3a5b9151677edc3ce50690b82b77bbaf Mon Sep 17 00:00:00 2001 From: Peter Spackman Date: Tue, 21 Jan 2025 18:57:46 +0800 Subject: [PATCH] Big refactor of interaction energy calculation --- CMakeLists.txt | 2 +- include/occ/core/log.h | 9 +- include/occ/interaction/ce_energy_model.h | 33 + include/occ/interaction/energy_model_base.h | 16 + .../lattice_convergence_settings.h | 20 + include/occ/interaction/lattice_energy.h | 36 + include/occ/interaction/pair_energy.h | 63 -- include/occ/interaction/pair_energy_store.h | 62 ++ .../occ/interaction/wavefunction_transform.h | 24 + include/occ/interaction/wolf.h | 38 +- include/occ/interaction/xtb_energy_model.h | 22 + include/occ/main/occ_cg.h | 2 +- include/occ/main/occ_elat.h | 2 +- pyproject.toml | 2 +- src/core/log.cpp | 86 ++- src/driver/crystal_growth.cpp | 30 +- src/interaction/CMakeLists.txt | 5 + src/interaction/ce_energy_model.cpp | 78 +++ src/interaction/lattice_energy.cpp | 171 +++++ src/interaction/pair_energy.cpp | 639 +----------------- src/interaction/pair_energy_store.cpp | 148 ++++ src/interaction/wavefunction_transform.cpp | 59 ++ src/interaction/wolf.cpp | 96 ++- src/interaction/xtb_energy_model.cpp | 40 ++ src/main/make_atomic_pair_potentials.cpp | 4 +- src/main/occ.cpp | 4 +- src/main/occ_elat.cpp | 55 +- src/occpy.cpp | 21 +- src/occpy/__init__.py | 16 +- tests/interaction_tests.cpp | 7 +- 30 files changed, 1017 insertions(+), 773 deletions(-) create mode 100644 include/occ/interaction/ce_energy_model.h create mode 100644 include/occ/interaction/energy_model_base.h create mode 100644 include/occ/interaction/lattice_convergence_settings.h create mode 100644 include/occ/interaction/lattice_energy.h create mode 100644 include/occ/interaction/pair_energy_store.h create mode 100644 include/occ/interaction/wavefunction_transform.h create mode 100644 include/occ/interaction/xtb_energy_model.h create mode 100644 src/interaction/ce_energy_model.cpp create mode 100644 src/interaction/lattice_energy.cpp create mode 100644 src/interaction/pair_energy_store.cpp create mode 100644 src/interaction/wavefunction_transform.cpp create mode 100644 src/interaction/xtb_energy_model.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index c54c3173a..7b42af0cb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.16.0 FATAL_ERROR) set(PROJECT_VERSION_MAJOR "0") set(PROJECT_VERSION_MINOR "6") -set(PROJECT_VERSION_PATCH "12") +set(PROJECT_VERSION_PATCH "13") project(occ VERSION "${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}.${PROJECT_VERSION_PATCH}" diff --git a/include/occ/core/log.h b/include/occ/core/log.h index c15bfb86e..b5ba50f06 100644 --- a/include/occ/core/log.h +++ b/include/occ/core/log.h @@ -1,12 +1,12 @@ #pragma once #include +#include namespace occ::log { using spdlog::critical; using spdlog::debug; using spdlog::error; using spdlog::info; -using spdlog::set_level; using spdlog::trace; using spdlog::warn; @@ -19,7 +19,10 @@ using spdlog::level::trace; using spdlog::level::warn; } // namespace level -void setup_logging(const std::string &verbosity); -void setup_logging(int verbosity); +void set_log_level(const std::string &verbosity); +void set_log_level(spdlog::level::level_enum level); +void set_log_level(int verbosity); + +void set_log_file(const std::string &filename); } // namespace occ::log diff --git a/include/occ/interaction/ce_energy_model.h b/include/occ/interaction/ce_energy_model.h new file mode 100644 index 000000000..006153bb3 --- /dev/null +++ b/include/occ/interaction/ce_energy_model.h @@ -0,0 +1,33 @@ +#pragma once +#include +#include + +namespace occ::interaction { + +class CEEnergyModel : public EnergyModelBase { +public: + CEEnergyModel(const crystal::Crystal &crystal, + const std::vector &wfns_a, + const std::vector &wfns_b = {}); + + void set_model_name(const std::string &model_name) { + m_model_name = model_name; + } + + CEEnergyComponents compute_energy(const core::Dimer &dimer) override; + Mat3N compute_electric_field(const core::Dimer &dimer) override; + const std::vector &partial_charges() const override; + double coulomb_scale_factor() const override; + +private: + Wavefunction prepare_wavefunction(const core::Molecule &mol, + const Wavefunction &wfn) const; + + crystal::Crystal m_crystal; + std::string m_model_name{"ce-b3lyp"}; + std::vector m_wavefunctions_a; + std::vector m_wavefunctions_b; + mutable std::vector m_partial_charges; +}; + +} // namespace occ::interaction diff --git a/include/occ/interaction/energy_model_base.h b/include/occ/interaction/energy_model_base.h new file mode 100644 index 000000000..509b6d3e6 --- /dev/null +++ b/include/occ/interaction/energy_model_base.h @@ -0,0 +1,16 @@ +#pragma once +#include +#include + +namespace occ::interaction { + +class EnergyModelBase { +public: + virtual ~EnergyModelBase() = default; + virtual CEEnergyComponents compute_energy(const core::Dimer &dimer) = 0; + virtual Mat3N compute_electric_field(const core::Dimer &dimer) = 0; + virtual const std::vector &partial_charges() const = 0; + virtual double coulomb_scale_factor() const = 0; +}; + +} // namespace occ::interaction diff --git a/include/occ/interaction/lattice_convergence_settings.h b/include/occ/interaction/lattice_convergence_settings.h new file mode 100644 index 000000000..0c8af432c --- /dev/null +++ b/include/occ/interaction/lattice_convergence_settings.h @@ -0,0 +1,20 @@ +#pragma once +#include + +namespace occ::interaction { + +struct LatticeConvergenceSettings { + double min_radius{3.8}; // angstroms + double max_radius{30.0}; // angstroms + double radius_increment{3.8}; // angstroms + double energy_tolerance{1.0}; // kj/mol + bool wolf_sum{false}; + bool crystal_field_polarization{false}; + std::string model_name{"ce-b3lyp"}; + std::string crystal_filename{""}; + std::string output_json_filename{""}; + bool spherical_basis{false}; + std::string charge_string; +}; + +} // namespace occ::interaction diff --git a/include/occ/interaction/lattice_energy.h b/include/occ/interaction/lattice_energy.h new file mode 100644 index 000000000..70499be2a --- /dev/null +++ b/include/occ/interaction/lattice_energy.h @@ -0,0 +1,36 @@ +#pragma once +#include +#include +#include +#include + +namespace occ::interaction { + +struct LatticeEnergyResult { + double lattice_energy{0.0}; + occ::crystal::CrystalDimers dimers; + std::vector energy_components; +}; + +class LatticeEnergyCalculator { +public: + LatticeEnergyCalculator(std::unique_ptr model, + const crystal::Crystal &crystal, + const std::string &basename, + LatticeConvergenceSettings settings); + + LatticeEnergyResult compute(); + +private: + void initialize_wolf_sum(); + bool is_converged(double current, double previous) const; + void report_convergence_statistics(const CEEnergyComponents &total); + + std::unique_ptr m_energy_model; + crystal::Crystal m_crystal; + std::string m_basename; + LatticeConvergenceSettings m_settings; + std::unique_ptr m_wolf_sum; +}; + +} // namespace occ::interaction diff --git a/include/occ/interaction/pair_energy.h b/include/occ/interaction/pair_energy.h index 4a8d21325..e4717c60c 100644 --- a/include/occ/interaction/pair_energy.h +++ b/include/occ/interaction/pair_energy.h @@ -4,14 +4,9 @@ #include #include #include -#include namespace occ::interaction { -using occ::core::Dimer; -using occ::crystal::Crystal; -using occ::qm::Wavefunction; - struct PairEnergy { struct Monomer { occ::qm::Wavefunction wfn; @@ -28,64 +23,6 @@ struct PairEnergy { CEEnergyComponents energy; }; -struct PairEnergyStore { - enum class Kind { JSON, Xyz, Memory }; - - bool save(int id, const Dimer &d, const CEEnergyComponents &); - bool load(int id, const Dimer &d, CEEnergyComponents &); - std::string dimer_filename(int id, const Dimer &d); - - Kind kind{Kind::JSON}; - std::string name; -}; - -bool load_dimer_energy(const std::string &, CEEnergyComponents &); -bool write_xyz_dimer(const std::string &, const Dimer &, - std::optional energies = {}); - -CEEnergyComponents ce_model_energy(const Dimer &dimer, - const std::vector &wfns_a, - const std::vector &wfns_b, - const Crystal &crystal); - -std::vector -ce_model_energies(const Crystal &crystal, const std::vector &dimers, - const std::vector &wfns_a, - const std::vector &wfns_b, - const std::string &basename = "dimer"); - -struct LatticeConvergenceSettings { - double min_radius{3.8}; // angstroms - double max_radius{30.0}; // angstroms - double radius_increment{3.8}; // angstroms - double energy_tolerance{1.0}; // kj/mol - bool wolf_sum{false}; - bool crystal_field_polarization{false}; - std::string model_name{"ce-b3lyp"}; - std::string crystal_filename{""}; - std::string output_json_filename{""}; - bool spherical_basis{false}; - std::string charge_string; -}; - -struct LatticeEnergyResult { - double lattice_energy{0.0}; - occ::crystal::CrystalDimers dimers; - std::vector energy_components; -}; - -LatticeEnergyResult -converged_lattice_energies(const Crystal &crystal, - const std::vector &wfns_a, - const std::vector &wfns_b, - const std::string &basename = "crystal_dimer", - const LatticeConvergenceSettings conv = {}); - -LatticeEnergyResult -converged_xtb_lattice_energies(const Crystal &crystal, - const std::string &basename = "crystal_dimer", - const LatticeConvergenceSettings conv = {}); - Wavefunction load_wavefunction(const std::string &filename); } // namespace occ::interaction diff --git a/include/occ/interaction/pair_energy_store.h b/include/occ/interaction/pair_energy_store.h new file mode 100644 index 000000000..066681858 --- /dev/null +++ b/include/occ/interaction/pair_energy_store.h @@ -0,0 +1,62 @@ +#pragma once +#include +#include +#include +#include +#include + +namespace occ::interaction { + +class PairEnergyStoreBase { +public: + virtual ~PairEnergyStoreBase() = default; + virtual bool save(int id, const core::Dimer &d, + const CEEnergyComponents &e) = 0; + virtual bool load(int id, const core::Dimer &d, CEEnergyComponents &e) = 0; + virtual std::string dimer_filename(int id, const core::Dimer &d) const = 0; +}; + +class MemoryPairEnergyStore : public PairEnergyStoreBase { +public: + bool save(int id, const core::Dimer &d, const CEEnergyComponents &e) override; + bool load(int id, const core::Dimer &d, CEEnergyComponents &e) override; + std::string dimer_filename(int id, const core::Dimer &d) const override; + +private: + std::unordered_map energy_store; +}; + +class FileSystemPairEnergyStore : public PairEnergyStoreBase { +public: + enum class Format { JSON, XYZ }; + explicit FileSystemPairEnergyStore(std::string path, + Format format = Format::XYZ) + : base_path(std::move(path)), format(format) {} + + bool save(int id, const core::Dimer &d, const CEEnergyComponents &e) override; + bool load(int id, const core::Dimer &d, CEEnergyComponents &e) override; + std::string dimer_filename(int id, const core::Dimer &d) const override; + +private: + std::string base_path; + Format format; +}; + +class PairEnergyStore { +public: + enum class Kind { JSON, XYZ, Memory }; + + explicit PairEnergyStore(Kind kind = Kind::XYZ, std::string name = ""); + + bool save(int id, const core::Dimer &d, const CEEnergyComponents &e); + bool load(int id, const core::Dimer &d, CEEnergyComponents &e); + std::string dimer_filename(int id, const core::Dimer &d) const; + + Kind kind{Kind::XYZ}; + std::string name; + +private: + std::unique_ptr store; +}; + +} // namespace occ::interaction diff --git a/include/occ/interaction/wavefunction_transform.h b/include/occ/interaction/wavefunction_transform.h new file mode 100644 index 000000000..f4dc4c8ce --- /dev/null +++ b/include/occ/interaction/wavefunction_transform.h @@ -0,0 +1,24 @@ +#pragma once +#include +#include + +namespace occ::interaction::transform { + +struct TransformResult { + Mat3 rotation{Mat3::Identity()}; + Vec3 translation{Vec3::Zero()}; + double rmsd{0.0}; +}; + +class WavefunctionTransformer { +public: + static TransformResult calculate_transform(const qm::Wavefunction &wfn, + const core::Molecule &mol, + const crystal::Crystal &crystal); + +private: + static double compute_position_rmsd(const Mat3N &pos_ref, + const Mat3N &pos_transformed); +}; + +} // namespace occ::interaction::transform diff --git a/include/occ/interaction/wolf.h b/include/occ/interaction/wolf.h index 30d8d4ac0..bb05bb5b9 100644 --- a/include/occ/interaction/wolf.h +++ b/include/occ/interaction/wolf.h @@ -1,15 +1,45 @@ #pragma once #include +#include +#include namespace occ::interaction { -struct WolfParams { - double cutoff{16.0}; // Angstroms; - double eta{0.2}; // 1/Angstroms +struct WolfParameters { + double cutoff{16.0}; + double alpha{0.2}; }; double wolf_coulomb_energy(double qi, const Vec3 &pi, Eigen::Ref qj, Eigen::Ref pj, - const WolfParams ¶ms = {}); + const WolfParameters ¶ms); + +class WolfSum { +public: + WolfSum(WolfParameters params = {}); + + void initialize(const crystal::Crystal &crystal, + const EnergyModelBase &energy_model); + + double compute_correction( + const std::vector &charge_energies, + const std::vector &model_energies) const; + + Mat3N &electric_field_for_molecule(size_t idx) { + return m_electric_field_values[idx]; + } + + inline const auto &asymmetric_charges() const { return m_asym_charges; } + +private: + void compute_self_energies(const crystal::Crystal &crystal); + void compute_wolf_energies(const crystal::Crystal &crystal); + + WolfParameters m_params; + Vec m_asym_charges; + std::vector m_charge_self_energies; + std::vector m_electric_field_values; + double m_total_energy{0.0}; +}; } // namespace occ::interaction diff --git a/include/occ/interaction/xtb_energy_model.h b/include/occ/interaction/xtb_energy_model.h new file mode 100644 index 000000000..c2cf7c0c1 --- /dev/null +++ b/include/occ/interaction/xtb_energy_model.h @@ -0,0 +1,22 @@ +#pragma once +#include +#include + +namespace occ::interaction { + +class XTBEnergyModel : public EnergyModelBase { +public: + explicit XTBEnergyModel(const crystal::Crystal &crystal); + + CEEnergyComponents compute_energy(const core::Dimer &dimer) override; + Mat3N compute_electric_field(const core::Dimer &dimer) override; + const std::vector &partial_charges() const override; + double coulomb_scale_factor() const override { return 1.0; } + +private: + crystal::Crystal m_crystal; + std::vector m_monomer_energies; + std::vector m_partial_charges; +}; + +} // namespace occ::interaction diff --git a/include/occ/main/occ_cg.h b/include/occ/main/occ_cg.h index 5ebb1edff..f00840963 100644 --- a/include/occ/main/occ_cg.h +++ b/include/occ/main/occ_cg.h @@ -2,7 +2,7 @@ #include #include #include -#include +#include namespace occ::main { diff --git a/include/occ/main/occ_elat.h b/include/occ/main/occ_elat.h index 7ffca0a47..3e3ef8a50 100644 --- a/include/occ/main/occ_elat.h +++ b/include/occ/main/occ_elat.h @@ -1,6 +1,6 @@ #pragma once #include -#include +#include namespace occ::main { diff --git a/pyproject.toml b/pyproject.toml index faff2f399..d1ad89278 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build" [project] name = "occpy" -version = "0.6.12" +version = "0.6.13" description = "A library for quantum chemistry" readme = "README.md" requires-python = ">=3.10" diff --git a/src/core/log.cpp b/src/core/log.cpp index da710d393..5bd42a510 100644 --- a/src/core/log.cpp +++ b/src/core/log.cpp @@ -1,51 +1,75 @@ +#include #include #include +#include namespace occ::log { +namespace { +std::shared_ptr current_logger = spdlog::default_logger(); -void setup_logging(const std::string &verbosity) { - auto level = occ::log::level::info; +spdlog::level::level_enum verbosity_to_level(const std::string &verbosity) { std::string level_lower = occ::util::to_lower_copy(verbosity); if (level_lower == "debug") - level = occ::log::level::trace; - else if (level_lower == "normal") - level = occ::log::level::info; - else if (level_lower == "verbose") - level = occ::log::level::debug; - else if (level_lower == "minimal") - level = occ::log::level::warn; - else if (level_lower == "silent") - level = occ::log::level::critical; - occ::log::set_level(level); - spdlog::set_level(level); - // store the last 32 debug messages in a buffer - spdlog::enable_backtrace(32); - spdlog::set_pattern("%v"); + return spdlog::level::trace; + if (level_lower == "verbose") + return spdlog::level::debug; + if (level_lower == "minimal") + return spdlog::level::warn; + if (level_lower == "silent") + return spdlog::level::critical; + return spdlog::level::info; // default for "normal" and unknown values } -void setup_logging(int verbosity) { - auto level = occ::log::level::info; +spdlog::level::level_enum verbosity_to_level(int verbosity) { switch (verbosity) { case 4: - level = occ::log::level::trace; - break; + return spdlog::level::trace; case 3: - level = occ::log::level::debug; - break; + return spdlog::level::debug; case 1: - level = occ::log::level::warn; - break; + return spdlog::level::warn; case 0: - level = occ::log::level::critical; - break; + return spdlog::level::critical; default: - level = occ::log::level::info; - break; + return spdlog::level::info; } - occ::log::set_level(level); - spdlog::set_level(level); - // store the last 32 debug messages in a buffer +} +} // namespace + +void set_log_level(spdlog::level::level_enum level) { + current_logger->set_level(level); + spdlog::set_pattern("%v"); + spdlog::enable_backtrace(32); +} + +void set_log_level(const std::string &verbosity) { + auto level = verbosity_to_level(verbosity); + current_logger->set_level(level); + spdlog::set_pattern("%v"); spdlog::enable_backtrace(32); +} + +void set_log_level(int verbosity) { + auto level = verbosity_to_level(verbosity); + current_logger->set_level(level); spdlog::set_pattern("%v"); + spdlog::enable_backtrace(32); } + +void set_log_file(const std::string &filename) { + try { + auto file_logger = spdlog::basic_logger_mt("occ_logger", filename, + true); // true = truncate + current_logger = file_logger; + file_logger->set_level(current_logger->level()); + spdlog::set_default_logger(current_logger); + } catch (const spdlog::spdlog_ex &ex) { + spdlog::warn( + "Failed to create file logger: {}. Using existing logger instead.", + ex.what()); + } + spdlog::set_pattern("%v"); + spdlog::enable_backtrace(32); +} + } // namespace occ::log diff --git a/src/driver/crystal_growth.cpp b/src/driver/crystal_growth.cpp index 7b42040ba..ea2427d4d 100644 --- a/src/driver/crystal_growth.cpp +++ b/src/driver/crystal_growth.cpp @@ -1,9 +1,16 @@ #include #include #include +#include +#include +#include #include namespace fs = std::filesystem; +using occ::interaction::CEEnergyModel; +using occ::interaction::LatticeConvergenceSettings; +using occ::interaction::LatticeEnergyCalculator; +using occ::interaction::XTBEnergyModel; namespace occ::driver { @@ -349,16 +356,22 @@ void CEModelCrystalGrowthCalculator::converge_lattice_energy() { occ::log::info("Computing crystal interactions using {} wavefunctions", wfn_choice); - occ::interaction::LatticeConvergenceSettings convergence_settings; + LatticeConvergenceSettings convergence_settings; convergence_settings.model_name = opts.energy_model; convergence_settings.max_radius = opts.outer_radius; convergence_settings.wolf_sum = opts.use_wolf_sum; convergence_settings.crystal_field_polarization = opts.use_crystal_polarization; - auto result = occ::interaction::converged_lattice_energies( - m_crystal, inner_wavefunctions(), outer_wavefunctions(), opts.basename, - convergence_settings); + auto energy_model = std::make_unique( + m_crystal, inner_wavefunctions(), outer_wavefunctions()); + energy_model->set_model_name(opts.energy_model); + + LatticeEnergyCalculator calculator(std::move(energy_model), m_crystal, + opts.basename, convergence_settings); + + auto result = calculator.compute(); + m_full_dimers = result.dimers; m_dimer_energies = result.energy_components; @@ -570,8 +583,13 @@ void XTBCrystalGrowthCalculator::converge_lattice_energy() { m_full_dimers = m_crystal.symmetry_unique_dimers(opts.outer_radius); std::vector energies; - auto result = occ::interaction::converged_xtb_lattice_energies( - m_crystal, opts.basename, convergence_settings); + + LatticeEnergyCalculator calculator( + std::make_unique(m_crystal), m_crystal, opts.basename, + convergence_settings); + + auto result = calculator.compute(); + m_full_dimers = result.dimers; for (const auto &e : result.energy_components) { diff --git a/src/interaction/CMakeLists.txt b/src/interaction/CMakeLists.txt index d3c287272..97e58e5c7 100644 --- a/src/interaction/CMakeLists.txt +++ b/src/interaction/CMakeLists.txt @@ -1,12 +1,17 @@ add_library(occ_interaction + "${CMAKE_CURRENT_SOURCE_DIR}/ce_energy_model.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/coulomb.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/disp.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/interaction_json.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/lattice_energy.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/pair_energy.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/pair_energy_store.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/pairinteraction.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/pair_potential.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/polarization.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/wavefunction_transform.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/wolf.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/xtb_energy_model.cpp" ${OCC_INTERACTION_INCLUDE_FILES} ) diff --git a/src/interaction/ce_energy_model.cpp b/src/interaction/ce_energy_model.cpp new file mode 100644 index 000000000..c632aa064 --- /dev/null +++ b/src/interaction/ce_energy_model.cpp @@ -0,0 +1,78 @@ +#include +#include +#include + +namespace occ::interaction { + +CEEnergyModel::CEEnergyModel(const crystal::Crystal &crystal, + const std::vector &wfns_a, + const std::vector &wfns_b) + : m_crystal(crystal), m_wavefunctions_a(wfns_a), + m_wavefunctions_b(wfns_b.empty() ? wfns_a : wfns_b) {} + +Wavefunction +CEEnergyModel::prepare_wavefunction(const core::Molecule &mol, + const Wavefunction &wfn) const { + + auto transform = transform::WavefunctionTransformer::calculate_transform( + wfn, mol, m_crystal); + + Wavefunction result = wfn; + result.apply_transformation(transform.rotation, transform.translation); + + Mat3N pos_mol = mol.positions(); + Mat3N pos_transformed = result.positions() * units::BOHR_TO_ANGSTROM; + + occ::log::debug("Transformed wavefunction positions RMSD = {}", + (pos_transformed - pos_mol).norm()); + + return result; +} + +CEEnergyComponents CEEnergyModel::compute_energy(const core::Dimer &dimer) { + core::Molecule mol_A = dimer.a(); + core::Molecule mol_B = dimer.b(); + + Wavefunction A = m_wavefunctions_a[mol_A.asymmetric_molecule_idx()]; + Wavefunction B = m_wavefunctions_b[mol_B.asymmetric_molecule_idx()]; + + A = prepare_wavefunction(mol_A, A); + B = prepare_wavefunction(mol_B, B); + + auto model = ce_model_from_string(m_model_name); + CEModelInteraction interaction(model); + + auto result = interaction(A, B); + result.is_computed = true; + + occ::log::debug("Finished model energy calculation"); + return result; +} + +Mat3N CEEnergyModel::compute_electric_field(const core::Dimer &dimer) { + core::Molecule mol_A = dimer.a(); + core::Molecule mol_B = dimer.b(); + + Wavefunction B = m_wavefunctions_b[mol_B.asymmetric_molecule_idx()]; + B = prepare_wavefunction(mol_B, B); + + Mat3N pos_A_bohr = mol_A.positions() * units::ANGSTROM_TO_BOHR; + return B.electric_field(pos_A_bohr); +} + +const std::vector &CEEnergyModel::partial_charges() const { + if (m_partial_charges.empty()) { + m_partial_charges.reserve(m_wavefunctions_a.size()); + for (const auto &wfn : m_wavefunctions_a) { + m_partial_charges.push_back(occ::qm::chelpg_charges(wfn)); + } + } + return m_partial_charges; +} + +double CEEnergyModel::coulomb_scale_factor() const { + auto model = ce_model_from_string(m_model_name); + return model.coulomb; +} + +} // namespace occ::interaction diff --git a/src/interaction/lattice_energy.cpp b/src/interaction/lattice_energy.cpp new file mode 100644 index 000000000..f722c5e6e --- /dev/null +++ b/src/interaction/lattice_energy.cpp @@ -0,0 +1,171 @@ +#include +#include +#include +#include +#include +#include +#include + +namespace occ::interaction { + +LatticeEnergyCalculator::LatticeEnergyCalculator( + std::unique_ptr model, const crystal::Crystal &crystal, + const std::string &basename, LatticeConvergenceSettings settings) + : m_energy_model(std::move(model)), m_crystal(crystal), + m_basename(basename), m_settings(settings) { + + if (m_settings.wolf_sum) { + initialize_wolf_sum(); + } +} + +void LatticeEnergyCalculator::initialize_wolf_sum() { + WolfParameters params; + params.cutoff = m_settings.max_radius; + params.alpha = 0.2; // Default alpha value + + m_wolf_sum = std::make_unique(params); + m_wolf_sum->initialize(m_crystal, *m_energy_model); +} + +bool LatticeEnergyCalculator::is_converged(double current, + double previous) const { + return std::abs(current - previous) <= m_settings.energy_tolerance; +} + +void LatticeEnergyCalculator::report_convergence_statistics( + const CEEnergyComponents &total) { + + occ::log::debug("Total polarization term: {:.3f}", + total.polarization_kjmol() * 0.5); + occ::log::debug("Total coulomb term: {:.3f}", total.coulomb_kjmol() * 0.5); + occ::log::debug("Total dispersion term: {:.3f}", total.dispersion_kjmol() * 0.5); + occ::log::debug("Total repulsion term: {:.3f}", total.repulsion_kjmol() * 0.5); + occ::log::debug("Total exchange term: {:.3f}", total.exchange_kjmol() * 0.5); +} + +CEEnergyComponents +compute_cycle_energy(const std::vector &energies, + const crystal::CrystalDimers &all_dimers) { + CEEnergyComponents total; + const auto &mol_neighbors = all_dimers.molecule_neighbors; + + // Loop over all molecules and their neighbors + for (const auto &n : mol_neighbors) { + for (const auto &[dimer, unique_idx] : n) { + const auto &e = energies[unique_idx]; + if (e.is_computed) { + total += e; + } + } + } + return total; +} + +LatticeEnergyResult LatticeEnergyCalculator::compute() { + occ::timing::StopWatch sw; + double lattice_energy = 0.0; + double previous_lattice_energy = 0.0; + double current_radius = + std::max(m_settings.radius_increment, m_settings.min_radius); + size_t cycle = 1; + + auto all_dimers = m_crystal.symmetry_unique_dimers(m_settings.max_radius); + const auto &dimers = all_dimers.unique_dimers; + std::vector converged_energies(dimers.size()); + std::vector charge_energies(dimers.size()); + + occ::log::info("Found {} symmetry unique dimers within max radius {:.3f}", + dimers.size(), m_settings.max_radius); + + do { + previous_lattice_energy = lattice_energy; + + // Create energy store for this cycle + PairEnergyStore store{PairEnergyStore::Kind::XYZ, + fmt::format("{}_dimers", m_basename)}; + + // Setup progress tracking + size_t dimers_to_compute = 0; + for (size_t i = 0; i < dimers.size(); i++) { + if (dimers[i].nearest_distance() <= current_radius && + !converged_energies[i].is_computed) { + dimers_to_compute++; + } + } + occ::core::ProgressTracker progress(dimers_to_compute); + + // Compute energies for this radius + size_t current_dimer = 0; + size_t computed_dimers = 0; + + for (const auto &dimer : dimers) { + sw.start(); + const auto &a = dimer.a(); + const auto &b = dimer.b(); + std::string dimer_name = dimer.name(); + + if (dimer.nearest_distance() <= current_radius && + !converged_energies[current_dimer].is_computed) { + + if (!store.load(current_dimer, dimer, + converged_energies[current_dimer])) { + // Compute new energy + converged_energies[current_dimer] = + m_energy_model->compute_energy(dimer); + store.save(current_dimer, dimer, converged_energies[current_dimer]); + + if (m_settings.wolf_sum) { + charge_energies[current_dimer] = + coulomb_interaction_energy_asym_charges( + dimer, m_wolf_sum->asymmetric_charges()); + } + } + + progress.update(computed_dimers, dimers_to_compute, + fmt::format("E[{}|{}]: {}", a.asymmetric_molecule_idx(), + b.asymmetric_molecule_idx(), dimer_name)); + computed_dimers++; + } + current_dimer++; + sw.stop(); + } + + CEEnergyComponents total = + compute_cycle_energy(converged_energies, all_dimers); + + lattice_energy = 0.5 * total.total_kjmol(); + if (m_settings.wolf_sum) { + double wolf_correction = + m_wolf_sum->compute_correction(charge_energies, converged_energies); + lattice_energy = + m_energy_model->coulomb_scale_factor() * wolf_correction + + 0.5 * total.total_kjmol(); + } + + report_convergence_statistics(total); + occ::log::info("Cycle {} lattice energy: {}", cycle, lattice_energy); + + cycle++; + current_radius += m_settings.radius_increment; + } while (!is_converged(lattice_energy, previous_lattice_energy)); + + // Prepare final result + for (size_t i = 0; i < converged_energies.size(); i++) { + const auto &e = converged_energies[i]; + auto &dimer = all_dimers.unique_dimers[i]; + + if (e.is_computed) { + dimer.set_interaction_energy(e.coulomb_kjmol(), "Coulomb"); + dimer.set_interaction_energy(e.exchange_kjmol(), "Exchange"); + dimer.set_interaction_energy(e.repulsion_kjmol(), "Repulsion"); + dimer.set_interaction_energy(e.dispersion_kjmol(), "Dispersion"); + dimer.set_interaction_energy(e.polarization_kjmol(), "Polarization"); + dimer.set_interaction_energy(e.total_kjmol(), "Total"); + } + } + + return {lattice_energy, all_dimers, converged_energies}; +} + +} // namespace occ::interaction diff --git a/src/interaction/pair_energy.cpp b/src/interaction/pair_energy.cpp index b35005e47..b62a33c09 100644 --- a/src/interaction/pair_energy.cpp +++ b/src/interaction/pair_energy.cpp @@ -1,5 +1,4 @@ #include -#include #include #include #include @@ -7,65 +6,23 @@ #include #include #include -#include #include -#include #include #include #include #include -#include -#include -#include -#include -#include - -namespace fs = std::filesystem; namespace occ::interaction { -using occ::core::Dimer; -using occ::core::Molecule; -using occ::crystal::Crystal; using occ::qm::Wavefunction; -using occ::units::BOHR_TO_ANGSTROM; - -Wavefunction load_wavefunction(const std::string &filename) { - namespace fs = std::filesystem; - using occ::util::to_lower; - std::string ext = fs::path(filename).extension().string(); - to_lower(ext); - if (ext == ".fchk") { - occ::log::debug("Loading Gaussian fchk file from {}", filename); - using occ::io::FchkReader; - FchkReader fchk(filename); - return Wavefunction(fchk); - } else if (ext == ".molden" || ext == ".input") { - occ::log::debug("Loading molden file from {}", filename); - using occ::io::MoldenReader; - MoldenReader molden(filename); - occ::log::debug("Wavefunction has {} atoms", molden.atoms().size()); - return Wavefunction(molden); - } else if (ext == ".json") { - occ::log::debug("Loading OCC JSON wavefunction from {}", filename); - occ::io::JsonWavefunctionReader reader(filename); - return reader.wavefunction(); - } else if (ext == ".orca.json") { - occ::log::debug("Loading Orca JSON file from {}", filename); - occ::io::OrcaJSONReader json(filename); - return Wavefunction(json); - } - throw std::runtime_error( - "Unknown file extension when reading wavefunction: " + ext); -} PairEnergy::PairEnergy(const occ::io::OccInput &input) { if (input.pair.source_a == "none" || input.pair.source_b == "none") { throw "Both monomers in a pair energy need a wavefunction set"; } - a.wfn = load_wavefunction(input.pair.source_a); - b.wfn = load_wavefunction(input.pair.source_b); + a.wfn = Wavefunction::load(input.pair.source_a); + b.wfn = Wavefunction::load(input.pair.source_b); a.rotation = input.pair.rotation_a; b.rotation = input.pair.rotation_b; a.translation = input.pair.translation_a * occ::units::ANGSTROM_TO_BOHR; @@ -118,596 +75,4 @@ CEEnergyComponents PairEnergy::compute() { return energy; } -auto calculate_transform(const Wavefunction &wfn, const Molecule &m, - const Crystal &c) { - using occ::crystal::SymmetryOperation; - - Mat3N pos_m = m.positions(); - Mat3 rotation; - Vec3 translation; - - ankerl::unordered_dense::set symops_tested; - - for (int i = 0; i < m.size(); i++) { - int sint = m.asymmetric_unit_symop()(i); - if (symops_tested.contains(sint)) - continue; - symops_tested.insert(sint); - SymmetryOperation symop(sint); - occ::Mat3N positions = wfn.positions() * BOHR_TO_ANGSTROM; - - rotation = - c.unit_cell().direct() * symop.rotation() * c.unit_cell().inverse(); - translation = (m.centroid() - (rotation * positions).rowwise().mean()) / - BOHR_TO_ANGSTROM; - Wavefunction tmp = wfn; - tmp.apply_transformation(rotation, translation); - - occ::Mat3N tmp_t = tmp.positions() * BOHR_TO_ANGSTROM; - double rmsd = (tmp_t - pos_m).norm(); - occ::log::debug("Test transform (symop={}) RMSD = {}\n", symop.to_string(), - rmsd); - if (rmsd < 1e-3) { - occ::log::debug("Symop found: RMSD = {:12.5f}"); - return std::make_pair(rotation, translation); - } - } - throw std::runtime_error( - "Unable to determine symmetry operation to transform wavefunction"); - return std::make_pair(rotation, translation); -} - -bool write_xyz_dimer(const std::string &filename, const Dimer &dimer, - std::optional energies) { - - using occ::core::Element; - auto output = fmt::output_file(filename, fmt::file::WRONLY | O_TRUNC | - fmt::file::CREATE); - const auto &pos = dimer.positions(); - const auto &nums = dimer.atomic_numbers(); - output.print("{}\n", nums.rows()); - if (energies) { - nlohmann::json j = *energies; - output.print("{}", j.dump()); - } - output.print("\n"); - for (size_t i = 0; i < nums.rows(); i++) { - output.print("{:5s} {:12.5f} {:12.5f} {:12.5f}\n", - Element(nums(i)).symbol(), pos(0, i), pos(1, i), pos(2, i)); - } - return true; -} - -using WavefunctionList = std::vector; - -class CEPairEnergyFunctor { -public: - CEPairEnergyFunctor(const Crystal &crystal, - const WavefunctionList &wavefunctions_a, - const WavefunctionList &wavefunctions_b = {}) - : m_crystal(crystal), m_wavefunctions_a(wavefunctions_a), - m_wavefunctions_b(wavefunctions_b) {} - - void set_model_name(const std::string &model_name) { - m_model_name = model_name; - } - - CEEnergyComponents operator()(const Dimer &dimer) { - Molecule mol_A = dimer.a(); - Molecule mol_B = dimer.b(); - Wavefunction A = m_wavefunctions_a[mol_A.asymmetric_molecule_idx()]; - Wavefunction B = (m_wavefunctions_b.size() > 0) - ? m_wavefunctions_a[mol_B.asymmetric_molecule_idx()] - : m_wavefunctions_b[mol_B.asymmetric_molecule_idx()]; - - auto transform_a = calculate_transform(A, mol_A, m_crystal); - A.apply_transformation(transform_a.first, transform_a.second); - - occ::Mat3N pos_A = mol_A.positions(); - occ::Mat3N pos_A_t = A.positions() * BOHR_TO_ANGSTROM; - occ::log::debug("Mol A transformed wavefunction positions RMSD = {}", - (pos_A_t - pos_A).norm()); - - auto transform_b = calculate_transform(B, mol_B, m_crystal); - B.apply_transformation(transform_b.first, transform_b.second); - - Mat3N pos_B = mol_B.positions(); - Mat3N pos_B_t = B.positions() * BOHR_TO_ANGSTROM; - occ::log::debug("Mol B transformed wavefunction positions RMSD = {}", - (pos_B_t - pos_B).norm()); - - auto model = occ::interaction::ce_model_from_string(m_model_name); - - CEModelInteraction interaction(model); - - auto interaction_energy = interaction(A, B); - interaction_energy.is_computed = true; - - occ::log::debug("Finished model energy"); - return interaction_energy; - } - - Mat3N electric_field(const Dimer &dimer) { - Molecule mol_A = dimer.a(); - Molecule mol_B = dimer.b(); - Wavefunction A = m_wavefunctions_a[mol_A.asymmetric_molecule_idx()]; - Wavefunction B = (m_wavefunctions_b.size() > 0) - ? m_wavefunctions_a[mol_B.asymmetric_molecule_idx()] - : m_wavefunctions_b[mol_B.asymmetric_molecule_idx()]; - - auto transform_b = calculate_transform(B, mol_B, m_crystal); - B.apply_transformation(transform_b.first, transform_b.second); - - occ::Mat3N pos_A = mol_A.positions(); - Mat3N pos_A_t = A.positions() * BOHR_TO_ANGSTROM; - - occ::Mat3N pos_A_bohr = mol_A.positions() * occ::units::ANGSTROM_TO_BOHR; - - Mat3N pos_B = mol_B.positions(); - Mat3N pos_B_t = B.positions() * BOHR_TO_ANGSTROM; - - return B.electric_field(pos_A_bohr); - } - - inline const auto &wavefunctions() const { return m_wavefunctions_a; } - - inline const auto &partial_charges() { - if (m_partial_charges.size() > 0) - return m_partial_charges; - - m_partial_charges = std::vector(m_wavefunctions_a.size()); - for (int i = 0; i < m_wavefunctions_a.size(); i++) { - m_partial_charges[i] = occ::qm::chelpg_charges(m_wavefunctions_a[i]); - } - return m_partial_charges; - } - - inline double coulomb_scale_factor() const { - auto model = occ::interaction::ce_model_from_string(m_model_name); - return model.coulomb; - } - -private: - Crystal m_crystal; - std::string m_model_name{"ce-b3lyp"}; - std::vector m_wavefunctions_a; - std::vector m_wavefunctions_b; - std::vector m_partial_charges; -}; - -class XTBPairEnergyFunctor { -public: - XTBPairEnergyFunctor(const Crystal &crystal) : m_crystal(crystal) { - for (const auto &mol : crystal.symmetry_unique_molecules()) { - occ::xtb::XTBCalculator calc(mol); - m_monomer_energies.push_back(calc.single_point_energy()); - - m_partial_charges.push_back(calc.partial_charges()); - } - } - - CEEnergyComponents operator()(const Dimer &dimer) { - Molecule mol_A = dimer.a(); - Molecule mol_B = dimer.b(); - - occ::xtb::XTBCalculator calc_AB(dimer); - double e_a = m_monomer_energies[mol_A.asymmetric_molecule_idx()]; - double e_b = m_monomer_energies[mol_B.asymmetric_molecule_idx()]; - double e_ab = calc_AB.single_point_energy(); - CEEnergyComponents result; - result.total = e_ab - e_a - e_b; - result.is_computed = true; - return result; - } - - inline const auto &partial_charges() { return m_partial_charges; } - inline const auto &monomer_energies() { return m_monomer_energies; } - inline double coulomb_scale_factor() const { return 1.0; } - -private: - Crystal m_crystal; - std::vector m_monomer_energies; - std::vector m_partial_charges; -}; - -int compute_coulomb_energies_radius(const std::vector &dimers, - const Vec &asym_charges, double radius, - std::vector &charge_energies) { - occ::timing::StopWatch sw; - if (charge_energies.size() < dimers.size()) { - charge_energies.resize(dimers.size()); - } - - size_t current_dimer{0}; - size_t computed_dimers{0}; - for (const auto &dimer : dimers) { - if (dimer.nearest_distance() > radius) { - current_dimer++; - continue; - } - sw.start(); - - charge_energies[current_dimer] = - occ::interaction::coulomb_interaction_energy_asym_charges(dimer, - asym_charges); - current_dimer++; - } - occ::log::debug("Finished calculating {} unique dimer coulomb energies", - computed_dimers); - return computed_dimers; -} - -Mat3N compute_point_charge_efield(const Dimer &dimer, const Vec &asym_charges) { - return occ::interaction::coulomb_efield_asym_charges(dimer, asym_charges) - .first; -} - -template -int compute_dimer_energies_radius( - EnergyModel &energy_model, const Crystal &crystal, - const std::vector &dimers, const std::string &basename, - double radius, std::vector &dimer_energies) { - - using occ::crystal::SymmetryOperation; - occ::timing::StopWatch sw; - if (dimer_energies.size() < dimers.size()) { - dimer_energies.resize(dimers.size()); - } - - size_t current_dimer{0}; - size_t computed_dimers{0}; - size_t dimers_to_compute{0}; - for (size_t i = 0; i < dimers.size(); i++) { - const auto &dimer = dimers[i]; - if (dimer.nearest_distance() > radius || dimer_energies[i].is_computed) - continue; - dimers_to_compute++; - } - - PairEnergyStore store{PairEnergyStore::Kind::Xyz, - fmt::format("{}_dimers", basename)}; - occ::core::ProgressTracker progress(dimers_to_compute); - - for (const auto &dimer : dimers) { - auto tprev = sw.read(); - sw.start(); - - const auto &a = dimer.a(); - const auto &b = dimer.b(); - const auto asym_idx_a = a.asymmetric_molecule_idx(); - const auto asym_idx_b = b.asymmetric_molecule_idx(); - const auto shift_b = dimer.b().cell_shift(); - std::string dimer_name = dimer.name(); - - CEEnergyComponents &dimer_energy = dimer_energies[current_dimer]; - std::string dimer_energy_file( - fmt::format("{}_dimer_{}_energies.xyz", basename, current_dimer)); - - if (dimer.nearest_distance() > radius || dimer_energy.is_computed) { - current_dimer++; - continue; - } - if (store.load(current_dimer, dimer, dimer_energy)) { - progress.update( - computed_dimers, dimers_to_compute, - fmt::format("Load [{}|{}]: {}", asym_idx_a, asym_idx_b, dimer_name)); - computed_dimers++; - current_dimer++; - continue; - } - - progress.update( - computed_dimers, dimers_to_compute, - fmt::format("E[{}|{}]: {}", asym_idx_a, asym_idx_b, dimer_name)); - occ::log::debug( - "{} ({}[{}] - {}[{} + ({},{},{})]), Rc = {: 5.2f}", current_dimer, - asym_idx_a, SymmetryOperation(a.asymmetric_unit_symop()(0)).to_string(), - asym_idx_b, SymmetryOperation(b.asymmetric_unit_symop()(0)).to_string(), - shift_b[0], shift_b[1], shift_b[2], dimer.center_of_mass_distance()); - - std::cout << std::flush; - dimer_energy = energy_model(dimer); - sw.stop(); - occ::log::debug("Took {:.3f} seconds", sw.read() - tprev); - store.save(current_dimer, dimer, dimer_energy); - computed_dimers++; - current_dimer++; - } - progress.clear(); - occ::log::info( - "Finished calculating {} unique dimer interaction energies in " - "{:%H:%M:%S}", - computed_dimers, - std::chrono::round(progress.time_taken())); - return computed_dimers; -} - -std::vector -ce_model_energies(const Crystal &crystal, const std::vector &dimers, - const std::vector &wfns_a, - const std::vector &wfns_b, - const std::string &basename) { - std::vector result; - CEPairEnergyFunctor f(crystal, wfns_a, wfns_b); - compute_dimer_energies_radius(f, crystal, dimers, basename, - std::numeric_limits::max(), result); - return result; -} - -bool load_dimer_energy(const std::string &filename, - CEEnergyComponents &energies) { - if (!fs::exists(filename)) - return false; - occ::log::debug("Load dimer energies from {}", filename); - std::ifstream file(filename); - std::string line; - std::getline(file, line); - std::getline(file, line); - energies = nlohmann::json::parse(line); - energies.is_computed = true; - return true; -} - -struct WolfSumAccelerator { - occ::interaction::WolfParams parameters{16.0, 0.2}; - Vec asym_charges; - std::vector charge_self_energies; - std::vector electric_field_values; - double energy{0.0}; - - template - void setup(const Crystal &crystal, EnergyModel &energy_model) { - auto asym_mols = crystal.symmetry_unique_molecules(); - asym_charges = Vec(crystal.asymmetric_unit().size()); - charge_self_energies = std::vector(asym_mols.size()); - // vector per unique molecule or wfn. - - const auto &partial_charges = energy_model.partial_charges(); - for (int i = 0; i < partial_charges.size(); i++) { - const auto &asymmetric_atom_indices = asym_mols[i].asymmetric_unit_idx(); - const auto &charge_vector = partial_charges[i]; - occ::log::info("Charges used in wolf for molecule {}", i); - for (int j = 0; j < charge_vector.rows(); j++) { - occ::log::info("Atom {}: {:12.5f}", j, charge_vector(j)); - asym_charges(asymmetric_atom_indices(j)) = charge_vector(j); - } - } - - int asym_idx = 0; - auto surrounds = - crystal.asymmetric_unit_atom_surroundings(parameters.cutoff); - Mat3N asym_cart = crystal.to_cartesian(crystal.asymmetric_unit().positions); - Vec asym_wolf(surrounds.size()); - for (const auto &s : surrounds) { - double qi = asym_charges(asym_idx); - Vec3 pi = asym_cart.col(asym_idx); - Vec qj(s.size()); - for (int j = 0; j < qj.rows(); j++) { - qj(j) = asym_charges(s.asym_idx(j)); - } - asym_wolf(asym_idx) = occ::interaction::wolf_coulomb_energy( - qi, pi, qj, s.cart_pos, parameters) * - units::AU_TO_KJ_PER_MOL; - asym_idx++; - } - for (int i = 0; i < asym_mols.size(); i++) { - const auto &mol = asym_mols[i]; - electric_field_values.push_back(Mat3N::Zero(3, mol.size())); - charge_self_energies[i] = - occ::interaction::coulomb_self_energy_asym_charges(mol, asym_charges); - for (int j = 0; j < mol.size(); j++) { - energy += asym_wolf(mol.asymmetric_unit_idx()(j)); - } - } - - occ::log::debug("Wolf energy ({} asymmetric atoms): {}\n", asym_idx, - asym_wolf.sum()); - - occ::log::debug("Wolf energy ({} asymmetric molecules): {}\n", - asym_mols.size(), energy); - } - - void initialize(const Crystal &crystal) {} -}; - -template -LatticeEnergyResult -converged_lattice_energies(EnergyModel &energy_model, const Crystal &crystal, - const std::string &basename, - const LatticeConvergenceSettings conv) { - - crystal::CrystalDimers converged_dimers; - std::vector converged_energies; - std::vector charge_energies; - double lattice_energy{0.0}, previous_lattice_energy{0.0}; - double current_radius = std::max(conv.radius_increment, conv.min_radius); - size_t cycle{1}; - - auto all_dimers = crystal.symmetry_unique_dimers(conv.max_radius); - const auto &asym_mols = crystal.symmetry_unique_molecules(); - - occ::log::info("Found {} symmetry unique dimers within max radius {:.3f}\n", - all_dimers.unique_dimers.size(), conv.max_radius); - occ::log::info("Lattice convergence settings:"); - occ::log::info("Start radius {: 8.4f} angstrom", conv.min_radius); - occ::log::info("Max radius {: 8.4f} angstrom", conv.max_radius); - occ::log::info("Radius increment {: 8.4f} angstrom", conv.radius_increment); - occ::log::info("Energy tolerance {: 8.4f} kJ/mol", conv.energy_tolerance); - - WolfSumAccelerator wolf; - if (conv.wolf_sum) { - wolf.parameters.cutoff = conv.max_radius; - wolf.setup(crystal, energy_model); - } - - do { - previous_lattice_energy = lattice_energy; - const auto &dimers = all_dimers.unique_dimers; - compute_dimer_energies_radius(energy_model, crystal, dimers, basename, - current_radius, converged_energies); - if (conv.wolf_sum) { - compute_coulomb_energies_radius(dimers, wolf.asym_charges, current_radius, - charge_energies); - } - - const auto &mol_neighbors = all_dimers.molecule_neighbors; - CEEnergyComponents total; - size_t mol_idx{0}; - double epol_total{0}; - double ecoul_real{0}; - double ecoul_exact_real{0}; - double ecoul_self{0}; - double coulomb_scale_factor = energy_model.coulomb_scale_factor(); - for (const auto &n : mol_neighbors) { - CEEnergyComponents molecule_total; - size_t dimer_idx{0}; - - Mat3N efield = Mat3N::Zero(3, asym_mols[mol_idx].size()); - - if (conv.wolf_sum && conv.crystal_field_polarization) { - wolf.electric_field_values[mol_idx].setZero(); - occ::log::debug("Field total for dimer {} =\n{}\n", dimer_idx, - format_matrix(wolf.electric_field_values[mol_idx])); - } - - for (const auto &[dimer, unique_idx] : n) { - const auto &e = converged_energies[unique_idx]; - - if (e.is_computed) { - molecule_total += e; - epol_total += e.polarization_kjmol(); - if (conv.wolf_sum) { - double e_charge = charge_energies[unique_idx]; - ecoul_exact_real += 0.5 * e.coulomb_kjmol(); - ecoul_real += 0.5 * e_charge * occ::units::AU_TO_KJ_PER_MOL; - if (conv.crystal_field_polarization) { - wolf.electric_field_values[mol_idx] += - compute_point_charge_efield(dimer, wolf.asym_charges); - } - } - if (conv.crystal_field_polarization) { - if constexpr (std::is_same::value) { - efield += energy_model.electric_field(dimer); - } - } - } - dimer_idx++; - } - occ::log::debug("epol total = {}", epol_total); - total += molecule_total; - if (conv.wolf_sum) { - ecoul_self += - wolf.charge_self_energies[mol_idx] * units::AU_TO_KJ_PER_MOL; - - if (conv.crystal_field_polarization) { - auto &electric_field = wolf.electric_field_values[mol_idx]; - occ::log::debug("Field total =\n{}", format_matrix(electric_field)); - occ::log::debug("Field total =\n{}", format_matrix(efield)); - if constexpr (std::is_same::value) { - const auto &wfn_a = energy_model.wavefunctions()[mol_idx]; - double e_pol_chg = occ::interaction::polarization_energy( - wfn_a.xdm_polarizabilities, electric_field); - occ::log::debug("Crystal polarizability (chg): {}", - e_pol_chg * occ::units::AU_TO_KJ_PER_MOL); - } - } - } - if (conv.crystal_field_polarization) { - if constexpr (std::is_same::value) { - const auto &wfn_a = energy_model.wavefunctions()[mol_idx]; - double e_pol_qm = occ::interaction::polarization_energy( - wfn_a.xdm_polarizabilities, efield); - occ::log::debug("Crystal polarizability (qm): {}", - e_pol_qm * occ::units::AU_TO_KJ_PER_MOL); - } - } - mol_idx++; - } - lattice_energy = 0.5 * total.total_kjmol(); - if (conv.wolf_sum) { - occ::log::info("Charge-charge intramolecular: {}", ecoul_self); - occ::log::info("Charge-charge real space: {}", ecoul_real); - occ::log::info("Wolf energy: {}", wolf.energy); - occ::log::info("Coulomb (exact) real: {}", ecoul_exact_real); - occ::log::info("Wolf - intra: {}", wolf.energy - ecoul_self); - occ::log::info("Wolf corrected Coulomb total: {}", - wolf.energy - ecoul_self - ecoul_real + ecoul_exact_real); - lattice_energy = - coulomb_scale_factor * (wolf.energy - ecoul_self - ecoul_real) + - 0.5 * total.total_kjmol(); - } - occ::log::info("Cycle {} lattice energy: {}", cycle, lattice_energy); - occ::log::debug("Total polarization term: {:.3f}", - 0.5 * total.polarization_kjmol()); - occ::log::debug("Total coulomb term: {:.3f}", 0.5 * total.coulomb_kjmol()); - occ::log::debug("Total dispersion term: {:.3f}", - 0.5 * total.dispersion_kjmol()); - occ::log::debug("Total repulsion term: {:.3f}", - 0.5 * total.repulsion_kjmol()); - occ::log::debug("Total exchange term: {:.3f}", - 0.5 * total.exchange_kjmol()); - occ::log::debug("Wolf correction: {:.3f}", - (wolf.energy - ecoul_self - ecoul_real)); - cycle++; - current_radius += conv.radius_increment; - } while (std::abs(lattice_energy - previous_lattice_energy) > - conv.energy_tolerance); - converged_dimers = all_dimers; - for (int i = 0; i < converged_energies.size(); i++) { - if (i > converged_dimers.unique_dimers.size()) { - throw std::runtime_error("Invalid dimer index for unique dimers when " - "propagating dimer energies"); - } - auto &dimer = converged_dimers.unique_dimers[i]; - const auto &e = converged_energies[i]; - dimer.set_interaction_energy(e.coulomb_kjmol(), "Coulomb"); - dimer.set_interaction_energy(e.exchange_kjmol(), "Exchange"); - dimer.set_interaction_energy(e.repulsion_kjmol(), "Repulsion"); - dimer.set_interaction_energy(e.dispersion_kjmol(), "Dispersion"); - dimer.set_interaction_energy(e.polarization_kjmol(), "Polarization"); - dimer.set_interaction_energy(e.total_kjmol(), "Total"); - } - return {lattice_energy, converged_dimers, converged_energies}; -} - -LatticeEnergyResult converged_lattice_energies( - const Crystal &crystal, const std::vector &wfns_a, - const std::vector &wfns_b, const std::string &basename, - const LatticeConvergenceSettings conv) { - - CEPairEnergyFunctor energy_model(crystal, wfns_a, wfns_b); - energy_model.set_model_name(conv.model_name); - return converged_lattice_energies(energy_model, crystal, basename, conv); -} - -LatticeEnergyResult -converged_xtb_lattice_energies(const Crystal &crystal, - const std::string &basename, - const LatticeConvergenceSettings conv) { - - XTBPairEnergyFunctor energy_model(crystal); - return converged_lattice_energies(energy_model, crystal, basename, conv); -} - -std::string PairEnergyStore::dimer_filename(int id, const Dimer &d) { - return fmt::format("dimer_{}.xyz", id); -} - -bool PairEnergyStore::save(int id, const Dimer &d, - const CEEnergyComponents &e) { - fs::path parent(name); - if (!fs::exists(parent)) { - fs::create_directories(parent); - } - return write_xyz_dimer((parent / fs::path(dimer_filename(id, d))).string(), d, - e); -} - -bool PairEnergyStore::load(int id, const Dimer &d, CEEnergyComponents &e) { - fs::path parent(name); - return load_dimer_energy((parent / fs::path(dimer_filename(id, d))).string(), - e); -} - } // namespace occ::interaction diff --git a/src/interaction/pair_energy_store.cpp b/src/interaction/pair_energy_store.cpp new file mode 100644 index 000000000..892e0e1d4 --- /dev/null +++ b/src/interaction/pair_energy_store.cpp @@ -0,0 +1,148 @@ +#include +#include +#include +#include +#include + +namespace fs = std::filesystem; +using occ::core::Dimer; + +namespace occ::interaction { + +bool load_dimer_energy(const std::string &filename, + CEEnergyComponents &energies) { + if (!fs::exists(filename)) + return false; + occ::log::debug("Load dimer energies from {}", filename); + std::ifstream file(filename); + std::string line; + std::getline(file, line); + std::getline(file, line); + energies = nlohmann::json::parse(line); + energies.is_computed = true; + return true; +} + +bool write_xyz_dimer(const std::string &filename, const Dimer &dimer, + std::optional energies) { + + using occ::core::Element; + auto output = fmt::output_file(filename, fmt::file::WRONLY | O_TRUNC | + fmt::file::CREATE); + const auto &pos = dimer.positions(); + const auto &nums = dimer.atomic_numbers(); + output.print("{}\n", nums.rows()); + if (energies) { + nlohmann::json j = *energies; + output.print("{}", j.dump()); + } + output.print("\n"); + for (size_t i = 0; i < nums.rows(); i++) { + output.print("{:5s} {:12.5f} {:12.5f} {:12.5f}\n", + Element(nums(i)).symbol(), pos(0, i), pos(1, i), pos(2, i)); + } + return true; +} + +bool MemoryPairEnergyStore::save(int id, const core::Dimer &, + const CEEnergyComponents &e) { + energy_store[id] = e; + return true; +} + +bool MemoryPairEnergyStore::load(int id, const core::Dimer &, + CEEnergyComponents &e) { + auto it = energy_store.find(id); + if (it == energy_store.end()) + return false; + e = it->second; + e.is_computed = true; + return true; +} + +std::string MemoryPairEnergyStore::dimer_filename(int id, + const core::Dimer &) const { + return fmt::format("memory_dimer_{}", id); +} + +bool FileSystemPairEnergyStore::save(int id, const core::Dimer &d, + const CEEnergyComponents &e) { + fs::path parent(base_path); + if (!fs::exists(parent)) { + fs::create_directories(parent); + } + + auto filepath = parent / fs::path(dimer_filename(id, d)); + + if (format == Format::JSON) { + nlohmann::json j; + j["id"] = id; + j["energy"] = e; + std::ofstream file(filepath); + file << j.dump(2); + return true; + } else { + return write_xyz_dimer(filepath.string(), d, e); + } +} + +bool FileSystemPairEnergyStore::load(int id, const core::Dimer &d, + CEEnergyComponents &e) { + fs::path parent(base_path); + auto filepath = parent / fs::path(dimer_filename(id, d)); + + if (!fs::exists(filepath)) + return false; + + if (format == Format::JSON) { + std::ifstream file(filepath); + nlohmann::json j; + file >> j; + e = j["energy"].get(); + e.is_computed = true; + return true; + } else { + return load_dimer_energy(filepath.string(), e); + } +} + +std::string +FileSystemPairEnergyStore::dimer_filename(int id, const core::Dimer &) const { + return fmt::format("dimer_{}.{}", id, + format == Format::JSON ? "json" : "xyz"); +} + +// PairEnergyStore implementation +PairEnergyStore::PairEnergyStore(Kind kind, std::string name) + : kind(kind), name(std::move(name)) { + switch (kind) { + case Kind::Memory: + store = std::make_unique(); + break; + case Kind::JSON: + store = std::make_unique( + this->name, FileSystemPairEnergyStore::Format::JSON); + break; + case Kind::XYZ: + store = std::make_unique( + this->name, FileSystemPairEnergyStore::Format::XYZ); + break; + } +} + +bool PairEnergyStore::save(int id, const core::Dimer &d, + const CEEnergyComponents &e) { + return store->save(id, d, e); +} + +bool PairEnergyStore::load(int id, const core::Dimer &d, + CEEnergyComponents &e) { + return store->load(id, d, e); +} + +std::string PairEnergyStore::dimer_filename(int id, + const core::Dimer &d) const { + return store->dimer_filename(id, d); +} + +} // namespace occ::interaction diff --git a/src/interaction/wavefunction_transform.cpp b/src/interaction/wavefunction_transform.cpp new file mode 100644 index 000000000..155e1b428 --- /dev/null +++ b/src/interaction/wavefunction_transform.cpp @@ -0,0 +1,59 @@ +#include +#include +#include + +namespace occ::interaction::transform { + +double +WavefunctionTransformer::compute_position_rmsd(const Mat3N &pos_ref, + const Mat3N &pos_transformed) { + return (pos_transformed - pos_ref).norm(); +} + +TransformResult +WavefunctionTransformer::calculate_transform(const qm::Wavefunction &wfn, + const core::Molecule &mol, + const crystal::Crystal &crystal) { + + using occ::crystal::SymmetryOperation; + TransformResult result; + + Mat3N pos_mol = mol.positions(); + ankerl::unordered_dense::set symops_tested; + + for (int i = 0; i < mol.size(); i++) { + int sint = mol.asymmetric_unit_symop()(i); + if (symops_tested.contains(sint)) + continue; + + symops_tested.insert(sint); + SymmetryOperation symop(sint); + Mat3N positions = wfn.positions() * units::BOHR_TO_ANGSTROM; + + result.rotation = crystal.unit_cell().direct() * symop.rotation() * + crystal.unit_cell().inverse(); + + result.translation = + (mol.centroid() - (result.rotation * positions).rowwise().mean()) / + units::BOHR_TO_ANGSTROM; + + qm::Wavefunction tmp = wfn; + tmp.apply_transformation(result.rotation, result.translation); + + Mat3N transformed_pos = tmp.positions() * units::BOHR_TO_ANGSTROM; + result.rmsd = compute_position_rmsd(pos_mol, transformed_pos); + + occ::log::debug("Test transform (symop={}) RMSD = {}\n", symop.to_string(), + result.rmsd); + + if (result.rmsd < 1e-3) { + occ::log::debug("Symop found: RMSD = {:12.5f}", result.rmsd); + return result; + } + } + + throw std::runtime_error( + "Unable to determine symmetry operation to transform wavefunction"); +} + +} // namespace occ::interaction::transform diff --git a/src/interaction/wolf.cpp b/src/interaction/wolf.cpp index fc8262b33..83c8a4f22 100644 --- a/src/interaction/wolf.cpp +++ b/src/interaction/wolf.cpp @@ -1,6 +1,8 @@ #include #include +#include #include +#include #include #include @@ -8,10 +10,10 @@ namespace occ::interaction { double wolf_coulomb_energy(double qi, const Vec3 &pi, Eigen::Ref qj, Eigen::Ref pj, - const WolfParams ¶ms) { + const WolfParameters ¶ms) { using occ::constants::sqrt_pi; using std::erfc; - double eta = params.eta / occ::units::ANGSTROM_TO_BOHR; + double eta = params.alpha / occ::units::ANGSTROM_TO_BOHR; double rc = params.cutoff * occ::units::ANGSTROM_TO_BOHR; double trc = erfc(eta * rc) / rc; @@ -24,4 +26,94 @@ double wolf_coulomb_energy(double qi, const Vec3 &pi, Eigen::Ref qj, return 0.5 * pair_term - self_term; } +WolfSum::WolfSum(WolfParameters params) : m_params(params) {} + +void WolfSum::initialize(const crystal::Crystal &crystal, + const EnergyModelBase &energy_model) { + + auto asym_mols = crystal.symmetry_unique_molecules(); + m_asym_charges = Vec(crystal.asymmetric_unit().size()); + m_charge_self_energies.resize(asym_mols.size()); + m_electric_field_values.clear(); + m_total_energy = 0.0; + + const auto &partial_charges = energy_model.partial_charges(); + for (size_t i = 0; i < partial_charges.size(); i++) { + const auto &asymmetric_atom_indices = asym_mols[i].asymmetric_unit_idx(); + const auto &charge_vector = partial_charges[i]; + + occ::log::info("Charges used in wolf for molecule {}", i); + for (int j = 0; j < charge_vector.rows(); j++) { + occ::log::info("Atom {}: {:12.5f}", j, charge_vector(j)); + m_asym_charges(asymmetric_atom_indices(j)) = charge_vector(j); + } + } + + compute_self_energies(crystal); + compute_wolf_energies(crystal); +} + +void WolfSum::compute_self_energies(const crystal::Crystal &crystal) { + auto asym_mols = crystal.symmetry_unique_molecules(); + for (size_t i = 0; i < asym_mols.size(); i++) { + const auto &mol = asym_mols[i]; + m_electric_field_values.push_back(Mat3N::Zero(3, mol.size())); + m_charge_self_energies[i] = + coulomb_self_energy_asym_charges(mol, m_asym_charges); + } +} + +void WolfSum::compute_wolf_energies(const crystal::Crystal &crystal) { + auto surrounds = crystal.asymmetric_unit_atom_surroundings(m_params.cutoff); + Mat3N asym_cart = crystal.to_cartesian(crystal.asymmetric_unit().positions); + Vec asym_wolf(surrounds.size()); + + for (size_t asym_idx = 0; asym_idx < surrounds.size(); asym_idx++) { + const auto &s = surrounds[asym_idx]; + double qi = m_asym_charges(asym_idx); + Vec3 pi = asym_cart.col(asym_idx); + + Vec qj(s.size()); + for (int j = 0; j < qj.rows(); j++) { + qj(j) = m_asym_charges(s.asym_idx(j)); + } + + asym_wolf(asym_idx) = + wolf_coulomb_energy(qi, pi, qj, s.cart_pos, m_params) * + units::AU_TO_KJ_PER_MOL; + } + + auto asym_mols = crystal.symmetry_unique_molecules(); + for (const auto &mol : asym_mols) { + for (int j = 0; j < mol.size(); j++) { + m_total_energy += asym_wolf(mol.asymmetric_unit_idx()(j)); + } + } + + occ::log::debug("Wolf energy ({} asymmetric molecules): {}\n", + asym_mols.size(), m_total_energy); +} + +double WolfSum::compute_correction( + const std::vector &charge_energies, + const std::vector &model_energies) const { + + double ecoul_real = 0.0; + double ecoul_exact_real = 0.0; + double ecoul_self = 0.0; + + for (size_t i = 0; i < charge_energies.size(); i++) { + if (model_energies[i].is_computed) { + ecoul_exact_real += 0.5 * model_energies[i].coulomb_kjmol(); + ecoul_real += 0.5 * charge_energies[i] * units::AU_TO_KJ_PER_MOL; + } + } + + for (double self_energy : m_charge_self_energies) { + ecoul_self += self_energy * units::AU_TO_KJ_PER_MOL; + } + + return (m_total_energy - ecoul_self - ecoul_real + ecoul_exact_real); +} + } // namespace occ::interaction diff --git a/src/interaction/xtb_energy_model.cpp b/src/interaction/xtb_energy_model.cpp new file mode 100644 index 000000000..0cf90ed74 --- /dev/null +++ b/src/interaction/xtb_energy_model.cpp @@ -0,0 +1,40 @@ +#include +#include + +namespace occ::interaction { + +XTBEnergyModel::XTBEnergyModel(const crystal::Crystal &crystal) + : m_crystal(crystal) { + + for (const auto &mol : crystal.symmetry_unique_molecules()) { + occ::xtb::XTBCalculator calc(mol); + m_monomer_energies.push_back(calc.single_point_energy()); + m_partial_charges.push_back(calc.partial_charges()); + } +} + +CEEnergyComponents XTBEnergyModel::compute_energy(const core::Dimer &dimer) { + core::Molecule mol_A = dimer.a(); + core::Molecule mol_B = dimer.b(); + + occ::xtb::XTBCalculator calc_AB(dimer); + double e_a = m_monomer_energies[mol_A.asymmetric_molecule_idx()]; + double e_b = m_monomer_energies[mol_B.asymmetric_molecule_idx()]; + double e_ab = calc_AB.single_point_energy(); + + CEEnergyComponents result; + result.total = e_ab - e_a - e_b; + result.is_computed = true; + return result; +} + +Mat3N XTBEnergyModel::compute_electric_field(const core::Dimer &) { + // XTB doesn't provide electric field calculations + return Mat3N::Zero(3, 1); +} + +const std::vector &XTBEnergyModel::partial_charges() const { + return m_partial_charges; +} + +} // namespace occ::interaction diff --git a/src/main/make_atomic_pair_potentials.cpp b/src/main/make_atomic_pair_potentials.cpp index 0ec764a49..6ca0de324 100644 --- a/src/main/make_atomic_pair_potentials.cpp +++ b/src/main/make_atomic_pair_potentials.cpp @@ -80,7 +80,7 @@ CEEnergyComponents compute_pair_energy(Wavefunction wfn_a, Wavefunction wfn_b, int main(int argc, char *argv[]) { occ::timing::start(occ::timing::category::global); occ::timing::start(occ::timing::category::io); - occ::log::setup_logging(2); + occ::log::set_log_level(2); std::string symbol_a, symbol_b; int min_charge_cation{0}, min_charge_anion{0}; int max_charge_cation{0}, max_charge_anion{0}; @@ -116,7 +116,7 @@ int main(int argc, char *argv[]) { // logging verbosity auto *verbosity_option = app.add_flag_function( "--verbosity{2}", - [](int verbosity) { occ::log::setup_logging(verbosity); }, + [](int verbosity) { occ::log::set_log_level(verbosity); }, "logging verbosity {0=silent,1=minimal,2=normal,3=verbose,4=debug}"); verbosity_option->default_val(2); verbosity_option->run_callback_for_default(); diff --git a/src/main/occ.cpp b/src/main/occ.cpp index 39eec8864..c6c3f574e 100644 --- a/src/main/occ.cpp +++ b/src/main/occ.cpp @@ -16,7 +16,7 @@ int main(int argc, char *argv[]) { occ::timing::start(occ::timing::category::global); occ::timing::start(occ::timing::category::io); - occ::log::setup_logging(2); + occ::log::set_log_level(2); CLI::App app("occ - A program for quantum chemistry"); app.allow_config_extras(CLI::config_extras_mode::error); @@ -38,7 +38,7 @@ int main(int argc, char *argv[]) { // logging verbosity auto *verbosity_option = app.add_flag_function( "--verbosity{2}", - [](int verbosity) { occ::log::setup_logging(verbosity); }, + [](int verbosity) { occ::log::set_log_level(verbosity); }, "logging verbosity {0=silent,1=minimal,2=normal,3=verbose,4=debug}"); verbosity_option->default_val(2); verbosity_option->run_callback_for_default(); diff --git a/src/main/occ_elat.cpp b/src/main/occ_elat.cpp index 247d2aef1..e9aa486f3 100644 --- a/src/main/occ_elat.cpp +++ b/src/main/occ_elat.cpp @@ -7,10 +7,14 @@ #include #include #include +#include #include +#include #include +#include #include #include +#include #include #include #include @@ -21,8 +25,13 @@ namespace fs = std::filesystem; using occ::crystal::Crystal; +using occ::crystal::SymmetryDimerLabeller; using occ::interaction::CEEnergyComponents; +using occ::interaction::CEEnergyModel; using occ::interaction::LatticeConvergenceSettings; +using occ::interaction::LatticeEnergyCalculator; +using occ::interaction::LatticeEnergyResult; +using occ::interaction::XTBEnergyModel; using occ::qm::Wavefunction; inline Crystal read_crystal(const std::string &filename) { @@ -104,7 +113,6 @@ inline void write_elat_json(const std::string &basename, const std::string &model, const occ::crystal::Crystal &crystal, const occ::crystal::CrystalDimers &dimers) { - nlohmann::json j; j["result_type"] = "elat"; j["title"] = basename; @@ -113,23 +121,39 @@ inline void write_elat_json(const std::string &basename, j["has_permutation_symmetry"] = true; const auto &uc_atoms = crystal.unit_cell_atoms(); + auto dimer_labeller = SymmetryDimerLabeller(crystal); + dimer_labeller.connection = "-"; + dimer_labeller.format.fmt_string = "{}"; j["pairs"] = {}; for (const auto &mol_pairs : dimers.molecule_neighbors) { nlohmann::json m; for (const auto &[dimer, unique_idx] : mol_pairs) { + const auto &unique_dimer = dimers.unique_dimers[unique_idx]; + if (unique_dimer.interaction_energy() == 0.0) + continue; + nlohmann::json d; nlohmann::json e; - const auto &unique_dimer = dimers.unique_dimers[unique_idx]; + + // Label generation + auto label = dimer_labeller(dimer); + d["Label"] = label; + d["Unique Index"] = unique_idx; + + // Energy components const auto &energies = unique_dimer.interaction_energies(); - if (energies.at("total") == 0.0) - continue; - e["unique_dimer_index"] = unique_idx; for (const auto &[k, v] : energies) { e[k] = v; } d["energies"] = e; + // Nearest neighbor calculation based on distance threshold + bool is_nearest = dimer.nearest_distance() <= + 4.0; // You may want to make this threshold configurable + d["Nearest Neighbor"] = is_nearest; + + // Unit cell atom offsets nlohmann::json offsets_a = {}; { const auto &a = dimer.a(); @@ -141,6 +165,7 @@ inline void write_elat_json(const std::string &basename, a_uc_shift(2, i)}); } } + nlohmann::json offsets_b = {}; { const auto &b = dimer.b(); @@ -153,10 +178,12 @@ inline void write_elat_json(const std::string &basename, } } d["uc_atom_offsets"] = {offsets_a, offsets_b}; + m.push_back(d); } j["pairs"].push_back(m); } + std::ofstream dest(fmt::format("{}_elat_results.json", basename)); dest << j.dump(2); } @@ -199,18 +226,26 @@ void calculate_lattice_energy(const LatticeConvergenceSettings settings) { occ::log::info("Calculating symmetry unique dimers"); occ::crystal::CrystalDimers crystal_dimers; std::vector energies; - occ::interaction::LatticeEnergyResult lattice_energy_result; + + std::unique_ptr energy_model; + if (settings.model_name == "xtb") { - lattice_energy_result = - converged_xtb_lattice_energies(c, basename, settings); + energy_model = std::make_unique(c); } else { wfns = occ::main::calculate_wavefunctions( basename, molecules, settings.model_name, settings.spherical_basis); occ::main::compute_monomer_energies(basename, wfns, settings.model_name); - lattice_energy_result = occ::interaction::converged_lattice_energies( - c, wfns, wfns, basename, settings); + + auto ce_model = std::make_unique(c, wfns, wfns); + ce_model->set_model_name(settings.model_name); + energy_model = std::move(ce_model); } + LatticeEnergyCalculator calculator(std::move(energy_model), c, basename, + settings); + + LatticeEnergyResult lattice_energy_result = calculator.compute(); + const auto &dimers = lattice_energy_result.dimers.unique_dimers; if (dimers.size() < 1) { occ::log::error("No dimers found using neighbour radius {:.3f}", diff --git a/src/occpy.cpp b/src/occpy.cpp index 471a007c0..98ac91f01 100644 --- a/src/occpy.cpp +++ b/src/occpy.cpp @@ -4,6 +4,8 @@ #include "python/dft_bindings.h" #include "python/qm_bindings.h" #include +#include +#include #include #include #include @@ -22,7 +24,22 @@ NB_MODULE(_occpy, m) { auto qm = register_qm_bindings(m); auto dft = register_dft_bindings(m); - m.def("setup_logging", [](int v) { occ::log::setup_logging(v); }); + nb::enum_(m, "LogLevel") + .value("TRACE", spdlog::level::level_enum::trace) + .value("DEBUG", spdlog::level::level_enum::debug) + .value("INFO", spdlog::level::level_enum::info) + .value("WARN", spdlog::level::level_enum::warn) + .value("ERROR", spdlog::level::level_enum::err) + .value("CRITICAL", spdlog::level::level_enum::critical) + .value("OFF", spdlog::level::level_enum::off); + + m.def("set_log_level", nb::overload_cast(occ::log::set_log_level)); + m.def("set_log_level", + nb::overload_cast(occ::log::set_log_level)); + m.def("set_log_level", + nb::overload_cast(occ::log::set_log_level)); + + m.def("set_log_file", occ::log::set_log_file); m.def("set_num_threads", [](int n) { occ::parallel::set_num_threads(n); }); m.def("set_data_directory", [](const std::string &s) { occ::qm::override_basis_set_directory(s); }); @@ -36,6 +53,6 @@ NB_MODULE(_occpy, m) { #ifdef VERSION_INFO m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); #else - m.attr("__version__") = "0.6.12"; + m.attr("__version__") = "0.6.13"; #endif } diff --git a/src/occpy/__init__.py b/src/occpy/__init__.py index 2d6953af3..fa802f2df 100644 --- a/src/occpy/__init__.py +++ b/src/occpy/__init__.py @@ -1,14 +1,20 @@ from pathlib import Path import site import warnings -from ._occpy import setup_logging, set_data_directory, set_num_threads +from ._occpy import ( + set_log_level, + set_log_file, + set_data_directory, + set_num_threads, + LogLevel, +) from .core import Atom, Element, Molecule, Dimer from .crystal import Crystal from .qm import AOBasis, HartreeFock, Shell, Wavefunction from .dft import DFT -# Set up logging first -setup_logging(0) +# Set up logging first, only log critical errors +set_log_level(LogLevel.CRITICAL) # Get site-packages directory and construct path to data _site_packages = Path(site.getsitepackages()[0]) @@ -32,11 +38,13 @@ "Dimer", "Element", "HartreeFock", + "LogLevel", "Molecule", "qm", "set_data_directory", "set_num_threads", - "setup_logging", + "set_log_file", + "set_log_level", "Shell", "Wavefunction", ] diff --git a/tests/interaction_tests.cpp b/tests/interaction_tests.cpp index 0f9bddefc..4f9e6d087 100644 --- a/tests/interaction_tests.cpp +++ b/tests/interaction_tests.cpp @@ -33,10 +33,10 @@ TEST_CASE("NaCl wolf sum", "[interaction,wolf]") { charges << 1.0, -1.0; double radius = 16.0; - double eta = 0.2; + double alpha = 0.2; auto surrounds = nacl.asymmetric_unit_atom_surroundings(radius); - occ::interaction::WolfParams params{radius, eta}; + occ::interaction::WolfParameters params{radius, alpha}; double wolf_energy = 0.0; occ::Mat3N cart_pos_asym = nacl.to_cartesian(nacl.asymmetric_unit().positions); @@ -47,7 +47,8 @@ TEST_CASE("NaCl wolf sum", "[interaction,wolf]") { for (int j = 0; j < qj.rows(); j++) { qj(j) = charges(surrounds[i].asym_idx(j)); } - wolf_energy += wolf_coulomb_energy(qi, pi, qj, surrounds[i].cart_pos); + wolf_energy += + wolf_coulomb_energy(qi, pi, qj, surrounds[i].cart_pos, params); } fmt::print("Wolf energy (NaCl): {}\n", wolf_energy); REQUIRE(wolf_energy == Catch::Approx(-0.3302735899347252));