diff --git a/3rdparty/CMakeLists.txt b/3rdparty/CMakeLists.txt index 9f2e1505b..a9bf207d7 100644 --- a/3rdparty/CMakeLists.txt +++ b/3rdparty/CMakeLists.txt @@ -1,27 +1,27 @@ CPMAddPackage( NAME fmt GITHUB_REPOSITORY "fmtlib/fmt" - VERSION "11.0.1" - GIT_TAG "11.0.1" + VERSION "11.0.2" + GIT_TAG "11.0.2" ) CPMAddPackage( NAME spdlog GITHUB_REPOSITORY "gabime/spdlog" OPTIONS "SPDLOG_FMT_EXTERNAL ON" - VERSION "1.14.1" + VERSION "1.15.0" ) CPMAddPackage( NAME unordered_dense GITHUB_REPOSITORY "martinus/unordered_dense" - VERSION "4.1.2" + VERSION "4.5.0" ) CPMAddPackage( NAME scnlib GITHUB_REPOSITORY "eliaskosunen/scnlib" - VERSION "2.0.2" + VERSION "4.0.1" OPTIONS "SIMDUTF_TESTS OFF" "ENABLE_FULL OFF" ) @@ -29,7 +29,7 @@ CPMAddPackage( CPMAddPackage( NAME CLI11 GITHUB_REPOSITORY "CLIUtils/CLI11" - VERSION "2.4.1" + VERSION "2.4.2" ) if("${USE_MLX}") diff --git a/cmake/CPM.cmake b/cmake/CPM.cmake index baf2d8c34..47e2e872a 100644 --- a/cmake/CPM.cmake +++ b/cmake/CPM.cmake @@ -2,8 +2,8 @@ # # SPDX-FileCopyrightText: Copyright (c) 2019-2023 Lars Melchior and contributors -set(CPM_DOWNLOAD_VERSION 0.40.2) -set(CPM_HASH_SUM "c8cdc32c03816538ce22781ed72964dc864b2a34a310d3b7104812a5ca2d835d") +set(CPM_DOWNLOAD_VERSION 0.40.5) +set(CPM_HASH_SUM "c46b876ae3b9f994b4f05a4c15553e0485636862064f1fcc9d8b4f832086bc5d") if(CPM_SOURCE_CACHE) set(CPM_DOWNLOAD_LOCATION "${CPM_SOURCE_CACHE}/cpm/CPM_${CPM_DOWNLOAD_VERSION}.cmake") diff --git a/include/occ/cg/crystal_growth_energies.h b/include/occ/cg/crystal_growth_energies.h deleted file mode 100644 index 412a143d8..000000000 --- a/include/occ/cg/crystal_growth_energies.h +++ /dev/null @@ -1,39 +0,0 @@ -#pragma once -#include -#include -#include -#include - -namespace occ::cg { - -using PairEnergies = std::vector; -using EnergyComponents = ankerl::unordered_dense::map; -using Energies = std::vector; - -struct EnergyTotal { - double crystal_energy{0.0}; - double interaction_energy{0.0}; - double solution_term{0.0}; -}; - -struct DimerSolventTerm { - double ab{0.0}; - double ba{0.0}; - double total{0.0}; -}; - -struct CrystalGrowthDimer { - core::Dimer dimer; - int unique_dimer_index{-1}; - double interaction_energy{0.0}; - DimerSolventTerm solvent_term{}; - double crystal_contribution{0.0}; - bool nearest_neighbor{false}; -}; - -struct CrystalGrowthResult { - std::vector> pair_energies; - std::vector total_energies; -}; - -} // namespace occ::cg diff --git a/include/occ/cg/interaction_mapper.h b/include/occ/cg/interaction_mapper.h index c0e895509..64aa7b72b 100644 --- a/include/occ/cg/interaction_mapper.h +++ b/include/occ/cg/interaction_mapper.h @@ -1,5 +1,5 @@ #pragma once -#include +#include #include #include @@ -20,13 +20,13 @@ class InteractionMapper { std::vector map_interactions(const std::vector &solution_terms, - const std::vector &interaction_energies_vec); + const std::vector &interaction_energies_vec); private: void map_molecule_interactions( size_t mol_idx, const Molecule &mol, const std::vector &solution_terms, - const std::vector &interaction_energies_vec, + const std::vector &interaction_energies_vec, std::vector &solution_terms_uc); void map_neighbor_interactions( @@ -34,12 +34,12 @@ class InteractionMapper { std::vector &unit_cell_neighbors, const std::vector &asymmetric_neighbors, - const Energies &interaction_energies); + const DimerResults &interaction_energies); void map_single_dimer(size_t mol_idx, size_t neighbor_idx, Dimer &dimer, const std::vector &asymmetric_neighbors, - const Energies &interaction_energies); + const DimerResults &interaction_energies); ankerl::unordered_dense::set build_related_set(const std::vector &related) const; @@ -57,7 +57,7 @@ class InteractionMapper { void update_dimer_properties( Dimer &dimer, size_t interaction_id, const std::vector &asym_dimers, - const Energies &energies) const; + const DimerResults &energies) const; void log_neighbor_info(size_t mol_idx, size_t uc_neighbors_size, size_t asym_neighbors_size, diff --git a/include/occ/cg/result_types.h b/include/occ/cg/result_types.h index 8496f2371..29fb29e8a 100644 --- a/include/occ/cg/result_types.h +++ b/include/occ/cg/result_types.h @@ -1,16 +1,21 @@ #pragma once #include +#include +#include #include #include namespace occ::cg { + +using PairEnergies = std::vector; using CGEnergyComponents = ankerl::unordered_dense::map; +using CGEnergies = std::vector; namespace components { constexpr const char *total = "Total"; constexpr const char *crystal_total = "Crystal Total"; -constexpr const char *crystal_nn = "Crystal (reassigned)"; +constexpr const char *crystal_nn = "Crystal (redistributed)"; constexpr const char *solvation_ab = "Solvation (A->B)"; constexpr const char *solvation_ba = "Solvation (B->A)"; constexpr const char *solvation_total = "Solvation Total"; @@ -22,9 +27,22 @@ constexpr const char *exchange = "Exchange"; } // namespace components +struct DimerSolventTerm { + double ab{0.0}; + double ba{0.0}; + double total{0.0}; +}; + +struct EnergyTotal { + double crystal_energy{0.0}; + double interaction_energy{0.0}; + double solution_term{0.0}; +}; + struct DimerResult { + occ::core::Dimer dimer; bool is_nearest_neighbor{false}; - size_t unique_idx{0}; + int unique_idx{0}; CGEnergyComponents energy_components{{components::total, 0.0}, {components::crystal_total, 0.0}, @@ -38,6 +56,8 @@ struct DimerResult { bool has_energy_component(const std::string &key) const; }; +using DimerResults = std::vector; + struct MoleculeResult { CGEnergyComponents energy_components{{components::total, 0.0}, {components::crystal_total, 0.0}, @@ -46,6 +66,8 @@ struct MoleculeResult { std::vector dimer_results; bool has_inversion_symmetry{true}; + cg::EnergyTotal total; + double total_energy() const; void add_dimer_result(const DimerResult &dimer); diff --git a/include/occ/core/dimer.h b/include/occ/core/dimer.h index d722f9a27..e97d2d183 100644 --- a/include/occ/core/dimer.h +++ b/include/occ/core/dimer.h @@ -27,6 +27,9 @@ class Dimer { BA = false /**< The B molecule data is first, then A */ }; + // Default (empty) constructor. + Dimer() = default; + /** * Constructor from two Molecule objects * @@ -200,7 +203,7 @@ class Dimer { * \param e a double representing the interaction energy. * \param key a string labelling the component of the interaction energy. */ - void set_interaction_energy(double e, const std::string &key = "total"); + void set_interaction_energy(double e, const std::string &key = "Total"); /** * Get the stored interaction energy for this Dimer. @@ -209,7 +212,7 @@ class Dimer { * * \returns a double representing the interaction energy (default = 0.0) */ - double interaction_energy(const std::string &key = "total") const; + double interaction_energy(const std::string &key = "Total") const; inline void set_interaction_energies( const ankerl::unordered_dense::map &e) { diff --git a/include/occ/driver/crystal_growth.h b/include/occ/driver/crystal_growth.h index 48f24c571..824d00325 100644 --- a/include/occ/driver/crystal_growth.h +++ b/include/occ/driver/crystal_growth.h @@ -1,6 +1,6 @@ #pragma once -#include #include +#include #include #include #include @@ -15,8 +15,6 @@ #include #include #include -#include -#include namespace occ::driver { @@ -109,8 +107,7 @@ class CrystalGrowthCalculator { virtual void init_monomer_energies() = 0; virtual void converge_lattice_energy() = 0; - virtual std::tuple> - process_neighbors_for_symmetry_unique_molecule( + virtual cg::MoleculeResult process_neighbors_for_symmetry_unique_molecule( int i, const std::string &molname) = 0; virtual cg::CrystalGrowthResult evaluate_molecular_surroundings() = 0; @@ -127,8 +124,8 @@ class CrystalGrowthCalculator { cg::PairEnergies m_dimer_energies; crystal::CrystalDimers m_nearest_dimers; std::vector m_solvation_breakdowns; - std::vector m_interaction_energies; - std::vector m_crystal_interaction_energies; + std::vector m_interaction_energies; + std::vector m_crystal_interaction_energies; std::vector m_solution_terms; private: @@ -146,8 +143,7 @@ class CEModelCrystalGrowthCalculator : public CrystalGrowthCalculator { cg::CrystalGrowthResult evaluate_molecular_surroundings() override; - std::tuple> - process_neighbors_for_symmetry_unique_molecule( + cg::MoleculeResult process_neighbors_for_symmetry_unique_molecule( int i, const std::string &molname) override; }; @@ -162,8 +158,7 @@ class XTBCrystalGrowthCalculator : public CrystalGrowthCalculator { cg::CrystalGrowthResult evaluate_molecular_surroundings() override; - std::tuple> - process_neighbors_for_symmetry_unique_molecule( + cg::MoleculeResult process_neighbors_for_symmetry_unique_molecule( int i, const std::string &molname) override; std::vector m_gas_phase_energies; diff --git a/include/occ/io/crystalgrower.h b/include/occ/io/crystalgrower.h index 53d8b097c..d5982e88c 100644 --- a/include/occ/io/crystalgrower.h +++ b/include/occ/io/crystalgrower.h @@ -1,35 +1,42 @@ #pragma once +#include #include #include -#include namespace occ::io::crystalgrower { class StructureWriter { - public: - StructureWriter(const std::string &filename); - StructureWriter(std::ostream &); - - void write(const occ::crystal::Crystal &, - const occ::crystal::CrystalDimers &); - void write_net(const occ::crystal::Crystal &, - const occ::crystal::CrystalDimers &); - - private: - std::ofstream m_owned_destination; - std::ostream &m_dest; +public: + StructureWriter(const std::string &filename); + StructureWriter(std::ostream &); + + void write(const occ::crystal::Crystal &, + const occ::crystal::CrystalDimers &); + void write_net(const occ::crystal::Crystal &, + const occ::crystal::CrystalDimers &); + +private: + std::ofstream m_owned_destination; + std::ostream &m_dest; }; class NetWriter { - public: - NetWriter(const std::string &filename); - NetWriter(std::ostream &); +public: + using InteractionLabels = + ankerl::unordered_dense::map; + NetWriter(const std::string &filename); + NetWriter(std::ostream &); + + void write(const occ::crystal::Crystal &, + const occ::crystal::CrystalDimers &); - void write(const occ::crystal::Crystal &, - const occ::crystal::CrystalDimers &); + inline const InteractionLabels &interaction_labels() const { + return m_interaction_labels; + }; - private: - std::ofstream m_owned_destination; - std::ostream &m_dest; +private: + std::ofstream m_owned_destination; + std::ostream &m_dest; + InteractionLabels m_interaction_labels; }; } // namespace occ::io::crystalgrower diff --git a/include/occ/main/occ_cg.h b/include/occ/main/occ_cg.h index f8428c187..5ebb1edff 100644 --- a/include/occ/main/occ_cg.h +++ b/include/occ/main/occ_cg.h @@ -1,6 +1,6 @@ #pragma once #include -#include +#include #include #include diff --git a/src/cg/interaction_mapper.cpp b/src/cg/interaction_mapper.cpp index 062b2f317..17ec00694 100644 --- a/src/cg/interaction_mapper.cpp +++ b/src/cg/interaction_mapper.cpp @@ -29,7 +29,7 @@ InteractionMapper::InteractionMapper(const Crystal &crystal, std::vector InteractionMapper::map_interactions( const std::vector &solution_terms, - const std::vector &interaction_energies_vec) { + const std::vector &interaction_energies_vec) { const auto &uc_molecules = m_crystal.unit_cell_molecules(); std::vector solution_terms_uc(m_uc_dimers.molecule_neighbors.size()); @@ -44,7 +44,7 @@ std::vector InteractionMapper::map_interactions( void InteractionMapper::map_molecule_interactions( size_t mol_idx, const Molecule &mol, const std::vector &solution_terms, - const std::vector &interaction_energies_vec, + const std::vector &interaction_energies_vec, std::vector &solution_terms_uc) { size_t asym_idx = mol.asymmetric_molecule_idx(); solution_terms_uc[mol_idx] = solution_terms[asym_idx]; @@ -65,7 +65,7 @@ void InteractionMapper::map_neighbor_interactions( std::vector &unit_cell_neighbors, const std::vector &asymmetric_neighbors, - const Energies &interaction_energies) { + const DimerResults &interaction_energies) { for (size_t j = 0; j < unit_cell_neighbors.size(); j++) { auto &[dimer, unique_idx] = unit_cell_neighbors[j]; map_single_dimer(mol_idx, j, dimer, asymmetric_neighbors, @@ -77,7 +77,7 @@ void InteractionMapper::map_single_dimer( size_t mol_idx, size_t neighbor_idx, Dimer &dimer, const std::vector &asymmetric_neighbors, - const Energies &interaction_energies) { + const DimerResults &interaction_energies) { const auto dimer_index = m_mapping_table.canonical_dimer_index(m_mapping_table.dimer_index(dimer)); const auto &related = m_mapping_table.symmetry_related_dimers(dimer_index); @@ -138,10 +138,12 @@ void InteractionMapper::handle_unmatched_dimer( void InteractionMapper::update_dimer_properties( Dimer &dimer, size_t interaction_id, const std::vector &asym_dimers, - const Energies &energies) const { + const DimerResults &energies) const { - dimer.set_property("asymmetric_dimer", fmt::format("dimer_{}", asym_dimers[interaction_id].unique_index)); - dimer.set_interaction_energies(energies[interaction_id]); + dimer.set_property( + "asymmetric_dimer", + fmt::format("dimer_{}", asym_dimers[interaction_id].unique_index)); + dimer.set_interaction_energies(energies[interaction_id].energy_components); dimer.set_interaction_id(interaction_id); } @@ -159,7 +161,7 @@ void InteractionMapper::log_neighbor_info(size_t mol_idx, void InteractionMapper::log_dimer_info(size_t neighbor_idx, const Dimer &dimer) const { - const auto& uc_mols = m_crystal.unit_cell_molecules(); + const auto &uc_mols = m_crystal.unit_cell_molecules(); auto idx_b = dimer.b().unit_cell_molecule_idx(); auto shift_b = dimer.b().cell_shift() - uc_mols[idx_b].cell_shift(); double rc = dimer.centroid_distance(); diff --git a/src/driver/crystal_growth.cpp b/src/driver/crystal_growth.cpp index abb231384..d2e05f95b 100644 --- a/src/driver/crystal_growth.cpp +++ b/src/driver/crystal_growth.cpp @@ -371,7 +371,7 @@ void CEModelCrystalGrowthCalculator::converge_lattice_energy() { } } -std::tuple> +cg::MoleculeResult CEModelCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( int i, const std::string &molname) { @@ -411,8 +411,8 @@ CEModelCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( occ::log::warn(std::string(95, '=')); size_t j = 0; - cg::EnergyTotal total; - std::vector dimer_energy_results; + cg::MoleculeResult dimer_energy_results; + auto &total = dimer_energy_results.total; size_t num_neighbors = std::accumulate( crystal_contributions.begin(), crystal_contributions.end(), 0, @@ -428,8 +428,8 @@ CEModelCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( for (const auto &[dimer, unique_idx] : full_neighbors) { const auto &e = m_dimer_energies[unique_idx]; if (!e.is_computed) { - interactions.push_back(cg::EnergyComponents{{"total", 0.0}}); - interactions_crystal.push_back(cg::EnergyComponents{{"total", 0.0}}); + interactions.push_back(cg::DimerResult{dimer, false, unique_idx}); + interactions_crystal.push_back(cg::DimerResult{dimer, false, unique_idx}); j++; continue; } @@ -467,33 +467,41 @@ CEModelCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( if (is_nearest_neighbor) { total.interaction_energy += interaction_energy; - interactions.push_back(cg::EnergyComponents{ - {"nn", e_nn}, - {"coulomb", e_coul}, - {"polarization", e_pol}, - {"repulsion", e_rep}, - {"exchange", e_exch}, - {"dispersion", e_disp}, - {"crystal", e_crys}, - {"solvent_ab", solvent_term.ab}, - {"solvent_ba", solvent_term.ba}, - {"solvent", solvent_term.total}, - {"total", interaction_energy}, - }); - - interactions_crystal.push_back(cg::EnergyComponents{ - {"nn", e_nn}, - {"coulomb", e_coul}, - {"polarization", e_pol}, - {"repulsion", e_rep}, - {"exchange", e_exch}, - {"dispersion", e_disp}, - {"crystal", e_crys}, - {"total", e_crys + e_nn}, - }); + interactions.push_back(cg::DimerResult{ + dimer, + true, + unique_idx, + { + {cg::components::crystal_nn, e_nn}, + {cg::components::coulomb, e_coul}, + {cg::components::polarization, e_pol}, + {cg::components::repulsion, e_rep}, + {cg::components::exchange, e_exch}, + {cg::components::dispersion, e_disp}, + {cg::components::crystal_total, e_crys}, + {cg::components::solvation_ab, solvent_term.ab}, + {cg::components::solvation_ba, solvent_term.ba}, + {cg::components::solvation_total, solvent_term.total}, + {cg::components::total, interaction_energy}, + }}); + + interactions_crystal.push_back( + cg::DimerResult{dimer, + true, + unique_idx, + { + {cg::components::crystal_nn, e_nn}, + {cg::components::coulomb, e_coul}, + {cg::components::polarization, e_pol}, + {cg::components::repulsion, e_rep}, + {cg::components::exchange, e_exch}, + {cg::components::dispersion, e_disp}, + {cg::components::crystal_total, e_crys}, + {cg::components::total, e_crys + e_nn}, + }}); } else { - interactions.push_back(cg::EnergyComponents{{"total", 0.0}}); - interactions_crystal.push_back(cg::EnergyComponents{{"total", 0.0}}); + interactions.push_back(cg::DimerResult{dimer, false, unique_idx}); + interactions_crystal.push_back(cg::DimerResult{dimer, false, unique_idx}); interaction_energy = 0; } @@ -508,13 +516,11 @@ CEModelCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( solvent_term.ba, solvent_term.total, crystal_contribution, interaction_energy); } - dimer_energy_results.emplace_back(cg::CrystalGrowthDimer{ - dimer, unique_idx, interaction_energy, solvent_term, - crystal_contribution, is_nearest_neighbor}); + dimer_energy_results.add_dimer_result(interactions.back()); j++; } m_solvation_breakdowns.push_back(solvation_breakdown); - return {total, dimer_energy_results}; + return dimer_energy_results; } cg::CrystalGrowthResult @@ -524,17 +530,16 @@ CEModelCrystalGrowthCalculator::evaluate_molecular_surroundings() { m_solution_terms = std::vector(m_molecules.size(), 0.0); for (size_t i = 0; i < m_molecules.size(); i++) { - auto [mol_total, mol_dimer_results] = - process_neighbors_for_symmetry_unique_molecule( - i, fmt::format("{}_{}_{}", opts.basename, i, opts.solvent)); + auto mol_dimer_results = process_neighbors_for_symmetry_unique_molecule( + i, fmt::format("{}_{}_{}", opts.basename, i, opts.solvent)); - result.pair_energies.push_back(mol_dimer_results); - result.total_energies.push_back(mol_total); + result.molecule_results.push_back(mol_dimer_results); - m_solution_terms[i] = mol_total.solution_term; - m_lattice_energies.push_back(mol_total.crystal_energy); - write_energy_summary(mol_total.crystal_energy, m_molecules[i], - mol_total.solution_term, mol_total.interaction_energy); + m_solution_terms[i] = mol_dimer_results.total.solution_term; + m_lattice_energies.push_back(mol_dimer_results.total.crystal_energy); + write_energy_summary(mol_dimer_results.total.crystal_energy, m_molecules[i], + mol_dimer_results.total.solution_term, + mol_dimer_results.total.interaction_energy); if (opts.write_debug_output_files) { // write neighbors file for molecule i @@ -611,18 +616,17 @@ XTBCrystalGrowthCalculator::evaluate_molecular_surroundings() { m_solution_terms = std::vector(m_molecules.size(), 0.0); for (size_t i = 0; i < m_molecules.size(); i++) { - auto [mol_total, mol_dimer_results] = - process_neighbors_for_symmetry_unique_molecule( - i, fmt::format("{}_{}_{}", opts.basename, i, opts.solvent)); - - result.pair_energies.push_back(mol_dimer_results); - result.total_energies.push_back(mol_total); - - m_solution_terms[i] = mol_total.solution_term; - m_lattice_energies.push_back(mol_total.crystal_energy); - occ::driver::write_energy_summary(mol_total.crystal_energy, m_molecules[i], - mol_total.solution_term, - mol_total.interaction_energy); + auto mol_dimer_results = process_neighbors_for_symmetry_unique_molecule( + i, fmt::format("{}_{}_{}", opts.basename, i, opts.solvent)); + + result.molecule_results.push_back(mol_dimer_results); + + m_solution_terms[i] = mol_dimer_results.total.solution_term; + m_lattice_energies.push_back(mol_dimer_results.total.crystal_energy); + occ::driver::write_energy_summary( + mol_dimer_results.total.crystal_energy, m_molecules[i], + mol_dimer_results.total.solution_term, + mol_dimer_results.total.interaction_energy); } return result; } @@ -673,7 +677,7 @@ void XTBCrystalGrowthCalculator::init_monomer_energies() { sw_solv.read()); } -std::tuple> +cg::MoleculeResult XTBCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( int i, const std::string &molname) { @@ -695,8 +699,9 @@ XTBCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( return x.is_nn ? a + 1 : a; }); - occ::cg::EnergyTotal total; - std::vector dimer_energy_results; + cg::MoleculeResult dimer_energy_results; + auto &total = dimer_energy_results.total; + total.solution_term = (m_solvated_energies[i] - m_gas_phase_energies[i]) * occ::units::AU_TO_KJ_PER_MOL; @@ -745,18 +750,26 @@ XTBCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( if (is_nearest_neighbor) { total.interaction_energy += interaction_energy; - interactions.push_back(cg::EnergyComponents{ - {"nn", crystal_contributions[j].energy}, - {"crystal", e}, - {"solv", solvent_term.total}, - {"total", interaction_energy}, - }); - - interactions_crystal.push_back(cg::EnergyComponents{ - {"crystal", e}, {"total", e + crystal_contributions[j].energy}}); + interactions.push_back(cg::DimerResult{ + dimer, + true, + unique_idx, + { + {cg::components::crystal_nn, crystal_contributions[j].energy}, + {cg::components::crystal_total, e}, + {cg::components::solvation_total, solvent_term.total}, + {cg::components::total, interaction_energy}, + }}); + + interactions_crystal.push_back(cg::DimerResult{ + dimer, + true, + unique_idx, + {{cg::components::crystal_total, e}, + {cg::components::total, e + crystal_contributions[j].energy}}}); } else { - interactions.push_back(cg::EnergyComponents{{"total", 0.0}}); - interactions_crystal.push_back(cg::EnergyComponents{{"total", 0.0}}); + interactions.push_back(cg::DimerResult{dimer, false, unique_idx}); + interactions_crystal.push_back(cg::DimerResult{dimer, false, unique_idx}); interaction_energy = 0; } @@ -774,12 +787,10 @@ XTBCrystalGrowthCalculator::process_neighbors_for_symmetry_unique_molecule( interaction_energy); } - dimer_energy_results.emplace_back(occ::cg::CrystalGrowthDimer{ - dimer, unique_idx, interaction_energy, solvent_term, - crystal_contribution, is_nearest_neighbor}); + dimer_energy_results.add_dimer_result(interactions.back()); j++; } - return {total, dimer_energy_results}; + return dimer_energy_results; } } // namespace occ::driver diff --git a/src/interaction/pair_energy.cpp b/src/interaction/pair_energy.cpp index f624dff9a..669e43a80 100644 --- a/src/interaction/pair_energy.cpp +++ b/src/interaction/pair_energy.cpp @@ -661,12 +661,12 @@ converged_lattice_energies(EnergyModel &energy_model, const Crystal &crystal, } 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"); + 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}; } @@ -710,4 +710,4 @@ bool PairEnergyStore::load(int id, const Dimer &d, CEEnergyComponents &e) { e); } -} // namespace occ::main +} // namespace occ::interaction diff --git a/src/io/crystalgrower.cpp b/src/io/crystalgrower.cpp index 8bbaf6a20..2b810facc 100644 --- a/src/io/crystalgrower.cpp +++ b/src/io/crystalgrower.cpp @@ -94,40 +94,6 @@ NetWriter::NetWriter(const std::string &filename) NetWriter::NetWriter(std::ostream &stream) : m_dest(stream) {} -struct FormulaIndex { - int count{1}; - char letter{'A'}; - bool operator==(const FormulaIndex &o) const { - return count == o.count && letter == o.letter; - } -}; - -struct FormulaIndexHash { - uint64_t operator()(const FormulaIndex &x) const noexcept { - return ankerl::unordered_dense::detail::wyhash::mix(x.count, x.letter); - } -}; - -std::vector -build_formula_indices_for_symmetry_unique_molecules(const Crystal &crystal) { - std::vector result; - ankerl::unordered_dense::map formula_count; - for (const auto &mol : crystal.symmetry_unique_molecules()) { - std::string formula = occ::core::chemical_formula(mol.elements()); - auto it = formula_count.find(formula); - FormulaIndex formula_index; - if (it != formula_count.end()) { - it->second.count++; - formula_index = it->second; - } else { - formula_index.letter += static_cast(formula_count.size()); - formula_count.insert({formula, formula_index}); - } - result.push_back(formula_index); - } - return result; -} - inline std::string get_asymmetric_dimer_name(const core::Dimer &n) { const auto &properties = n.properties(); std::string name = ""; @@ -139,6 +105,7 @@ inline std::string get_asymmetric_dimer_name(const core::Dimer &n) { } void NetWriter::write(const Crystal &crystal, const CrystalDimers &uc_dimers) { + m_interaction_labels.clear(); const auto &uc_molecules = crystal.unit_cell_molecules(); CrystalDimers uc_dimers_copy = uc_dimers; sort_neighbors(uc_dimers_copy); @@ -151,27 +118,12 @@ void NetWriter::write(const Crystal &crystal, const CrystalDimers &uc_dimers) { size_t uc_idx = 0; constexpr double max_de = 1e-4; - std::vector sym_formula_indices = - build_formula_indices_for_symmetry_unique_molecules(crystal); + ankerl::unordered_dense::map> f2e; - ankerl::unordered_dense::map, - FormulaIndexHash> - f2e; - - // TODO fix for multiple molecules in asymmetric unit - // 1A -> 1 is the conformer number, A is the compound i.e. chemical - // composition id for (const auto &mol : uc_molecules) { - FormulaIndex formula_index = - sym_formula_indices[mol.asymmetric_molecule_idx()]; - if (!f2e.contains(formula_index)) { - f2e.insert({formula_index, {}}); - } + std::string formula_index = mol.name(); std::vector &unique_interaction_energies = f2e[formula_index]; - occ::log::debug("uc mol {}, asym = {}, formula = {} {}", uc_idx, - mol.asymmetric_molecule_idx(), formula_index.count, - formula_index.letter); auto mol_neighbors = neighbors[uc_idx]; for (const auto &[n, unique_index] : neighbors[uc_idx]) { const auto uc_shift = n.b().cell_shift(); @@ -192,11 +144,14 @@ void NetWriter::write(const Crystal &crystal, const CrystalDimers &uc_dimers) { 1 + std::distance(unique_interaction_energies.begin(), match); } - auto dimer_name = get_asymmetric_dimer_name(n); + auto asym_dimer_filename = get_asymmetric_dimer_name(n); + auto dimer_label = dimer_labeller(n); + m_interaction_labels.insert( + {dimer_label, fmt::format("{}-{}", formula_index, interaction_idx)}); - fmt::print(m_dest, "{}:[{}{}][{}] {} R={:.3f}\n", interaction_idx, - formula_index.count, formula_index.letter, dimer_labeller(n), - dimer_name, n.centroid_distance()); + fmt::print(m_dest, "{}:[{}][{}] {} R={:.3f}\n", interaction_idx, + formula_index, dimer_label, asym_dimer_filename, + n.centroid_distance()); } for (double e : unique_interaction_energies) { fmt::print(m_dest, "{:7.3f}\n", e * occ::units::KJ_TO_KCAL); diff --git a/src/main/occ_cg.cpp b/src/main/occ_cg.cpp index 5b0852874..36937f833 100644 --- a/src/main/occ_cg.cpp +++ b/src/main/occ_cg.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -41,36 +42,38 @@ namespace fs = std::filesystem; using occ::crystal::Crystal; using occ::crystal::CrystalDimers; - using occ::driver::WavefunctionChoice; -void write_cg_structure_file(const std::string &filename, - const Crystal &crystal, - const CrystalDimers &uc_dimers) { +inline void write_cg_structure_file(const std::string &filename, + const Crystal &crystal, + const CrystalDimers &uc_dimers) { occ::log::info("Writing crystalgrower structure file to '{}'", filename); occ::io::crystalgrower::StructureWriter cg_structure_writer(filename); cg_structure_writer.write(crystal, uc_dimers); } -void write_cg_net_file(const std::string &filename, const Crystal &crystal, - const CrystalDimers &uc_dimers) { +inline auto write_cg_net_file(const std::string &filename, + const Crystal &crystal, + const CrystalDimers &uc_dimers) { occ::log::info("Writing crystalgrower net file to '{}'", filename); occ::io::crystalgrower::NetWriter cg_net_writer(filename); cg_net_writer.write(crystal, uc_dimers); + return cg_net_writer.interaction_labels(); } -void write_kmcpp_input_file(const std::string &filename, const Crystal &crystal, - const CrystalDimers &uc_dimers, - const std::vector &solution_terms) { +inline void write_kmcpp_input_file(const std::string &filename, + const Crystal &crystal, + const CrystalDimers &uc_dimers, + const std::vector &solution_terms) { occ::log::info("Writing kmcpp structure file to '{}'", filename); occ::io::kmcpp::InputWriter kmcpp_structure_writer(filename); kmcpp_structure_writer.write(crystal, uc_dimers, solution_terms); } -std::vector map_unique_interactions_to_uc_molecules( +inline std::vector map_unique_interactions_to_uc_molecules( const Crystal &crystal, const CrystalDimers &dimers, CrystalDimers &uc_dimers, const std::vector &solution_terms, - const std::vector &interaction_energies_vec, + const std::vector &interaction_energies_vec, bool inversion) { occ::cg::InteractionMapper mapper(crystal, dimers, uc_dimers, inversion); @@ -169,9 +172,12 @@ inline void write_wulff(const std::string &filename, occ::io::write_ply_mesh(filename, mesh, {}, false); } -void write_cg_dimers(const occ::driver::CrystalGrowthCalculatorOptions &opts, - const occ::crystal::Crystal &crystal, - const occ::cg::CrystalGrowthResult &result) { +inline void write_cg_dimers( + const occ::driver::CrystalGrowthCalculatorOptions &opts, + const occ::crystal::Crystal &crystal, + const occ::cg::CrystalGrowthResult &result, + const occ::io::crystalgrower::NetWriter::InteractionLabels &cg_labels) { + nlohmann::json j; j["result_type"] = "cg"; j["title"] = opts.basename; @@ -182,7 +188,8 @@ void write_cg_dimers(const occ::driver::CrystalGrowthCalculatorOptions &opts, j["crystal"] = crystal; j["totals_per_molecule"] = {}; - for (const auto &mol_total : result.total_energies) { + for (const auto &mol_result : result.molecule_results) { + const auto &mol_total = mol_result.total; nlohmann::json e; e["crystal_energy"] = mol_total.crystal_energy; e["interaction_energy"] = mol_total.interaction_energy; @@ -192,24 +199,34 @@ void write_cg_dimers(const occ::driver::CrystalGrowthCalculatorOptions &opts, const auto &uc_atoms = crystal.unit_cell_atoms(); + auto dimer_labeller = occ::crystal::SymmetryDimerLabeller(crystal); + dimer_labeller.connection = "-"; + dimer_labeller.format.fmt_string = "{}"; + j["pairs"] = {}; - for (const auto &mol_pairs : result.pair_energies) { + for (const auto &mol_result : result.molecule_results) { nlohmann::json m; - for (const auto &cg_dimer : mol_pairs) { + for (const auto &dimer_result : mol_result.dimer_results) { + const auto &dimer = dimer_result.dimer; nlohmann::json d; nlohmann::json e; - e["unique_dimer_index"] = cg_dimer.unique_dimer_index; - e["interaction_energy"] = cg_dimer.interaction_energy; - e["crystal_contribution"] = cg_dimer.crystal_contribution; - e["is_nearest_neighbor"] = cg_dimer.nearest_neighbor; - e["solvent_ab"] = cg_dimer.solvent_term.ab; - e["solvent_ba"] = cg_dimer.solvent_term.ba; - e["solvent_total"] = cg_dimer.solvent_term.total; + auto label = dimer_labeller(dimer); + std::string cg_id{"??"}; + const auto kv = cg_labels.find(label); + if (kv != cg_labels.end()) { + cg_id = kv->second; + } + for (const auto &[k, v] : dimer_result.energy_components) { + e[k] = v; + } + d["Nearest Neighbor"] = dimer_result.is_nearest_neighbor; + d["Unique Index"] = dimer_result.unique_idx; + d["Crystalgrower Identifier"] = cg_id; d["energies"] = e; nlohmann::json offsets_a = {}; { - const auto &a = cg_dimer.dimer.a(); + const auto &a = dimer.a(); const auto &a_uc_idx = a.unit_cell_idx(); const auto &a_uc_shift = a.unit_cell_shift(); for (int i = 0; i < a_uc_idx.rows(); i++) { @@ -220,7 +237,7 @@ void write_cg_dimers(const occ::driver::CrystalGrowthCalculatorOptions &opts, } nlohmann::json offsets_b = {}; { - const auto &b = cg_dimer.dimer.b(); + const auto &b = dimer.b(); const auto &b_uc_idx = b.unit_cell_idx(); const auto &b_uc_shift = b.unit_cell_shift(); for (int i = 0; i < b_uc_idx.rows(); i++) { @@ -307,6 +324,13 @@ occ::cg::CrystalGrowthResult run_cg_impl(CGConfig const &config) { calc.crystal(), calc.full_dimers(), uc_dimers, calc.solution_terms(), calc.interaction_energies(), !opts.use_asymmetric_partition); + // TODO tidy this up, but for now just do the same thing for crystal energies + // too so we get vacuum surface energies + auto vacuum_terms_uc = map_unique_interactions_to_uc_molecules( + calc.crystal(), calc.full_dimers(), uc_dimers_vacuum, + calc.solution_terms(), calc.crystal_interaction_energies(), + !opts.use_asymmetric_partition); + if (config.write_kmcpp_file) { write_kmcpp_input_file(fmt::format("{}_kmcpp.json", basename), calc.crystal(), uc_dimers, solution_terms_uc); @@ -329,8 +353,27 @@ occ::cg::CrystalGrowthResult run_cg_impl(CGConfig const &config) { surface_energies); write_wulff(fmt::format("{}_vacuum.ply", basename, config.solvent), vacuum_surface_energies); - j["vacuum"] = calc.crystal_interaction_energies(); - j["solvated"] = calc.interaction_energies(); + + // TODO refactor this + nlohmann::json vacuum_energies; + for (const auto &mol : calc.crystal_interaction_energies()) { + nlohmann::json tmp; + for (const auto &v : mol) { + tmp.push_back(v.energy_components); + } + vacuum_energies.push_back(tmp); + } + nlohmann::json solvated_energies; + for (const auto &mol : calc.interaction_energies()) { + nlohmann::json tmp; + for (const auto &v : mol) { + tmp.push_back(v.energy_components); + } + solvated_energies.push_back(tmp); + } + j["vacuum"] = vacuum_energies; + j["solvated"] = solvated_energies; + std::string surf_energy_filename = fmt::format("{}_surface_energies.json", basename); occ::log::info("Writing surface energies to '{}'", surf_energy_filename); @@ -338,9 +381,11 @@ occ::cg::CrystalGrowthResult run_cg_impl(CGConfig const &config) { destination << j.dump(2); } - write_cg_dimers(opts, calc.crystal(), result); - write_cg_net_file(fmt::format("{}_{}_net.txt", basename, config.solvent), - calc.crystal(), uc_dimers); + auto cg_interaction_labels = + write_cg_net_file(fmt::format("{}_{}_net.txt", basename, config.solvent), + calc.crystal(), uc_dimers); + + write_cg_dimers(opts, calc.crystal(), result, cg_interaction_labels); // calc.dipole_correction(); return result; diff --git a/src/occpy.cpp b/src/occpy.cpp index 3c8362d1a..077fbd854 100644 --- a/src/occpy.cpp +++ b/src/occpy.cpp @@ -395,9 +395,10 @@ NB_MODULE(_occpy, m) { uc.cell_type(), uc.a(), uc.b(), uc.c()); }); - using occ::cg::CrystalGrowthDimer; using occ::cg::CrystalGrowthResult; + using occ::cg::DimerResult; using occ::cg::DimerSolventTerm; + using occ::cg::MoleculeResult; using occ::interaction::LatticeConvergenceSettings; using occ::main::CGConfig; @@ -428,17 +429,22 @@ NB_MODULE(_occpy, m) { .def_ro("ba", &DimerSolventTerm::ba) .def_ro("total", &DimerSolventTerm::total); - nb::class_(m, "CrystalGrowthDimer") - .def_ro("dimer", &CrystalGrowthDimer::dimer) - .def_ro("unique_dimer_index", &CrystalGrowthDimer::unique_dimer_index) - .def_ro("interaction_energy", &CrystalGrowthDimer::interaction_energy) - .def_ro("solvent_term", &CrystalGrowthDimer::solvent_term) - .def_ro("crystal_contribution", &CrystalGrowthDimer::crystal_contribution) - .def_ro("nearest_neighbor", &CrystalGrowthDimer::nearest_neighbor); + nb::class_(m, "DimerResult") + .def_ro("dimer", &DimerResult::dimer) + .def_ro("unique_idx", &DimerResult::unique_idx) + .def("total_energy", &DimerResult::total_energy) + .def("energy_component", &DimerResult::energy_component) + .def_ro("is_nearest_neighbor", &DimerResult::is_nearest_neighbor); + + nb::class_(m, "MoleculeResult") + .def_ro("dimer_results", &MoleculeResult::dimer_results) + .def_ro("total", &MoleculeResult::total) + .def_ro("has_inversion_symmetry", &MoleculeResult::has_inversion_symmetry) + .def("total_energy", &MoleculeResult::total_energy) + .def("energy_component", &MoleculeResult::energy_component); nb::class_(m, "CrystalGrowthResult") - .def_ro("pair_energies", &CrystalGrowthResult::pair_energies) - .def_ro("total_energies", &CrystalGrowthResult::total_energies); + .def_ro("molecule_results", &CrystalGrowthResult::molecule_results); nb::class_(m, "CrystalGrowthEnergyTotal") .def_ro("crystal", &occ::cg::EnergyTotal::crystal_energy)