// pybind11 bindings for the vibe-qc C++ core. #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "vibeqc/build_config.hpp" #include "vibeqc/molecule.hpp" #include "vibeqc/basis.hpp" #include "vibeqc/ao_eval.hpp" #include "vibeqc/bloch.hpp" #include "vibeqc/crystal.hpp" #include "vibeqc/dft_plus_u.hpp" #include "vibeqc/dispersion.hpp" #include "vibeqc/eeq_charges.hpp" #include "vibeqc/semiempirical/dftb0.hpp" #include "vibeqc/semiempirical/gradient.hpp" #include "vibeqc/semiempirical/periodic_dftb0.hpp" #include "vibeqc/semiempirical/periodic_scc_dftb.hpp" #include "vibeqc/semiempirical/kpoints_dftb0.hpp" #include "vibeqc/semiempirical/kpoints_scc_dftb.hpp" #include "vibeqc/semiempirical/repulsive_spline.hpp" #include "vibeqc/semiempirical/core/method_registry.hpp" #include "vibeqc/semiempirical/core/parameters.hpp" #include "vibeqc/semiempirical/methods/xtb/gfn2_parameters.hpp" #include "vibeqc/semiempirical/methods/xtb/gfn2_driver.hpp" #include "vibeqc/semiempirical/methods/xtb/periodic_gfn2.hpp" #include "vibeqc/semiempirical/methods/xtb/ugfn2.hpp" #include "vibeqc/semiempirical/methods/nddo/pm6_parameters.hpp" #include "vibeqc/semiempirical/methods/nddo/pm6_fock.hpp" #include "vibeqc/semiempirical/methods/nddo/omx_parameters.hpp" #include "vibeqc/semiempirical/methods/nddo/omx_fock.hpp" #include "vibeqc/semiempirical/methods/nddo/periodic_pm6.hpp" #include "vibeqc/semiempirical/methods/nddo/periodic_omx.hpp" #include "vibeqc/semiempirical/core/hamiltonian_builders.hpp" #include "vibeqc/semiempirical/basis.hpp" #include "vibeqc/semiempirical/scc_dftb.hpp" #include "vibeqc/semiempirical/parameters.hpp" #include "vibeqc/ecp.hpp" #include "vibeqc/ewald.hpp" #include "vibeqc/fft_poisson.hpp" #include "vibeqc/lattice_integrals.hpp" #include "vibeqc/lattice_sum.hpp" #include "vibeqc/multipole_moments_lattice.hpp" #include "vibeqc/spheropole.hpp" #include "vibeqc/jk_builder.hpp" #include "vibeqc/periodic.hpp" #include "vibeqc/periodic_fock.hpp" #include "vibeqc/periodic_jk_builder.hpp" #include "vibeqc/periodic_rhf.hpp" #include "vibeqc/periodic_scf.hpp" #include "vibeqc/periodic_xc.hpp" #include "vibeqc/thread_pool.hpp" #include "vibeqc/fock.hpp" #include "vibeqc/gradient.hpp" #include "vibeqc/periodic_gradient.hpp" #include "vibeqc/hessian_integrals.hpp" #include "vibeqc/grid.hpp" #include "vibeqc/guess.hpp" #include "vibeqc/spheropole.hpp" #include "vibeqc/init.hpp" #include "vibeqc/integrals.hpp" #include "vibeqc/cosx.hpp" #include "vibeqc/cosx_robust.hpp" #include "vibeqc/cosx_kernel.hpp" #include "vibeqc/grid_batch.hpp" #include "vibeqc/schwarz.hpp" #include "vibeqc/df.hpp" #include "vibeqc/aux_eri.hpp" #include "vibeqc/diis.hpp" #include "vibeqc/kdiis.hpp" #include "vibeqc/mp2.hpp" #include "vibeqc/ccsd.hpp" #include "vibeqc/ediis.hpp" #include "vibeqc/newton.hpp" #include "vibeqc/rhf.hpp" #include "vibeqc/rks.hpp" #include "vibeqc/uhf.hpp" #include "vibeqc/uks.hpp" #include "vibeqc/ump2.hpp" #include "vibeqc/xc.hpp" #include "vibeqc/xc_kernel.hpp" namespace py = pybind11; PYBIND11_MODULE(_vibeqc_core, m) { m.doc() = "vibeqc native core (C++ backend)"; // ----- Smoke / diagnostics ------------------------------------------- m.def("hello", []() { return std::string{"vibeqc core alive"}; }, "Smoke test: returns a static greeting."); m.def("libint_version", []() { return libint2::libint_version_string(false); }, "libint2 version string, e.g. \"2.13.1\"."); m.def("libxc_version", []() { return std::string{xc_version_string()}; }, "libxc version string, e.g. \"7.0.0\"."); m.def("blas_info", []() { py::dict info; info["libraries"] = std::string{ vibeqc::build_config::kBlasLibraries}; info["blas_enabled"] = vibeqc::build_config::kBlasEnabled; info["lapacke_enabled"] = vibeqc::build_config::kLapackeEnabled; return info; }, "Build-time BLAS/LAPACK linkage info: dict with keys " "'libraries' (raw BLAS_LIBRARIES string from CMake's FindBLAS — " "e.g. '-framework Accelerate' on macOS, '/usr/lib/libblas.so' " "on Linux netlib, empty when no BLAS was linked), " "'blas_enabled' (True iff EIGEN_USE_BLAS was set), and " "'lapacke_enabled' (True iff EIGEN_USE_LAPACKE was set; " "controls whether Eigen's dense solvers delegate to LAPACKE)."); m.def("libecpint_version", &vibeqc::libecpint_version, "libecpint version string. Phase 14 ECP integrals depend on " "this library; vibe-qc vendors v1.0.7 under " "third_party/libecpint/ alongside its pugixml + libcerf " "dependencies for byte-identical builds across machines."); m.def("fftw3_version", &vibeqc::fftw3_version, "FFTW3 version string, e.g. \"3.3.10-sse2-avx\". The " "FFT-Poisson long-range Hartree solver " "(cpp/src/fft_poisson.cpp) and the upcoming GAPW route " "both depend on FFTW3; surfaced through the banner per " "CLAUDE.md § 6."); py::class_(m, "ECPCenter", "Per-atom ECP placement: atomic number Z + Cartesian " "position (bohr).") .def(py::init<>()) .def(py::init([](int Z, std::array xyz) { return vibeqc::ECPCenter{Z, xyz}; }), py::arg("Z"), py::arg("xyz")) .def_readwrite("Z", &vibeqc::ECPCenter::Z) .def_readwrite("xyz", &vibeqc::ECPCenter::xyz); m.def("compute_ecp_matrix", &vibeqc::compute_ecp_matrix, py::arg("basis"), py::arg("ecp_centers"), py::arg("library_name") = std::string("ecp10mdf"), py::arg("share_dir") = std::string(""), py::call_guard(), "Compute the AO-basis ECP matrix V_ECP_{μν} = ⟨χ_μ|V_ECP|χ_ν⟩ " "via libecpint's built-in XML library (ecp10mdf, ecp28mdf, " "ecp46mdf, ecp60mdf, ecp78mdf, lanl2dz). Output is in " "spherical (real-solid-harmonic) basis. share_dir defaults " "to the vendored third_party/libecpint/install/share/libecpint."); m.def("get_num_threads", &vibeqc::omp_max_threads, "Maximum number of OpenMP threads the next parallel region " "will use."); m.def("set_num_threads", &vibeqc::set_num_threads, py::arg("n"), "Set the OpenMP thread count; ``n <= 0`` restores the " "OMP_NUM_THREADS-derived default. Returns the actual thread " "count that will be used."); // ----- Atom ---------------------------------------------------------- py::class_(m, "Atom") .def(py::init<>()) .def(py::init([](int Z, std::array xyz) { return vibeqc::Atom{Z, xyz}; }), py::arg("Z"), py::arg("xyz"), "Create an Atom with atomic number Z and position in bohr.") .def_readwrite("Z", &vibeqc::Atom::Z) .def_readwrite("xyz", &vibeqc::Atom::xyz) .def("__repr__", [](const vibeqc::Atom& a) { return "Atom(Z=" + std::to_string(a.Z) + ", xyz=[" + std::to_string(a.xyz[0]) + ", " + std::to_string(a.xyz[1]) + ", " + std::to_string(a.xyz[2]) + "])"; }); // ----- Molecule ------------------------------------------------------ py::class_(m, "Molecule") .def(py::init, int, int>(), py::arg("atoms"), py::arg("charge") = 0, py::arg("multiplicity") = 1) .def_property_readonly("atoms", &vibeqc::Molecule::atoms) .def_property_readonly("charge", &vibeqc::Molecule::charge) .def_property_readonly("multiplicity", &vibeqc::Molecule::multiplicity) .def("n_electrons", &vibeqc::Molecule::n_electrons) .def("nuclear_repulsion", &vibeqc::Molecule::nuclear_repulsion, "Point-charge nuclear repulsion energy in Hartree."); // ----- BasisSet ------------------------------------------------------ py::class_(m, "ShellInfo", "One contracted Gaussian shell: angular momentum, primitive " "exponents, contraction coefficients, and the atom it sits on. " "Produced by BasisSet.shells(); also constructible from Python so " "callers can build modrho-renormalised aux bases or any other " "custom-coefficient shell list and feed it to " "BasisSet(molecule, shells, name, coefficients_pre_normalized).") .def(py::init<>()) .def(py::init([](int atom_index, int l, bool pure, std::vector exponents, std::vector coefficients, std::array origin) { vibeqc::ShellInfo info; info.atom_index = atom_index; info.l = l; info.pure = pure; info.exponents = std::move(exponents); info.coefficients = std::move(coefficients); info.origin = origin; return info; }), py::arg("atom_index"), py::arg("l"), py::arg("pure"), py::arg("exponents"), py::arg("coefficients"), py::arg("origin"), "Construct a ShellInfo from explicit per-shell metadata. " "Field semantics match the read-only properties; see " "BasisSet(molecule, shells, ...) for how this is consumed.") .def_readwrite("atom_index", &vibeqc::ShellInfo::atom_index, "Index of the owning atom in Molecule.atoms() (0-based).") .def_readwrite("l", &vibeqc::ShellInfo::l, "Angular momentum: 0=s, 1=p, 2=d, 3=f, ...") .def_readwrite("pure", &vibeqc::ShellInfo::pure, "True: spherical harmonics (2L+1 AOs). " "False: Cartesian ((L+1)(L+2)/2 AOs).") .def_readwrite("exponents", &vibeqc::ShellInfo::exponents, "Primitive Gaussian exponents α_i (bohr^-2).") .def_readwrite("coefficients", &vibeqc::ShellInfo::coefficients, "Contraction coefficients c_i (same length as " "exponents). When read from BasisSet.shells() these " "include libint's primitive normalisation; pass " "coefficients_pre_normalized=True (default) when " "constructing a derived BasisSet from a modified " "shell list to keep them as-is.") .def_readwrite("origin", &vibeqc::ShellInfo::origin, "Cartesian origin (bohr).") .def("__repr__", [](const vibeqc::ShellInfo& s) { const char lch[] = "spdfghikl"; char lc = (s.l >= 0 && s.l < 9) ? lch[s.l] : '?'; return std::string{"ShellInfo("} + lc + ", atom=" + std::to_string(s.atom_index) + ", n_prim=" + std::to_string(s.exponents.size()) + (s.pure ? ", pure" : ", cart") + ")"; }); py::class_(m, "BasisSet") .def(py::init(), py::arg("molecule"), py::arg("name")) // Build from an explicit ShellInfo list. coefficients_pre_normalized // defaults to True (round-trip from .shells() preserves them as-is); // set False when feeding shells parsed from a Gaussian-style file // format that needs libint to embed primitive normalisation. .def(py::init&, const std::string&, bool>(), py::arg("molecule"), py::arg("shells"), py::arg("name") = std::string(""), py::arg("coefficients_pre_normalized") = true, "Build a BasisSet from explicit per-shell metadata. " "Each ShellInfo contributes one libint shell at the supplied " "origin with the given primitives + contraction " "coefficients. Use this to construct modrho-renormalised aux " "bases or any other custom-coefficient basis without " "round-tripping through a .g94 file. " "coefficients_pre_normalized=True (default) takes the " "coefficients as-given (the natural form for round-tripping " "shells extracted via .shells()); set False to have libint " "embed primitive normalisation factors (the Gaussian / .g94 " "convention).") .def_property_readonly("name", &vibeqc::BasisSet::name) .def_property_readonly("nbasis", &vibeqc::BasisSet::nbasis, "Number of basis functions (n_bf) — sum over shells of " "(2L+1) for spherical / (L+1)(L+2)/2 for Cartesian. " "Property (no parens). Was a method until v0.4 polish.") .def_property_readonly("nshells", &vibeqc::BasisSet::nshells, "Number of contracted shells. Property (no parens). Was " "a method until v0.4 polish.") .def("shells", &vibeqc::BasisSet::shells, "List of ShellInfo structs — one per (shell × contraction) " "pair, in libint's native shell order. Use this to export " "the basis to external formats (molden, NWChem, ...). " "Method, not property — returns a freshly-built list."); // ----- Integrals ----------------------------------------------------- // Eigen matrices are returned; pybind11 delivers them as NumPy arrays. // Release the GIL for the duration of the native compute — Python isn't // doing anything during integral evaluation. m.def("compute_overlap", &vibeqc::compute_overlap, py::arg("basis"), py::call_guard(), "AO overlap matrix S_{mu nu} = ."); m.def("compute_kinetic", &vibeqc::compute_kinetic, py::arg("basis"), py::call_guard(), "AO kinetic-energy matrix T_{mu nu}."); m.def("compute_nuclear", &vibeqc::compute_nuclear, py::arg("basis"), py::arg("molecule"), py::call_guard(), "AO nuclear-attraction matrix V_{mu nu}."); py::class_(m, "DipoleIntegrals", "x, y, z components of the dipole matrix <μ|r_c − O_c|ν>.") .def_readonly("x", &vibeqc::DipoleIntegrals::x) .def_readonly("y", &vibeqc::DipoleIntegrals::y) .def_readonly("z", &vibeqc::DipoleIntegrals::z); m.def("compute_dipole", &vibeqc::compute_dipole, py::arg("basis"), py::arg("origin") = std::array{0.0, 0.0, 0.0}, py::call_guard(), "Dipole-integral matrices <μ|r − O|ν> about ``origin`` (bohr). " "Used to compute dipole moments via " " = -tr(P·M) + Σ Z_A (R_A − O)."); // Compute the dense 4-index ERI tensor. We hand ownership of the // std::vector buffer to NumPy via a capsule so there's no copy on the // way back to Python (important as nbasis grows). m.def("compute_eri", [](const vibeqc::BasisSet& basis) { vibeqc::Eri4D eri; { py::gil_scoped_release unlock_gil; eri = vibeqc::compute_eri(basis); } const auto n = static_cast(eri.n); // Move the vector onto the heap so the capsule can outlive this // scope and own it through NumPy's lifetime. auto* owner = new std::vector(std::move(eri.data)); py::capsule keep_alive(owner, [](void* p) { delete static_cast*>(p); }); return py::array_t( {n, n, n, n}, {static_cast(n * n * n * sizeof(double)), static_cast(n * n * sizeof(double)), static_cast(n * sizeof(double)), static_cast(sizeof(double))}, owner->data(), keep_alive); }, py::arg("basis"), "AO electron-repulsion integral tensor (mu nu | lambda sigma), " "chemists' notation. Shape (n, n, n, n) with n = nbasis."); // Density-fitting integral kernels. Math + references in // cpp/include/vibeqc/df.hpp; the Python side // (vibeqc.density_fitting.DensityFitting) consumes these and owns // the Cholesky factorisation, B-tensor, and J / K / MO contractions. m.def("compute_2c_eri", &vibeqc::compute_2c_eri, py::arg("aux"), py::call_guard(), "Two-centre Coulomb metric V_PQ = (P|Q) on an auxiliary basis. " "Returns a symmetric positive-definite (n_aux, n_aux) matrix."); m.def("compute_2c_eri_gradient_weighted", &vibeqc::compute_2c_eri_gradient_weighted, py::arg("aux"), py::arg("molecule"), py::arg("omega"), py::call_guard(), "Per-atom Σ_PQ Ω_PQ ∂V_PQ/∂R for an arbitrary (n_aux, n_aux) " "weight Ω. The DF-J gradient passes -(1/2) γ γ^T; the DF-K " "gradient passes ω = (η^P : η^Q). Returns (n_atoms, 3)."); m.def("compute_2c_eri_gradient_contribution", &vibeqc::compute_2c_eri_gradient_contribution, py::arg("aux"), py::arg("molecule"), py::arg("gamma"), py::call_guard(), "Per-atom Σ_PQ γ_P γ_Q ∂V_PQ/∂R (rank-1 special case of " "compute_2c_eri_gradient_weighted; convenience for the J path)."); m.def("compute_3c_eri_gradient_weighted", &vibeqc::compute_3c_eri_gradient_weighted, py::arg("orbital"), py::arg("aux"), py::arg("molecule"), py::arg("W"), py::call_guard(), "Per-atom Σ_{P, μν} W^P_μν ∂(P|μν)/∂R for an arbitrary " "(n_aux, n_orb*n_orb) row-major weight tensor W. The DF-J " "gradient passes γ_P · D_μν; the DF-K gradient passes " "-2 Y^P_μν. Returns (n_atoms, 3)."); m.def("compute_3c_eri_gradient_contribution", &vibeqc::compute_3c_eri_gradient_contribution, py::arg("orbital"), py::arg("aux"), py::arg("molecule"), py::arg("D"), py::arg("gamma"), py::call_guard(), "Per-atom Σ_P γ_P · Σ_μν D_μν ∂(P|μν)/∂R (convenience for " "the J path; equivalent to " "compute_3c_eri_gradient_weighted with W^P_μν = γ_P · D_μν)."); m.def("compute_df_j_gradient", [](const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux, const vibeqc::Molecule& mol, const Eigen::MatrixXd& D) { vibeqc::DensityFitting df(orbital, aux); return df.compute_j_gradient(mol, D); }, py::arg("orbital"), py::arg("aux"), py::arg("molecule"), py::arg("D"), py::call_guard(), "Assemble the complete DF Coulomb-energy gradient at fixed D: " "Σ_P γ_P ∂ρ_P/∂R - (1/2) Σ_PQ γ_P γ_Q ∂V_PQ/∂R. For pure DFT " "(α_HF = 0) this is the *complete* DF analytic two-electron " "gradient. For HF / hybrid DFT pair with the K-gradient " "contribution."); m.def("compute_df_k_gradient", [](const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux, const vibeqc::Molecule& mol, const Eigen::MatrixXd& C_occ, double alpha_hf) { vibeqc::DensityFitting df(orbital, aux); return df.compute_k_gradient(mol, C_occ, alpha_hf); }, py::arg("orbital"), py::arg("aux"), py::arg("molecule"), py::arg("C_occ"), py::arg("alpha_hf") = 1.0, py::call_guard(), "Assemble the DF Exchange-energy gradient at fixed (C_occ, " "α_HF). E_K = -(α_HF/4) tr(D · K_DF) with D = 2 C_occ C_occ^T " "(closed shell). For pure DFT (α_HF = 0) returns zero."); m.def("compute_df_jk_gradient", [](const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux, const vibeqc::Molecule& mol, const Eigen::MatrixXd& D, const Eigen::MatrixXd& C_occ, double alpha_hf) { vibeqc::DensityFitting df(orbital, aux); return df.compute_jk_gradient(mol, D, C_occ, alpha_hf); }, py::arg("orbital"), py::arg("aux"), py::arg("molecule"), py::arg("D"), py::arg("C_occ"), py::arg("alpha_hf") = 1.0, py::call_guard(), "Combined DF J + α_HF·K gradient — the full DF analytic " "two-electron gradient for HF / hybrid DFT closed-shell. " "Folds both contractions into one 2c + 3c kernel call each, " "halving the libint derivative work compared to " "compute_df_j_gradient + compute_df_k_gradient."); m.def("build_cosx_q", &vibeqc::build_cosx_q, py::arg("basis"), py::arg("cosx_grid"), py::call_guard(), "Build the basis-only overlap-fit Q-junction matrix used by " "the COSX K corrector (Neese 2009 §2.4). Q = S · S_grid^-1 " "is a function of the AO basis and the cosx grid only — " "compute it once at SCF setup and pass it to every " "``compute_cosx_k`` call to skip the per-iter Q rebuild. " "Returns 0×0 on LDLT failure of S_grid (signal to fall back " "to the uncorrected K_naive)."); m.def("build_cosx_schwarz", &vibeqc::build_cosx_schwarz, py::arg("basis"), py::call_guard(), "Build the basis-only Schwarz upper-bound table " "Q(s1, s2) = sqrt(max |(μν|μν)|) for the COSX shell-pair " "screen. Symmetric (n_shells, n_shells). Pass to " "``compute_cosx_k`` as ``schwarz_cached`` to enable the " "tighter screen: drops pairs whose contribution to K is " "below tolerance before paying for the libint " "nuclear-attraction call. Basis-only and SCF-invariant."); m.def("compute_shell_radial_cutoffs", [](const vibeqc::BasisSet& basis, double tol) { return vibeqc::compute_shell_radial_cutoffs( basis.libint(), tol); }, py::arg("basis"), py::arg("tol") = 1e-12, "Per-shell radial extent (bohr): smallest r* such that " "|χ_μ(r)| < tol for any μ in the shell and |r - O_s| > r*. " "Conservative (Gaussian-extent) bound derived from each " "shell's slowest-decaying primitive. Pass to " "``compute_cosx_k`` as ``shell_cutoffs`` to enable the " "active-shell pre-prune in the K-build pair loop."); // --- GridBatch / GridBatches: per-batch view of the integration grid. // The structs are exposed read-only — they are infrastructure for // COSX-K / batched XC builders, not a user-facing API. Surfacing // them in Python is needed only to validate the batching from // ``tests/test_grid_batch.py``. py::class_(m, "GridBatch") .def_readonly("start_index", &vibeqc::GridBatch::start_index) .def_readonly("n_points", &vibeqc::GridBatch::n_points) .def_readonly("points", &vibeqc::GridBatch::points) .def_readonly("weights", &vibeqc::GridBatch::weights) .def_readonly("primary_shells", &vibeqc::GridBatch::primary_shells) .def_readonly("primary_bfs", &vibeqc::GridBatch::primary_bfs) .def_readonly("chi_primary", &vibeqc::GridBatch::chi_primary) .def_readonly("has_gradient", &vibeqc::GridBatch::has_gradient); py::class_(m, "GridBatches") .def_readonly("batches", &vibeqc::GridBatches::batches) .def_readonly("n_pts_total", &vibeqc::GridBatches::n_pts_total) .def_readonly("n_bf_total", &vibeqc::GridBatches::n_bf_total) .def_readonly("batch_size_hint", &vibeqc::GridBatches::batch_size_hint); m.def("build_grid_batches", &vibeqc::build_grid_batches, py::arg("basis"), py::arg("grid"), py::arg("batch_size"), py::arg("shell_cutoffs"), py::arg("need_gradient") = false, py::call_guard(), "Build a batched view of an integration grid with per-batch " "primary-shell pruning + a basis-only chi cache. Returns a " "``GridBatches`` carrying ``batches[i].chi_primary`` already " "evaluated. ``shell_cutoffs`` should come from " "``compute_shell_radial_cutoffs(basis, AO_TOL)``. SCF-" "invariant — build once at setup; reuse across iterations."); m.def("compute_robust_cosx_k", [](const Eigen::MatrixXd& D, const Eigen::MatrixXd& K_XvX, const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux) { vibeqc::DensityFitting df(orbital, aux); return vibeqc::compute_robust_cosx_k(D, K_XvX, df); }, py::arg("D"), py::arg("K_XvX"), py::arg("orbital"), py::arg("aux"), py::call_guard(), "Robust Dunlap-fit COSX exchange. K_robust = K_XvX + K_QvQ " "− K_XvQ − K_QvX where K_XvQ = Σ_P B^P · K_XvX · B^P. " "Constructs a C++ DensityFitting from orbital+aux internally " "for the B-tensor. (Neese 2009 §2.5)."); m.def("compute_cosx_k", [](const vibeqc::BasisSet& basis, const Eigen::MatrixXd& density, const vibeqc::Grid& cosx_grid, const Eigen::MatrixXd& q_cached, const Eigen::MatrixXd& schwarz_cached, const std::vector& shell_cutoffs, const vibeqc::GridBatches* grid_batches) { return vibeqc::compute_cosx_k( basis, density, cosx_grid, q_cached, schwarz_cached, shell_cutoffs, grid_batches); }, py::arg("basis"), py::arg("density"), py::arg("cosx_grid"), py::arg("q_cached") = Eigen::MatrixXd(), py::arg("schwarz_cached") = Eigen::MatrixXd(), py::arg("shell_cutoffs") = std::vector(), py::arg("grid_batches") = static_cast(nullptr), py::call_guard(), "Seminumerical exchange matrix K via Neese's chain-of-spheres " "(COSX) algorithm. K_{μν} ≈ Σ_g w_g χ_μ(r_g) (A(r_g)·D·χ(r_g))_ν " "where A(r_g) is the analytic 1/|r-r_g| Coulomb-attraction " "matrix at the grid point. Returns the symmetrised K. Pair " "with RI-J for the standard RIJCOSX hybrid-DFT acceleration. " "Optional basis-only / SCF-invariant caches — build once at " "setup: ``q_cached`` from ``build_cosx_q`` (skip Q rebuild); " "``schwarz_cached`` from ``build_cosx_schwarz`` (tighter " "shell-pair screen); ``shell_cutoffs`` from " "``compute_shell_radial_cutoffs`` (active-shell pre-prune)."); m.def("compute_cosx_k_gradient_contribution", &vibeqc::compute_cosx_k_gradient_contribution, py::arg("molecule"), py::arg("basis"), py::arg("density"), py::arg("cosx_grid"), py::arg("alpha_hf"), py::arg("q_cached") = Eigen::MatrixXd(), py::call_guard(), "Analytic nuclear gradient of the COSX K-energy term: " "dE_K/dR with E_K = -(α_hf/2) tr(D · K_cosx). Frozen " "grid / weights / Q-fit; per-atom error ~µHa/bohr on neutral " "organics with the default sparse grid. Pass the cached Q " "matrix from the K-build to skip redundant Q recomputation."); // --- Test-only entry points for the in-tree COSX nuclear- // attraction kernel. These exercise the kernel one shell pair at // a time so tests can pin it against libint at ULP. Not part of // the user-facing API — prefix with underscore. m.def("_cosx_nuclear_pair_custom", [](const vibeqc::BasisSet& basis, int s1, int s2, const std::array& C) { vibeqc::ensure_libint_initialized(); const auto& shells = basis.libint(); const int n_shells = static_cast(shells.size()); if (s1 < 0 || s1 >= n_shells || s2 < 0 || s2 >= n_shells) { throw std::out_of_range( "_cosx_nuclear_pair_custom: shell index out of range"); } auto boys = vibeqc::build_boys_table(shells.max_l()); auto cache = vibeqc::build_primitive_pair_cache(shells); const auto& pp = cache.pairs[ static_cast(s1) * n_shells + s2]; return vibeqc::cosx_nuclear_pair( shells[s1], shells[s2], pp, C, boys); }, py::arg("basis"), py::arg("s1"), py::arg("s2"), py::arg("C"), py::call_guard(), "Test-only: in-tree nuclear-attraction kernel on a single " "shell pair (s1, s2) with one pseudo-nucleus of charge q = −1 " "at C (bohr). Returns the (n1, n2) contracted block. Throws " "on non-(s, s) shells until the general Obara-Saika path " "lands (kernel arc commit 2)."); m.def("_cosx_nuclear_pair_libint", [](const vibeqc::BasisSet& basis, int s1, int s2, const std::array& C) { vibeqc::ensure_libint_initialized(); const auto& shells = basis.libint(); const int n_shells = static_cast(shells.size()); if (s1 < 0 || s1 >= n_shells || s2 < 0 || s2 >= n_shells) { throw std::out_of_range( "_cosx_nuclear_pair_libint: shell index out of range"); } libint2::Engine engine(libint2::Operator::nuclear, shells.max_nprim(), shells.max_l(), 0); std::vector>> charges{{-1.0, C}}; engine.set_params(charges); engine.compute(shells[s1], shells[s2]); const auto n1 = static_cast(shells[s1].size()); const auto n2 = static_cast(shells[s2].size()); Eigen::Matrix out = Eigen::Matrix::Zero(n1, n2); const double* block = engine.results()[0]; if (block != nullptr) { std::memcpy(out.data(), block, sizeof(double) * static_cast(n1) * static_cast(n2)); } return out; }, py::arg("basis"), py::arg("s1"), py::arg("s2"), py::arg("C"), py::call_guard(), "Test-only: libint Operator::nuclear reference for the same " "(s1, s2, C, q = −1) input as ``_cosx_nuclear_pair_custom``. " "Used by tests/test_cosx_kernel.py to pin the in-tree kernel " "to ULP against libint."); m.def("compute_3c_eri", [](const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux) { vibeqc::Eri3D T; { py::gil_scoped_release unlock_gil; T = vibeqc::compute_3c_eri(orbital, aux); } const auto na = static_cast(T.n_aux); const auto no = static_cast(T.n_orb); auto* owner = new std::vector(std::move(T.data)); py::capsule keep_alive(owner, [](void* p) { delete static_cast*>(p); }); return py::array_t( {na, no, no}, {static_cast(no * no * sizeof(double)), static_cast(no * sizeof(double)), static_cast(sizeof(double))}, owner->data(), keep_alive); }, py::arg("orbital"), py::arg("aux"), "Three-centre ERI tensor (P | mu nu) with shape " "(n_aux, n_orb, n_orb). Symmetric in (mu, nu); both off-diagonal " "positions are filled so callers iterate the full extent."); // ----- Periodic DF kernels ------------------------------------------ // Periodic 2c / 3c Coulomb integrals for the GDF Lpq build. See // cpp/include/vibeqc/aux_eri.hpp for the math contract; the Python // consumer is vibeqc.aux_basis.build_lpq_native. m.def("compute_2c_eri_lattice", &vibeqc::compute_2c_eri_lattice, py::arg("aux"), py::arg("system"), py::arg("opts"), py::call_guard(), "Periodic 2-centre Coulomb metric M_PQ = Σ_T (P_0 | Q_T) on an " "auxiliary basis. Image sum runs over direct-lattice cells with " "|g| ≤ opts.cutoff_bohr. Returns a symmetric (n_aux, n_aux) " "matrix. Bare lattice sum — for general convergence on " "non-charge-compensated aux bases use a charge-compensated aux " "basis (see vibeqc.aux_basis.make_aux_basis_set)."); m.def("compute_2c_eri_lattice_blocks", [](const vibeqc::BasisSet& aux, const vibeqc::PeriodicSystem& system, const vibeqc::LatticeSumOptions& opts) { vibeqc::Lattice2CEriBlockSet blocks; { py::gil_scoped_release unlock_gil; blocks = vibeqc::compute_2c_eri_lattice_blocks( aux, system, opts); } const auto nc = static_cast(blocks.cells.size()); const auto na = static_cast(blocks.n_aux); py::array_t indices({nc, py::ssize_t{3}}); py::array_t vectors({nc, py::ssize_t{3}}); py::array_t data({nc, na, na}); auto idx = indices.mutable_unchecked<2>(); auto vec = vectors.mutable_unchecked<2>(); auto arr = data.mutable_unchecked<3>(); for (py::ssize_t c = 0; c < nc; ++c) { const auto& cell = blocks.cells[static_cast(c)]; for (py::ssize_t a = 0; a < 3; ++a) { idx(c, a) = cell.index[static_cast(a)]; vec(c, a) = cell.r_cart[static_cast(a)]; } const auto& M = blocks.blocks[static_cast(c)]; for (py::ssize_t p = 0; p < na; ++p) { for (py::ssize_t q = 0; q < na; ++q) { arr(c, p, q) = M( static_cast(p), static_cast(q)); } } } return py::make_tuple(indices, vectors, data); }, py::arg("aux"), py::arg("system"), py::arg("opts"), "Cell-resolved periodic 2-centre Coulomb metric blocks. " "Returns (cell_indices, cell_vectors_bohr, blocks), where " "blocks[c, P, Q] = (P_0 | Q_Tc). Summing over c reproduces " "compute_2c_eri_lattice at Γ."); m.def("compute_3c_eri_lattice", [](const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux, const vibeqc::PeriodicSystem& system, const vibeqc::LatticeSumOptions& opts) { vibeqc::Eri3D T; { py::gil_scoped_release unlock_gil; T = vibeqc::compute_3c_eri_lattice(orbital, aux, system, opts); } const auto na = static_cast(T.n_aux); const auto no = static_cast(T.n_orb); auto* owner = new std::vector(std::move(T.data)); py::capsule keep_alive(owner, [](void* p) { delete static_cast*>(p); }); return py::array_t( {na, no, no}, {static_cast(no * no * sizeof(double)), static_cast(no * sizeof(double)), static_cast(sizeof(double))}, owner->data(), keep_alive); }, py::arg("orbital"), py::arg("aux"), py::arg("system"), py::arg("opts"), "Periodic 3-centre ERI tensor T_{P, μ, ν} = Σ_T (P_0 | μ_0 ν_T). " "Aux P and orbital μ stay anchored at the reference cell; the " "second orbital ν is shifted across direct-lattice translations " "with |g| ≤ opts.cutoff_bohr. Symmetric in (μ, ν). Shape " "(n_aux, n_orb, n_orb)."); m.def("compute_3c_eri_lattice_blocks", [](const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux, const vibeqc::PeriodicSystem& system, const vibeqc::LatticeSumOptions& opts) { vibeqc::Lattice3CEriBlockSet blocks; { py::gil_scoped_release unlock_gil; blocks = vibeqc::compute_3c_eri_lattice_blocks( orbital, aux, system, opts); } const auto nc = static_cast(blocks.cells.size()); const auto na = static_cast(blocks.n_aux); const auto no = static_cast(blocks.n_orb); py::array_t indices({nc, py::ssize_t{3}}); py::array_t vectors({nc, py::ssize_t{3}}); py::array_t data({nc, na, no, no}); auto idx = indices.mutable_unchecked<2>(); auto vec = vectors.mutable_unchecked<2>(); auto arr = data.mutable_unchecked<4>(); for (py::ssize_t c = 0; c < nc; ++c) { const auto& cell = blocks.cells[static_cast(c)]; for (py::ssize_t a = 0; a < 3; ++a) { idx(c, a) = cell.index[static_cast(a)]; vec(c, a) = cell.r_cart[static_cast(a)]; } const auto& T = blocks.blocks[static_cast(c)]; for (py::ssize_t p = 0; p < na; ++p) { for (py::ssize_t mu = 0; mu < no; ++mu) { for (py::ssize_t nu = 0; nu < no; ++nu) { arr(c, p, mu, nu) = T( static_cast(p), static_cast(mu), static_cast(nu)); } } } } return py::make_tuple(indices, vectors, data); }, py::arg("orbital"), py::arg("aux"), py::arg("system"), py::arg("opts"), "Cell-resolved periodic 3-centre ERI blocks. Returns " "(cell_indices, cell_vectors_bohr, blocks), where " "blocks[c, P, mu, nu] = (P_0 | mu_0 nu_Tc). Individual " "blocks are not symmetrized in (mu,nu); summing over c " "reproduces compute_3c_eri_lattice at Γ."); // RSGDF — short-range halves (Ye & Berkelbach, J. Chem. Phys. 154, // 131104 (2021), DOI 10.1063/5.0046617). Same kernels as above but // with the Coulomb operator replaced by erfc(ω r)/r. m.def("compute_2c_eri_lattice_sr", &vibeqc::compute_2c_eri_lattice_sr, py::arg("aux"), py::arg("system"), py::arg("opts"), py::arg("omega"), py::call_guard(), "RSGDF short-range 2-centre metric " "M^SR_{PQ}(ω) = Σ_T (P_0 | erfc(ωr)/r | Q_T). The lattice sum " "converges absolutely (no compensating-charge bookkeeping " "needed) because the erfc kernel decays as exp(-ω²r²). The " "long-range half is built in Python from analytical Gaussian " "Fourier transforms — see vibeqc.aux_basis.build_lpq_rsgdf."); m.def("compute_3c_eri_lattice_sr", [](const vibeqc::BasisSet& orbital, const vibeqc::BasisSet& aux, const vibeqc::PeriodicSystem& system, const vibeqc::LatticeSumOptions& opts, double omega) { vibeqc::Eri3D T; { py::gil_scoped_release unlock_gil; T = vibeqc::compute_3c_eri_lattice_sr( orbital, aux, system, opts, omega); } const auto na = static_cast(T.n_aux); const auto no = static_cast(T.n_orb); auto* owner = new std::vector(std::move(T.data)); py::capsule keep_alive(owner, [](void* p) { delete static_cast*>(p); }); return py::array_t( {na, no, no}, {static_cast(no * no * sizeof(double)), static_cast(no * sizeof(double)), static_cast(sizeof(double))}, owner->data(), keep_alive); }, py::arg("orbital"), py::arg("aux"), py::arg("system"), py::arg("opts"), py::arg("omega"), "RSGDF short-range 3-centre tensor " "T^SR_{P, μν}(ω) = Σ_T (P_0 | erfc(ωr)/r | μ_0 ν_T). " "Image-summed on the second orbital factor; converges " "absolutely thanks to the erfc decay. Shape (n_aux, n_orb, " "n_orb), symmetric in (μ, ν)."); // ----- RHF ----------------------------------------------------------- py::enum_(m, "InitialGuess") .value("AUTO", vibeqc::InitialGuess::AUTO) .value("HCORE", vibeqc::InitialGuess::HCORE) .value("SAD", vibeqc::InitialGuess::SAD) .value("SAP", vibeqc::InitialGuess::SAP) .value("PATOM", vibeqc::InitialGuess::PATOM) .value("HUECKEL", vibeqc::InitialGuess::HUECKEL) .value("MINAO", vibeqc::InitialGuess::MINAO) .value("READ", vibeqc::InitialGuess::READ); // ---- Engine helpers used by the Python periodic drivers ------------ // See python/vibeqc/guess.py for the Python-facing wrappers. These // private bindings call into GuessEngine with S/Hcore/JK = nullptr — // matching what the Python periodic drivers can supply today. The // future Fock-mode guesses (SAP, HUECKEL) will need an overload // that takes Hcore + a periodic JKBuilder; out of scope for v0.9.x // Commit 1, which is the scaffolding-only step. m.def("_guess_closed_shell_density", [](const vibeqc::Molecule& mol, const vibeqc::BasisSet& basis, int n_occ, vibeqc::InitialGuess kind, bool is_periodic) -> py::object { vibeqc::SystemHints hints; hints.is_periodic = is_periodic; auto g = vibeqc::GuessEngine::build_closed_shell( mol, basis, n_occ, kind, /*S=*/nullptr, /*Hcore=*/nullptr, /*jk=*/nullptr, hints); if (g.D.size() == 0) return py::none(); return py::cast(g.D); }, py::arg("mol"), py::arg("basis"), py::arg("n_occ"), py::arg("kind"), py::arg("is_periodic") = true, "Closed-shell initial density via GuessEngine. Returns None " "for HCORE (and any future Fock-mode kind reached without a " "JKBuilder)."); m.def("_guess_open_shell_density", [](const vibeqc::Molecule& mol, const vibeqc::BasisSet& basis, int n_alpha, int n_beta, vibeqc::InitialGuess kind, bool is_periodic) -> py::object { vibeqc::SystemHints hints; hints.is_periodic = is_periodic; hints.is_open_shell = true; auto g = vibeqc::GuessEngine::build_open_shell( mol, basis, n_alpha, n_beta, kind, /*S=*/nullptr, /*Hcore=*/nullptr, /*jk=*/nullptr, hints); if (g.D_alpha.size() == 0) return py::none(); return py::make_tuple(g.D_alpha, g.D_beta); }, py::arg("mol"), py::arg("basis"), py::arg("n_alpha"), py::arg("n_beta"), py::arg("kind"), py::arg("is_periodic") = true, "Open-shell per-spin initial densities via GuessEngine. " "Returns None for HCORE."); py::class_(m, "DIIS", "Pulay's Direct Inversion of the Iterative Subspace SCF " "accelerator (Chem. Phys. Lett. 73, 393 (1980); J. Comput. " "Chem. 3, 556 (1982)). Maintains a rolling history of " "(Fock, error) pairs; on each call solves the small (n+1)-" "dim Pulay linear system for coefficients c_i that minimise " "‖Σ c_i e_i‖ subject to Σ c_i = 1, and returns the " "extrapolated Fock Σ c_i F_i. Used internally by every " "molecular and periodic SCF driver — this Python binding " "exists so periodic Python SCF loops can call into the " "same canonical implementation instead of re-deriving " "Pulay's recurrence in numpy.") .def(py::init(), py::arg("max_subspace") = 8, "Construct an empty DIIS with the given history cap " "(must be ≥ 2; default 8).") .def("clear", &vibeqc::DIIS::clear, "Drop all (F, error) history. Used for level-shift " "warm-up restarts and similar SCF resets.") .def_property_readonly("subspace_size", &vibeqc::DIIS::subspace_size, "Current history size (0 immediately after construction " "or ``clear()``; capped at ``max_subspace``).") .def("extrapolate", &vibeqc::DIIS::extrapolate, py::arg("fock"), py::arg("error"), "Append ``(fock, error)`` to the history and return the " "Pulay-extrapolated Fock Σ_i c_i F_i. Returns ``fock`` " "unchanged before the history has at least two entries. " "``error`` is conventionally the commutator " "``F D S − S D F``, evaluated at the density that " "produced ``fock``."); py::enum_(m, "SCFAccelerator", "SCF Fock-extrapolation accelerator. " "DIIS = Pulay's commutator DIIS (historical default); " "KDIIS = Kollmar 2005 (orbital-rotation gradient error vector; " "ORCA's strict-convergence default); " "EDIIS = energy-DIIS (Kudin / Scuseria / Cancès, J. Chem. " "Phys. 116, 8255 (2002)); " "EDIIS_DIIS = production hybrid (Garza / Scuseria, J. Chem. " "Phys. 137, 054110 (2012)) — EDIIS while the commutator " "norm is above ``ediis_diis_switch_threshold``, then plain " "DIIS for the asymptotic regime; " "ADIIS = augmented-Roothaan-Hall DIIS (Hu & Yang, J. Chem. " "Phys. 132, 054109 (2010)) — EDIIS sibling, near-identical " "at the HF level.") .value("DIIS", vibeqc::SCFAccelerator::DIIS) .value("KDIIS", vibeqc::SCFAccelerator::KDIIS) .value("EDIIS", vibeqc::SCFAccelerator::EDIIS) .value("EDIIS_DIIS", vibeqc::SCFAccelerator::EDIIS_DIIS) .value("ADIIS", vibeqc::SCFAccelerator::ADIIS); py::enum_(m, "SCFMode", "SCF Fock-build mode. CONVENTIONAL = in-core 4-index ERI " "tensor (fast for ≤~150 BF, O(n_bf⁴) memory). DIRECT = " "on-the-fly Schwarz-screened libint quartet evaluation " "(the only path that survives at >250 BF). AUTO (default) " "selects on basis size (DIRECT when n_bf > " "scf_mode_auto_threshold). Orthogonal to density_fit + " "cosx — those supersede it when enabled.") .value("AUTO", vibeqc::SCFMode::AUTO) .value("CONVENTIONAL", vibeqc::SCFMode::CONVENTIONAL) .value("DIRECT", vibeqc::SCFMode::DIRECT); py::class_(m, "EDIIS", "Energy-DIIS extrapolator (Kudin / Scuseria / Cancès, " "J. Chem. Phys. 116, 8255 (2002)). Minimises a quadratic energy " "functional on the convex hull of stored (F, D) pairs subject " "to Σ c_i = 1, c_i ≥ 0 — a small simplex-constrained QP. The " "SCF drivers use this internally when scf_accelerator = EDIIS " "or EDIIS_DIIS; the class is exposed primarily for testing " "the coefficient solver against hand-constructed inputs.") .def(py::init(), py::arg("max_subspace") = 8) .def("clear", &vibeqc::EDIIS::clear) .def("subspace_size", &vibeqc::EDIIS::subspace_size) .def("last_coeffs", [](const vibeqc::EDIIS& e) { return e.last_coeffs(); }, "Coefficients from the most recent extrapolate() call. " "Empty before the first call. Sum to 1, all ≥ 0.") .def("extrapolate", py::overload_cast(&vibeqc::EDIIS::extrapolate), py::arg("fock"), py::arg("density"), py::arg("energy"), "Closed-shell extrapolation: store (F, D, energy) and " "return Σ_i c_i F_i. Returns F unchanged when n < 2.") .def("extrapolate_uhf", py::overload_cast(&vibeqc::EDIIS::extrapolate), py::arg("fock_alpha"), py::arg("fock_beta"), py::arg("density_alpha"), py::arg("density_beta"), py::arg("energy"), "Open-shell extrapolation: store the four matrices + total " "energy and return (Σ c_i F_α,i, Σ c_i F_β,i). The same " "coefficient set applies to both spins.") .def("extrapolate_blocks", py::overload_cast&, const std::vector&, double>(&vibeqc::EDIIS::extrapolate), py::arg("fock_blocks"), py::arg("density_blocks"), py::arg("energy"), "Block-vector extrapolation. Stores an (F, D, energy) " "iterate whose Fock and density are each a list of " "nbf×nbf real blocks and returns Σ_i c_i F_i block-by-" "block. The QP cross-term ⟨F_i | D_j⟩ = Σ_b (F_i[b] ⊙ " "D_j[b]).sum() sums over all blocks. Used by multi-k " "periodic SCF (blocks = real-space cells of a " "LatticeMatrixSet) and by the open-shell UHF kernel " "(blocks = α + β). Block layout (count, ordering, " "dimensions) must match across pushed iterates. Returns " "``fock_blocks`` unchanged when n < 2."); py::class_(m, "ADIIS", "Augmented-Roothaan-Hall DIIS extrapolator (Hu & Yang, " "J. Chem. Phys. 132, 054109 (2010)). Minimises the ARH " "energy model expanded about the most recent iterate, on the " "convex hull of stored (F, D) pairs subject to Σ c_i = 1, " "c_i ≥ 0 — the same simplex-constrained QP solver as EDIIS. " "Needs no per-iterate energy. The SCF drivers use this when " "scf_accelerator = ADIIS; the class is exposed primarily for " "testing the coefficient solver against hand-constructed " "inputs.") .def(py::init(), py::arg("max_subspace") = 8) .def("clear", &vibeqc::ADIIS::clear) .def("subspace_size", &vibeqc::ADIIS::subspace_size) .def("last_coeffs", [](const vibeqc::ADIIS& a) { return a.last_coeffs(); }, "Coefficients from the most recent extrapolate() call. " "Empty before the first call. Sum to 1, all ≥ 0.") .def("extrapolate", py::overload_cast( &vibeqc::ADIIS::extrapolate), py::arg("fock"), py::arg("density"), "Closed-shell extrapolation: store (F, D) and return " "Σ_i c_i F_i. Returns F unchanged when n < 2.") .def("extrapolate_uhf", py::overload_cast( &vibeqc::ADIIS::extrapolate), py::arg("fock_alpha"), py::arg("fock_beta"), py::arg("density_alpha"), py::arg("density_beta"), "Open-shell extrapolation: store the four matrices and " "return (Σ c_i F_α,i, Σ c_i F_β,i). The same coefficient " "set applies to both spins.") .def("extrapolate_blocks", py::overload_cast&, const std::vector&>( &vibeqc::ADIIS::extrapolate), py::arg("fock_blocks"), py::arg("density_blocks"), "Block-vector extrapolation. Stores an (F, D) iterate " "whose Fock and density are each a list of nbf×nbf real " "blocks and returns Σ_i c_i F_i block-by-block. The ARH " "trace cross-term ⟨D_i − D_n | F_j − F_n⟩ = Σ_b " "((D_i − D_n)[b] ⊙ (F_j − F_n)[b]).sum() sums over all " "blocks. Used by multi-k periodic SCF (blocks = real-" "space cells) and by the open-shell UHF kernel (blocks " "= α + β). Block layout must match across pushed " "iterates. Returns ``fock_blocks`` unchanged when " "n < 2."); py::class_(m, "KDIIS", "Kollmar's DIIS — orbital-rotation-gradient DIIS (C. Kollmar, " "Int. J. Quantum Chem. 105, 685 (2005); ORCA's strict-" "convergence default). Replaces the AO-basis commutator error " "e = F·D·S − S·D·F with the canonical-MO-basis occ-vir block " "g_{ai} = (C^T F C)_{ai}. The Pulay least-squares extrapolation " "is identical (Σ c_i F_i with Σ c_i = 1 minimising " "‖Σ c_i e_i‖_F); only the error metric differs. Particularly " "robust on transition-metal complexes and open-shell systems " "where the AO-commutator error is dominated by the wrong " "eigenmodes. The SCF drivers use this when " "scf_accelerator = KDIIS; the class is exposed primarily for " "testing and for Python periodic SCF loops to call into the " "canonical implementation instead of re-deriving the recurrence " "in numpy.") .def(py::init(), py::arg("max_subspace") = 8) .def("clear", &vibeqc::KDIIS::clear) .def_property_readonly("subspace_size", &vibeqc::KDIIS::subspace_size) .def("extrapolate", py::overload_cast(&vibeqc::KDIIS::extrapolate), py::arg("fock"), py::arg("mo_coeffs"), py::arg("mo_energies"), py::arg("n_occ"), "Closed-shell extrapolation. Builds g_{ai} = (C^T F C)_{ai} " "(occ-vir block in the canonical MO basis), pushes (F, g) " "onto the rolling history, solves the Pulay least-squares, " "and returns the extrapolated Fock in the AO basis. Returns " "``fock`` unchanged before the history has at least two " "entries.") .def("extrapolate_uhf", py::overload_cast(&vibeqc::KDIIS::extrapolate), py::arg("fock_alpha"), py::arg("fock_beta"), py::arg("mo_coeffs_alpha"), py::arg("mo_coeffs_beta"), py::arg("mo_energies_alpha"), py::arg("mo_energies_beta"), py::arg("n_alpha"), py::arg("n_beta"), "Open-shell extrapolation. Builds per-spin orbital-gradient " "errors, concatenates them into a single error matrix for " "the Pulay B-matrix solve, and returns the extrapolated " "(F_α, F_β) pair. The same coefficient set applies to both " "spins (spin-coupled extrapolation, matching the EDIIS/ADIIS " "open-shell convention)."); m.def("sad_density", &vibeqc::sad_density, py::arg("molecule"), py::arg("basis"), "Superposition of atomic densities (SAD) initial-guess density " "matrix. Builds an isolated-atom RHF SCF (with fractional Aufbau " "occupations) for each unique element in the molecule, then " "assembles the full molecular density matrix block-diagonally " "by AO range. The atomic SCFs are cached per (Z, basis_name) " "for the lifetime of the call. Returns a (n_bf, n_bf) numpy " "array. tr(D · S) ≈ n_electrons. Standard reference: " "Van Lenthe et al., J. Comput. Chem. 27, 926 (2006)."); py::class_(m, "RHFOptions") .def(py::init<>()) .def_readwrite("max_iter", &vibeqc::RHFOptions::max_iter) .def_readwrite("conv_tol_energy", &vibeqc::RHFOptions::conv_tol_energy) .def_readwrite("conv_tol_grad", &vibeqc::RHFOptions::conv_tol_grad) .def_readwrite("damping", &vibeqc::RHFOptions::damping) .def_readwrite("dynamic_damping", &vibeqc::RHFOptions::dynamic_damping, "Adaptive density-mixing α (Zerner-Hehenberger 1979). " "When True, ``damping`` is the initial α; the SCF " "loop adjusts it within " "[``dynamic_damping_min``, ``dynamic_damping_max``] " "based on energy decrease. Default False (static).") .def_readwrite("dynamic_damping_min", &vibeqc::RHFOptions::dynamic_damping_min, "Lower bound for dynamic-damping α. Default 0.0.") .def_readwrite("dynamic_damping_max", &vibeqc::RHFOptions::dynamic_damping_max, "Upper bound for dynamic-damping α. Default 0.95.") .def_readwrite("fock_mixing", &vibeqc::RHFOptions::fock_mixing, "Fock/Kohn-Sham matrix mixing fraction. This is " "the weight of the previous Fock matrix in the " "matrix diagonalised to form the next density. " "CRYSTAL FMIXING 30 corresponds to 0.30. " "Default 0.0 (off).") .def_readwrite("use_diis", &vibeqc::RHFOptions::use_diis) .def_readwrite("diis_start_iter", &vibeqc::RHFOptions::diis_start_iter) .def_readwrite("diis_subspace_size", &vibeqc::RHFOptions::diis_subspace_size) .def_readwrite("scf_accelerator", &vibeqc::RHFOptions::scf_accelerator, "SCF Fock-extrapolation accelerator. See " "SCFAccelerator. Default DIIS (back-compat); " "EDIIS_DIIS is the production hybrid " "recommended by Garza / Scuseria 2012 for " "transition metals and broken-symmetry cases.") .def_readwrite("ediis_diis_switch_threshold", &vibeqc::RHFOptions::ediis_diis_switch_threshold, "EDIIS → DIIS commutator-norm threshold for " "SCFAccelerator.EDIIS_DIIS. Use EDIIS while " "‖F D S − S D F‖_F is above this value; switch " "to plain DIIS once below. Default 1e-1.") .def_readwrite("initial_guess", &vibeqc::RHFOptions::initial_guess) .def_readwrite("linear_dep_threshold", &vibeqc::RHFOptions::linear_dep_threshold) .def_readwrite("ecp_centers", &vibeqc::RHFOptions::ecp_centers, "List[ECPCenter] — atoms carrying an effective " "core potential. Empty (default) = all-electron " "calculation. See Phase 14c.") .def_readwrite("ecp_library", &vibeqc::RHFOptions::ecp_library, "ECP XML library name (default 'ecp10mdf' when " "ecp_centers is non-empty and this is empty).") .def_readwrite("level_shift", &vibeqc::RHFOptions::level_shift, "Phase C1a-2 — Saunders-Hillier level shift " "(Hartree). Adds F + b·S − (b/2)·S·D·S to F " "before diagonalisation, raising virtual orbital " "eigenvalues by b. Inert at the converged " "density (SCF fixed point unchanged). Useful " "for small-HOMO-LUMO-gap molecules where DIIS " "oscillates. Default 0.0 (no shift). Typical " "values: 0.1 – 0.5 Hartree.") .def_readwrite("quadratic_fallback_iter", &vibeqc::RHFOptions::quadratic_fallback_iter, "Phase C1c — second-order SCF fallback " "activation iteration. When > 0, after this " "iteration RHF switches from 'diagonalize F' to " "a Newton step in MO space with diagonal-Hessian " "preconditioning. Use for small-gap systems " "where DIIS + level shift fail to converge. " "Default 0 (disabled).") .def_readwrite("quadratic_fallback_shift", &vibeqc::RHFOptions::quadratic_fallback_shift, "Damping shift λ in the orbital-rotation " "denominator (ε_a − ε_i + λ) of the C1c " "Newton step. Default 0.1 Hartree. Larger " "values produce more conservative steps.") .def_readwrite("quadratic_fallback_max_step", &vibeqc::RHFOptions::quadratic_fallback_max_step, "Trust-region cap on ‖κ‖ per C1c Newton " "step. Default 0.1. Larger values produce " "more aggressive rotations; smaller values " "are safer.") .def_readwrite("newton_threshold", &vibeqc::RHFOptions::newton_threshold, "Phase D2c — Newton (full orbital-Hessian Newton " "via preconditioned CG) activation threshold on " "‖F D S − S D F‖_F. When > 0 and the gradient " "drops below this value, the SCF replaces " "diagonalize-F with a Newton step on the " "orbital-rotation manifold. 1.0 (Fischer-Almlöf " "1992 convention) is a reasonable production " "default. Default 0.0 (disabled — SCF behaves " "exactly as before).") .def_readwrite("newton_opts", &vibeqc::RHFOptions::newton_opts, "Knobs for the Newton step (CG cap, tolerance, " "trust radius). See NewtonOptions.") .def_readwrite("soscf_threshold", &vibeqc::RHFOptions::soscf_threshold, "Phase D2d — Neese SOSCF (approximate second-" "order SCF: augmented-Hessian step with diagonal-" "dominant Hessian) activation threshold on " "‖F D S − S D F‖_F. Cheaper per step than Newton " "(no Fock build), linearly convergent. Mutually " "exclusive with Newton; if both thresholds are " "set, Newton wins. Default 0.0 (disabled).") .def_readwrite("soscf_opts", &vibeqc::RHFOptions::soscf_opts, "Knobs for the SOSCF step (trust radius). See " "SOSCFOptions.") .def_readwrite("trah_threshold", &vibeqc::RHFOptions::trah_threshold, "Phase D2e TRAH (Trust-Region Augmented Hessian; " "Helmich-Paris 2022) activation threshold on " "‖F D S − S D F‖_F. Full Hessian like Newton, " "with an adaptive trust radius driven by Powell's " "rho test. Mutual exclusion with Newton + SOSCF " "(priority: quadratic > Newton > TRAH > SOSCF). " "Default 0.0 (disabled).") .def_readwrite("trah_opts", &vibeqc::RHFOptions::trah_opts, "Knobs for the TRAH step (CG cap + trust-region " "schedule). See TRAHOptions.") .def_readwrite("density_fit", &vibeqc::RHFOptions::density_fit, "Enable density fitting (RI) for the J + K " "Fock build. Replaces the four-index ERI tensor " "with the Cholesky-factorised B-tensor " "(Whitten 1973, Eichkorn et al. 1995). Reduces " "the per-iter Fock cost from O(n^4) to " "O(n^2 · n_aux). Requires aux_basis to be set. " "Default False.") .def_readwrite("aux_basis", &vibeqc::RHFOptions::aux_basis, "Auxiliary basis name for density fitting. " "libint-recognised name, e.g. \"def2-svp-jk\", " "\"cc-pvtz-jkfit\". Use " "vibeqc.default_aux_basis_for(orbital_basis_name, " "kind=\"jk\") for autodetection. Empty + " "density_fit=True raises.") .def_readwrite("scf_mode", &vibeqc::RHFOptions::scf_mode, "SCF Fock-build mode (SCFMode). AUTO = pick by " "basis size (default). CONVENTIONAL = in-core " "4-index ERI. DIRECT = Schwarz-screened on-the-" "fly libint quartet evaluation. Orthogonal to " "density_fit + cosx.") .def_readwrite("scf_mode_auto_threshold", &vibeqc::RHFOptions::scf_mode_auto_threshold, "AUTO cutoff: pick DIRECT when n_bf > " "threshold, else CONVENTIONAL. Default 200.") .def_readwrite("schwarz_threshold", &vibeqc::RHFOptions::schwarz_threshold, "Per-quartet Schwarz skip bound for the DIRECT " "kernel (|⟨μν|λσ⟩| ≤ Q_μν · Q_λσ). Default " "1e-10 (ORCA convention). 0 disables screening.") .def_readwrite("incremental_fock", &vibeqc::RHFOptions::incremental_fock, "Enable Almlöf-style incremental ΔP Fock build " "for the DIRECT path. Caches D_prev + G_2e_prev; " "per iter, builds G_2e[ΔD] and adds to prev — " "3-10× total-SCF speedup at large basis. " "Default False.") .def_readwrite("incremental_fock_reset_freq", &vibeqc::RHFOptions::incremental_fock_reset_freq, "Full-rebuild frequency for incremental_fock — " "every N iterations the ΔP cache is discarded " "and rebuilt from scratch to bound floating-" "point drift. Default 8 (ORCA's DirectResetFreq).") .def_readwrite("schwarz_threshold_loose", &vibeqc::RHFOptions::schwarz_threshold_loose, "Two-phase Schwarz threshold for DIRECT mode: " "start at this looser cutoff (more aggressive " "screening), tighten to ``schwarz_threshold`` " "once SCF grad-norm < schwarz_threshold_tighten_at. " "Default 1e-7 (ORCA convention); set ≤ " "schwarz_threshold to disable.") .def_readwrite("schwarz_threshold_tighten_at", &vibeqc::RHFOptions::schwarz_threshold_tighten_at, "Gradient-norm cutoff that triggers the loose→" "tight Schwarz transition. Default 1e-3.") .def_readwrite("cosx", &vibeqc::RHFOptions::cosx, "Use chain-of-spheres exchange (COSX) for the K " "build. Pair with density_fit=True for the " "standard RIJCOSX hybrid-DFT acceleration. " "Default False.") .def_readwrite("cosx_grid", &vibeqc::RHFOptions::cosx_grid, "GridOptions for the COSX integration grid (only " "used when cosx=True). Default = same as the DFT " "XC defaults; tune for the speed-accuracy " "trade-off.") .def_readwrite("dft_plus_u_sites", &vibeqc::RHFOptions::dft_plus_u_sites, "Internal: list of _HubbardSiteCxx (Hartree). " "Populated by the vibeqc.run_rhf Python wrapper " "from the user-facing dft_plus_u=[HubbardSite(...)] " "kwarg. Empty by default — empty ≡ no +U.") .def_readwrite("dft_plus_u_ao_groups", &vibeqc::RHFOptions::dft_plus_u_ao_groups, "Internal: parallel array of AO-index lists, one " "per Hubbard site. Precomputed once at SCF setup " "(see python/vibeqc/dft_plus_u.py::ao_group_indices) " "so the SCF loop pays the shell walk zero times."); py::class_(m, "NewtonOptions", "Knobs for the Newton Newton step (Phase D2c). See " "cpp/include/vibeqc/newton.hpp for the math + references " "(Bacskay 1981, Fischer/Almlöf 1992, Neese 2000, " "Helmich-Paris 2022).") .def(py::init<>()) .def_readwrite("cg_max_iter", &vibeqc::NewtonOptions::cg_max_iter, "CG iteration cap. 50 is generous; Newton normally " "converges in 5-15 iters once active. Default 50.") .def_readwrite("cg_tol", &vibeqc::NewtonOptions::cg_tol, "Relative residual tolerance on the CG solver: " "‖A κ + g‖ ≤ tol · ‖g‖. Default 1e-4.") .def_readwrite("trust_radius", &vibeqc::NewtonOptions::trust_radius, "Trust-region cap on ‖κ‖_F. Default 0.3."); py::class_(m, "SOSCFOptions", "Knobs for the Neese SOSCF step (Phase D2d) — augmented-" "Hessian eigsolve with diagonal-dominant orbital Hessian. " "See cpp/include/vibeqc/soscf.hpp for the math + reference " "(Neese 2000, Chem. Phys. Lett. 325, 93). Distinct from " "NewtonOptions (D2c, full Hessian via CG).") .def(py::init<>()) .def_readwrite("trust_radius", &vibeqc::SOSCFOptions::trust_radius, "Trust-region cap on ‖κ‖_F. Default 0.3 — same " "default as NewtonOptions for consistency."); py::class_(m, "TRAHOptions", "Knobs for the TRAH (Trust-Region Augmented Hessian) step " "(Phase D2e — Helmich-Paris 2022). Same full-Hessian CG as " "Newton; differs in the adaptive trust radius driven by " "Powell's rho test (actual / model-predicted dE). See " "cpp/include/vibeqc/trah.hpp.") .def(py::init<>()) .def_readwrite("cg_max_iter", &vibeqc::TRAHOptions::cg_max_iter) .def_readwrite("cg_tol", &vibeqc::TRAHOptions::cg_tol) .def_readwrite("initial_trust_radius", &vibeqc::TRAHOptions::initial_trust_radius) .def_readwrite("max_trust_radius", &vibeqc::TRAHOptions::max_trust_radius) .def_readwrite("min_trust_radius", &vibeqc::TRAHOptions::min_trust_radius) .def_readwrite("rho_shrink", &vibeqc::TRAHOptions::rho_shrink) .def_readwrite("rho_expand", &vibeqc::TRAHOptions::rho_expand) .def_readwrite("trust_shrink_factor", &vibeqc::TRAHOptions::trust_shrink_factor) .def_readwrite("trust_expand_factor", &vibeqc::TRAHOptions::trust_expand_factor) .def_readwrite("ms_max_iter", &vibeqc::TRAHOptions::ms_max_iter, "More-Sorensen secular-root-find iteration cap for " "the boundary level shift lambda.") .def_readwrite("ms_tol", &vibeqc::TRAHOptions::ms_tol, "Relative tolerance on ||kappa(lambda)|| = Delta for " "the More-Sorensen secular root-find."); py::class_(m, "SCFIteration") .def(py::init<>()) .def(py::init([](int iter, double energy, double delta_e, double grad_norm, int diis_subspace, int newton_cg_iter, double trah_level_shift) { vibeqc::SCFIteration it; it.iter = iter; it.energy = energy; it.delta_e = delta_e; it.grad_norm = grad_norm; it.diis_subspace = diis_subspace; it.newton_cg_iter = newton_cg_iter; it.trah_level_shift = trah_level_shift; return it; }), py::arg("iter") = 0, py::arg("energy") = 0.0, py::arg("delta_e") = 0.0, py::arg("grad_norm") = 0.0, py::arg("diis_subspace") = 0, py::arg("newton_cg_iter") = 0, py::arg("trah_level_shift") = 0.0, "Construct an SCFIteration record. Used by Python SCF drivers " "(periodic Ewald RHF / UHF / UKS) so every backend's " "``scf_trace`` carries the same item type — instead of plain " "tuples on the Ewald path and SCFIteration on the C++ path.") .def_readonly("iter", &vibeqc::SCFIteration::iter) .def_readonly("energy", &vibeqc::SCFIteration::energy, "Total energy at this iteration (Hartree).") .def_readonly("delta_e", &vibeqc::SCFIteration::delta_e, "E[k] - E[k-1]; 0 on the first iteration.") .def_readonly("grad_norm", &vibeqc::SCFIteration::grad_norm, "||F D S - S D F||_Frobenius; zero at convergence.") .def_readonly("diis_subspace", &vibeqc::SCFIteration::diis_subspace, "DIIS history size used this iteration (0 if off).") .def_readonly("newton_cg_iter", &vibeqc::SCFIteration::newton_cg_iter, "Newton CG iterations consumed this step (0 unless " "the Newton Newton update was active this iter).") .def_readonly("trah_level_shift", &vibeqc::SCFIteration::trah_level_shift, "TRAH level shift λ applied this step. 0 unless the " "TRAH update took a trust-region boundary step; " "> 0 ⇒ the step landed on the ‖κ‖ = Δ boundary " "with κ(λ) = −(A + λI)⁻¹g.") .def("__repr__", [](const vibeqc::SCFIteration& s) { return std::string{"SCFIteration(iter="} + std::to_string(s.iter) + ", energy=" + std::to_string(s.energy) + ", dE=" + std::to_string(s.delta_e) + ", grad=" + std::to_string(s.grad_norm) + ")"; }); py::class_(m, "RHFResult") .def_readonly("energy", &vibeqc::RHFResult::energy, "Total HF energy (Hartree).") .def_readonly("e_electronic", &vibeqc::RHFResult::e_electronic, "Electronic energy = total - nuclear_repulsion.") .def_readonly("e_dft_plus_u", &vibeqc::RHFResult::e_dft_plus_u, "Dudarev DFT+U contribution to the total energy " "(Hartree). 0 when no Hubbard sites are configured " "in RHFOptions.dft_plus_u_sites.") .def_readonly("n_iter", &vibeqc::RHFResult::n_iter) .def_readonly("converged", &vibeqc::RHFResult::converged) .def_readonly("mo_energies", &vibeqc::RHFResult::mo_energies, "Orbital energies (Hartree), ascending.") .def_readonly("mo_coeffs", &vibeqc::RHFResult::mo_coeffs, "MO coefficient matrix: columns are MOs in the AO basis.") .def_readonly("density", &vibeqc::RHFResult::density, "Converged density matrix D = 2 C_occ C_occ^T.") .def_readonly("fock", &vibeqc::RHFResult::fock, "Converged Fock matrix F = Hcore + G(D).") .def_readonly("scf_trace", &vibeqc::RHFResult::scf_trace, "Per-iteration SCF trace (list of SCFIteration).") .def("__repr__", [](const vibeqc::RHFResult& r) { return std::string{"RHFResult(energy="} + std::to_string(r.energy) + ", n_iter=" + std::to_string(r.n_iter) + ", converged=" + (r.converged ? "True" : "False") + ")"; }); m.def("run_rhf", &vibeqc::run_rhf, py::arg("molecule"), py::arg("basis"), py::arg("options") = vibeqc::RHFOptions{}, py::call_guard(), "Run restricted Hartree-Fock SCF on a closed-shell molecule. " "Returns an RHFResult."); m.def("run_rhf_scf_with_jk", [](const vibeqc::BasisSet& basis, int n_electrons, const Eigen::MatrixXd& S, const Eigen::MatrixXd& Hcore, double E_nuc, const std::shared_ptr& jk, const vibeqc::RHFOptions& options, const Eigen::MatrixXd& initial_density) { if (!jk) throw std::invalid_argument( "run_rhf_scf_with_jk: jk_builder must not be None"); return vibeqc::run_rhf_scf_with_jk( basis, n_electrons, S, Hcore, E_nuc, *jk, options, initial_density); }, py::arg("basis"), py::arg("n_electrons"), py::arg("S"), py::arg("Hcore"), py::arg("E_nuc"), py::arg("jk_builder"), py::arg("options") = vibeqc::RHFOptions{}, py::arg("initial_density") = Eigen::MatrixXd{}, py::call_guard(), "Lower-level closed-shell RHF SCF entry point. Takes the AO " "overlap S, the one-electron Hamiltonian Hcore, the system's " "nuclear repulsion E_nuc, and a JKBuilder for the two-electron " "Fock piece. Drives the same SCF iteration loop as run_rhf " "(canonical orthogonalisation, DIIS, optional damping / level-" "shift / quadratic fallback). Use it to drive periodic-Γ " "RHF or any SCF where Hcore is built externally — pair with " "make_periodic_gamma_jk_builder for periodic-Γ inputs. " "Pass ``initial_density`` (closed-shell density matrix) to " "seed D directly — empty (default) falls back to " "diagonalising Hcore."); // ----- Gradient ------------------------------------------------------- py::class_(m, "GradientOptions") .def(py::init<>()) .def_readwrite("density_fit", &vibeqc::GradientOptions::density_fit, "Enable density fitting for the two-electron " "gradient. Replaces the four-index ERI-derivative " "path with the J + α_HF·K factorisation through " "the precomputed B-tensor (Weigend 2002). Hcore, " "overlap, and (for DFT) XC pieces are unchanged. " "Default False.") .def_readwrite("aux_basis", &vibeqc::GradientOptions::aux_basis, "Auxiliary basis name for density fitting. " "See RHFOptions.aux_basis.") .def_readwrite("cosx", &vibeqc::GradientOptions::cosx, "Use chain-of-spheres exchange (COSX) for the K " "piece of the two-electron gradient. Pair with " "density_fit=True for the standard RIJCOSX path. " "Default False.") .def_readwrite("cosx_grid", &vibeqc::GradientOptions::cosx_grid, "GridOptions for the COSX integration grid (only " "used when cosx=True). Sparser default than the " "DFT XC grid; see " "vibeqc.cosx.default_cosx_grid_options."); m.def("compute_gradient", &vibeqc::compute_gradient, py::arg("molecule"), py::arg("basis"), py::arg("result"), py::arg("options") = vibeqc::GradientOptions{}, py::call_guard(), "Analytic RHF nuclear gradient. Returns an (n_atoms, 3) Eigen " "matrix in Hartree/bohr. `result` must come from a converged " "run_rhf(...) call. Pass options.density_fit=True with " "options.aux_basis to use the DF analytic gradient (Weigend " "2002, PCCP 4, 4285)."); m.def("compute_gradient_uhf", &vibeqc::compute_gradient_uhf, py::arg("molecule"), py::arg("basis"), py::arg("result"), py::arg("options") = vibeqc::GradientOptions{}, py::call_guard(), "Analytic UHF nuclear gradient. Returns an (n_atoms, 3) matrix " "in Hartree/bohr. `result` must come from a converged " "run_uhf(...) call. Pass options.density_fit=True with " "options.aux_basis to use the DF analytic gradient (composes " "J on D_α+D_β with per-spin K at α_HF/2). Required to opt " "into the v0.7.x f-shell direct-gradient bug workaround for " "open-shell systems."); // compute_gradient_rks has a GridOptions default argument; see below, // registered after GridOptions. // ----- Numerical integration grid (DFT prerequisite) ------------------ py::enum_(m, "AngularScheme", "Angular-grid family. ``Lebedev`` (Lebedev-Laikov) is the modern " "default for ~half the points at the same XC accuracy; " "``ProductGaussLegendre`` is the legacy n_θ × n_φ grid kept for " "parity-matrix bit-reproducibility.") .value("Lebedev", vibeqc::AngularScheme::Lebedev) .value("ProductGaussLegendre", vibeqc::AngularScheme::ProductGaussLegendre); py::enum_(m, "AtomicPartition", "Atomic-partition cell-function family. ``Becke`` is the " "iterated-polynomial 1988 scheme (J. Chem. Phys. 88, 2547); " "``Stratmann`` is the 1996 piecewise polynomial with |μ| ≥ a " "hard cutoff (Chem. Phys. Lett. 257, 213) — modern default in " "PySCF / NWChem / Q-Chem, asymptotically linear per grid point.") .value("Becke", vibeqc::AtomicPartition::Becke) .value("Stratmann", vibeqc::AtomicPartition::Stratmann); py::enum_(m, "AngularPruning", "Angular-order pruning scheme for Lebedev grids. ``None`` (the " "default) uses ``lebedev_order`` at every radial shell. " "``NWChem`` is the 3-tier Murray-Handy-Laming 1993 / SG-1 " "(Gill-Johnson-Pople 1993) scheme — inner core and outer tail " "use a smaller angular order, the chemically-active middle " "band uses the full order. Typical 25-40 % point-count " "reduction with sub-µHa accuracy cost on first-row organics.") .value("None", vibeqc::AngularPruning::None) .value("NWChem", vibeqc::AngularPruning::NWChem); py::class_(m, "GridOptions") .def(py::init<>()) .def_readwrite("n_radial", &vibeqc::GridOptions::n_radial) // Angular grid selector + per-scheme parameters. .def_property("angular", [](const vibeqc::GridOptions& o) { return o.angular; }, [](vibeqc::GridOptions& o, const py::object& v) { // Accept either the enum directly or a case-insensitive // string ("lebedev" / "product" / "product_gauss_legendre") // — the string form is what most user code wants. if (py::isinstance(v)) { o.angular = v.cast(); return; } std::string s = v.cast(); std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c){ return std::tolower(c); }); if (s == "lebedev") { o.angular = vibeqc::AngularScheme::Lebedev; } else if (s == "product" || s == "product_gauss_legendre" || s == "productgausslegendre") { o.angular = vibeqc::AngularScheme::ProductGaussLegendre; } else { throw std::invalid_argument( "GridOptions.angular: unknown value '" + s + "' (use 'lebedev' or 'product')"); } }, "Angular-grid family — either an AngularScheme enum or a " "case-insensitive string ('lebedev' / 'product').") .def_readwrite("lebedev_order", &vibeqc::GridOptions::lebedev_order, "Lebedev algebraic order. Bundled tiers: " "11/17/23/29/35/41/47/53 (50/110/194/302/434/590/" "770/974 points). Used only when angular=Lebedev.") .def_readwrite("n_theta", &vibeqc::GridOptions::n_theta, "Gauss-Legendre θ nodes per shell — used only " "with angular=ProductGaussLegendre.") .def_readwrite("n_phi", &vibeqc::GridOptions::n_phi, "Uniform φ nodes per shell — used only with " "angular=ProductGaussLegendre.") // Atomic-partition selector + per-scheme parameters. .def_property("partition", [](const vibeqc::GridOptions& o) { return o.partition; }, [](vibeqc::GridOptions& o, const py::object& v) { if (py::isinstance(v)) { o.partition = v.cast(); return; } std::string s = v.cast(); std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c){ return std::tolower(c); }); if (s == "becke") { o.partition = vibeqc::AtomicPartition::Becke; } else if (s == "stratmann") { o.partition = vibeqc::AtomicPartition::Stratmann; } else { throw std::invalid_argument( "GridOptions.partition: unknown value '" + s + "' (use 'becke' or 'stratmann')"); } }, "Atomic-partition cell-function — either an AtomicPartition " "enum or a case-insensitive string ('becke' / 'stratmann').") .def_readwrite("becke_k", &vibeqc::GridOptions::becke_k, "Iterated-polynomial smoothing order — used only " "with partition=Becke.") .def_property("angular_pruning", [](const vibeqc::GridOptions& o) { return o.angular_pruning; }, [](vibeqc::GridOptions& o, const py::object& v) { if (py::isinstance(v)) { o.angular_pruning = v.cast(); return; } std::string s = v.cast(); std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c){ return std::tolower(c); }); if (s == "none") { o.angular_pruning = vibeqc::AngularPruning::None; } else if (s == "nwchem") { o.angular_pruning = vibeqc::AngularPruning::NWChem; } else { throw std::invalid_argument( "GridOptions.angular_pruning: unknown value '" + s + "' (use 'none' or 'nwchem')"); } }, "Angular-order pruning scheme for Lebedev grids — either " "an AngularPruning enum or a case-insensitive string " "('none' / 'nwchem'). Only effective when angular=Lebedev."); py::class_(m, "Grid") .def_readonly("points", &vibeqc::Grid::points, "(N, 3) grid points in bohr.") .def_readonly("weights", &vibeqc::Grid::weights, "(N,) integration weights (bohr^3).") .def_readonly("atom_of_point", &vibeqc::Grid::atom_of_point, "(N,) index of the atom that 'owns' each point."); m.def("build_grid", &vibeqc::build_grid, py::arg("molecule"), py::arg("options") = vibeqc::GridOptions{}, py::call_guard(), "Build the DFT numerical-integration grid for the given molecule."); m.def("build_grid_periodic", &vibeqc::build_grid_periodic, py::arg("grid_molecule"), py::arg("partition_atom_positions"), py::arg("options") = vibeqc::GridOptions{}, py::call_guard(), "Build a DFT integration grid with the periodic-Becke partition. " "Points are generated around ``grid_molecule`` atoms (home cell), " "but the Becke fuzzy-cell partition denominator includes every " "atom in ``partition_atom_positions`` (home + images). The first " "len(grid_molecule.atoms()) entries of partition_atom_positions " "must equal the home-cell atom positions in order. Reduces to " "build_grid when partition_atom_positions == grid_molecule."); // Depends on GridOptions being registered (default arg), so it lands // here rather than beside the other gradient bindings. m.def("compute_gradient_rks", &vibeqc::compute_gradient_rks, py::arg("molecule"), py::arg("basis"), py::arg("result"), py::arg("grid_options") = vibeqc::GridOptions{}, py::arg("options") = vibeqc::GradientOptions{}, py::call_guard(), "Analytic RKS (closed-shell DFT) nuclear gradient. Supports " "LDA, GGA, and hybrid functionals. Pass options.density_fit" "=True with options.aux_basis for the DF analytic gradient. " "Returns (n_atoms, 3) in Hartree/bohr."); m.def("compute_gradient_uks", &vibeqc::compute_gradient_uks, py::arg("molecule"), py::arg("basis"), py::arg("result"), py::arg("grid_options") = vibeqc::GridOptions{}, py::arg("options") = vibeqc::GradientOptions{}, py::call_guard(), "Analytic UKS (open-shell DFT) nuclear gradient. Supports LDA, " "GGA, and hybrid functionals. Returns (n_atoms, 3) in Hartree/" "bohr. Pass options.density_fit=True with options.aux_basis " "for the DF analytic gradient (Fix D workaround for the " "v0.7.x f-shell direct-gradient bug; matches " "compute_gradient_uhf's per-spin J/K composition + XC Pulay " "piece unchanged)."); // ---- AO Fock-builder primitives ---- // Exposed so the analytic Hessian (Phase 17b-3) can build F at a // displaced geometry with a frozen reference density (no SCF), and // so external code can probe Coulomb / exchange contributions. m.def("build_fock_g", [](const py::array_t& eri, const Eigen::MatrixXd& D) { if (eri.ndim() != 4) { throw std::invalid_argument( "build_fock_g: eri must be a 4-D ndarray"); } vibeqc::Eri4D eri_obj; eri_obj.n = static_cast(eri.shape(0)); eri_obj.data.assign(eri.data(), eri.data() + eri.size()); return vibeqc::build_fock_g(eri_obj, D); }, py::arg("eri"), py::arg("D"), py::call_guard(), "Closed-shell Fock 2-electron contribution G[D] = J[D] − ½ K[D]. " "``eri`` is the (n, n, n, n) AO integral tensor from compute_eri."); m.def("build_coulomb", [](const py::array_t& eri, const Eigen::MatrixXd& D) { vibeqc::Eri4D eri_obj; eri_obj.n = static_cast(eri.shape(0)); eri_obj.data.assign(eri.data(), eri.data() + eri.size()); return vibeqc::build_coulomb(eri_obj, D); }, py::arg("eri"), py::arg("D"), py::call_guard(), "Coulomb matrix J[D]_μν = Σ_λσ (μν|λσ) D_λσ."); m.def("build_exchange", [](const py::array_t& eri, const Eigen::MatrixXd& D) { vibeqc::Eri4D eri_obj; eri_obj.n = static_cast(eri.shape(0)); eri_obj.data.assign(eri.data(), eri.data() + eri.size()); return vibeqc::build_exchange(eri_obj, D); }, py::arg("eri"), py::arg("D"), py::call_guard(), "Exchange matrix K[D]_μν = Σ_λσ (μλ|νσ) D_λσ."); // ---- Per-contribution gradient bindings ---- // These are the building blocks the high-level `compute_gradient_*` // wrap. Exposed so the analytic-Hessian assembly (Phase 17b-3) can // FD them with frozen weight matrices, and so cross-check tests can // verify each piece against PySCF independently. m.def("nuclear_repulsion_gradient", &vibeqc::nuclear_repulsion_gradient, py::arg("molecule"), py::call_guard(), "Closed-form ∂E_nuc/∂R per atom. Returns (n_atoms, 3) in Hartree/bohr."); m.def("overlap_gradient_contribution", &vibeqc::overlap_gradient_contribution, py::arg("basis"), py::arg("molecule"), py::arg("W"), py::call_guard(), "−Σ_μν W ∂S/∂R contraction (per atom). Pass the energy-" "weighted density W. Returns (n_atoms, 3) in Hartree/bohr."); m.def("one_electron_gradient_contribution", &vibeqc::one_electron_gradient_contribution, py::arg("basis"), py::arg("molecule"), py::arg("D"), py::call_guard(), "Σ_μν D ∂(T+V)/∂R contraction (per atom). Pass the closed-" "shell total density D. Returns (n_atoms, 3) in Hartree/bohr."); m.def("two_electron_gradient_contribution", &vibeqc::two_electron_gradient_contribution, py::arg("basis"), py::arg("molecule"), py::arg("D"), py::arg("alpha_hf") = 1.0, py::call_guard(), "Σ_μνλσ Γ ∂(μν|λσ)/∂R contraction (per atom). For plain HF " "use alpha_hf = 1; for hybrid DFT pass the functional's HF-" "exchange fraction. Returns (n_atoms, 3) in Hartree/bohr."); m.def("two_electron_gradient_contribution_uhf", &vibeqc::two_electron_gradient_contribution_uhf, py::arg("basis"), py::arg("molecule"), py::arg("D_alpha"), py::arg("D_beta"), py::arg("alpha_hf") = 1.0, py::call_guard(), "UHF/UKS counterpart of two_electron_gradient_contribution."); // CPCM solvation-gradient helper (v0.9.1). Analytic per-center // derivatives of the electron-density / external-point-charge // interaction — see vibeqc/gradient.hpp for the math. py::class_( m, "ExternalChargeGradient", "Per-center derivatives of E_ext = Σ_i q_i Tr(D·M_i) where " "M_i,μν = ⟨μ|1/|r−s_i||ν⟩. ``atom_grad`` (n_atoms, 3) carries " "the bra+ket basis-center derivatives; ``point_grad`` " "(n_points, 3) carries the per-external-charge derivatives. " "Both in Hartree/bohr with the sign of dE_ext/dR.") .def_readonly("atom_grad", &vibeqc::ExternalChargeGradient::atom_grad) .def_readonly("point_grad", &vibeqc::ExternalChargeGradient::point_grad); m.def("compute_external_charge_density_gradient", &vibeqc::compute_external_charge_density_gradient, py::arg("basis"), py::arg("molecule"), py::arg("D"), py::arg("charges"), py::arg("positions"), py::call_guard(), "Analytic per-center gradient of the electron density's " "interaction with external point charges — the CPCM " "electronic-ESP-derivative kernel. ``charges`` and " "``positions`` are the apparent surface charges q_i and " "cavity points s_i; ``D`` is the SCF density. Returns an " "ExternalChargeGradient (atom_grad + point_grad). Drives " "libint's nuclear-attraction gradient engine in one pass — " "the closed-form replacement for the integral-FD path in " "vibeqc.solvation.cpcm_gradient."); // ---- Phase G1a — periodic gradient primitives ------------------- // Lattice-summed analogues of the molecular gradient primitives // above, exposed individually so callers (and tests) can FD each // piece in isolation against the molecular reference. m.def("nuclear_repulsion_gradient_per_cell", &vibeqc::nuclear_repulsion_gradient_per_cell, py::arg("system"), py::arg("opts"), py::call_guard(), "Per-cell nuclear-repulsion gradient. Returns (n_atoms, 3) " "in Hartree/bohr. Same lattice-sum convention as " "nuclear_repulsion_per_cell."); m.def("overlap_lattice_gradient_contribution", &vibeqc::overlap_lattice_gradient_contribution, py::arg("basis"), py::arg("system"), py::arg("W_set"), py::arg("opts"), py::call_guard(), "−Σ_g W(g) · ∂S(g)/∂R contraction (per atom). Pass the real-" "space energy-weighted density on the same cell list as the " "matching compute_overlap_lattice call. Returns (n_atoms, 3)."); m.def("kinetic_lattice_gradient_contribution", &vibeqc::kinetic_lattice_gradient_contribution, py::arg("basis"), py::arg("system"), py::arg("D_set"), py::arg("opts"), py::call_guard(), "Σ_g D(g) · ∂T(g)/∂R contraction (per atom). Pass the real-" "space density on the same cell list as compute_kinetic_lattice. " "Returns (n_atoms, 3)."); m.def("nuclear_lattice_gradient_contribution", &vibeqc::nuclear_lattice_gradient_contribution, py::arg("basis"), py::arg("system"), py::arg("D_set"), py::arg("opts"), py::call_guard(), "Σ_g D(g) · ∂V(g)/∂R contraction (per atom). Includes the " "Σ_h-over-nuclear-images derivative — basis-center derivatives " "(2) plus one buffer per nuclear point charge in the lattice " "sum. Returns (n_atoms, 3)."); m.def("eri_lattice_gradient_contribution", &vibeqc::eri_lattice_gradient_contribution, py::arg("basis"), py::arg("system"), py::arg("D_set"), py::arg("opts"), py::arg("alpha_hf") = 1.0, py::call_guard(), "Periodic 2-electron ERI gradient: Σ_{g_a,g_b,g_c} Σ_μνλσ " "Γ(g_a,g_b,g_c) · ∂(μ_0 ν_{g_a} | λ_{g_b} σ_{g_c})/∂R with " "Γ = (1/2) D(g_a)_μν D(g_c-g_b)_λσ - (α_HF/4) D(g_b)_μλ " "D(g_c-g_a)_νσ. For pure DFT pass alpha_hf = 0; for plain " "RHF pass 1; for hybrid DFT pass the functional's HF-exchange " "fraction. Returns (n_atoms, 3)."); // ---- Spheropole kernel (EXT EL-SPHEROPOLE) ------------------------- m.def("compute_spheropole_kernel_lattice", &vibeqc::compute_spheropole_kernel_lattice, py::arg("prim_basis"), py::arg("system"), py::arg("opts"), py::arg("orig_nbf"), py::arg("c_map"), py::call_guard(), "Compute the EXT EL-SPHEROPOLE kernel K(g) for every lattice " "cell. Uses a primitive-split basis (one primitive per shell). " "Applies the CJAT33 per-(l_a,l_b) kernel with d^std scaling " "and contraction-map back to the original AO basis. Returns " "a LatticeMatrixSet. The global TOTCHR/(2V·π^{3/2}) factor " "is NOT included — apply after tracing with density."); // ---- Phase 17b-2 — Hessian skeleton-derivative contractions ---- m.def("compute_overlap_hessian_contribution", &vibeqc::compute_overlap_hessian_contribution, py::arg("basis"), py::arg("molecule"), py::arg("W"), py::call_guard(), "Skeleton −Σ W ∂²S/∂R_α∂R_β contraction. Returns a (3N, 3N) " "symmetric matrix in Hartree/bohr². Pass the energy-weighted " "density W = 2 C_occ · diag(ε_occ) · C_occ^T (closed-shell " "convention). Used internally by the analytic RHF Hessian."); m.def("compute_kinetic_nuclear_hessian_contribution", &vibeqc::compute_kinetic_nuclear_hessian_contribution, py::arg("basis"), py::arg("molecule"), py::arg("D"), py::call_guard(), "Skeleton Σ D ∂²(T+V)/∂R_α∂R_β contraction. Returns a (3N, 3N) " "symmetric matrix in Hartree/bohr². Pass the closed-shell total " "density D = 2 C_occ · C_occ^T. Includes nuclear-attraction " "second derivatives w.r.t. both basis-function and nuclear " "centers."); m.def("compute_eri_hessian_contribution", &vibeqc::compute_eri_hessian_contribution, py::arg("basis"), py::arg("molecule"), py::arg("D"), py::arg("alpha_hf") = 1.0, py::call_guard(), "Skeleton Σ Γ ∂²(μν|λσ)/∂R_α∂R_β contraction with the closed-" "shell two-particle density Γ = (1/2) D D − (α_HF/4) D D' " "(exchange permutation). For plain HF use α_HF = 1; for hybrid " "DFT pass the functional's HF-exchange fraction."); m.def("compute_eri_hessian_contribution_uhf", &vibeqc::compute_eri_hessian_contribution_uhf, py::arg("basis"), py::arg("molecule"), py::arg("D_alpha"), py::arg("D_beta"), py::arg("alpha_hf") = 1.0, py::call_guard(), "UHF / UKS counterpart of compute_eri_hessian_contribution. " "Two-particle density: Γ = (1/2)(D_α+D_β)(D_α+D_β) − " "(α_HF/2)(D_α D_α + D_β D_β) (per-spin exchange)."); m.def("nuclear_repulsion_hessian", &vibeqc::nuclear_repulsion_hessian, py::arg("molecule"), py::call_guard(), "Closed-form ∂²E_nuc/∂R_α∂R_β for the nuclear-repulsion " "energy. Returns a (3N, 3N) symmetric matrix."); m.def("evaluate_ao", &vibeqc::evaluate_ao, py::arg("basis"), py::arg("points"), py::call_guard(), "Evaluate χ_μ(r_g) for every basis function and every grid point. " "Returns an (n_points, n_basis) matrix."); m.def("evaluate_bloch_ao", &vibeqc::evaluate_bloch_ao, py::arg("basis"), py::arg("points"), py::arg("k_cart"), py::arg("lattice_translations"), py::call_guard(), "Bloch-summed AO matrix χ_μ^k(r_g) = Σ_T exp(i k·T) χ_μ(r_g − T) " "on a flat (n_points, 3) array of Cartesian grid points (bohr). " "``lattice_translations`` is an (n_T, 3) array of Cartesian " "lattice shifts T (bohr) — usually obtained by stacking the " "``r_cart`` of cells from ``direct_lattice_cells(system, cutoff)``. " "Returns a complex (n_points, n_basis) matrix. Contracting against " "a column of the multi-k SCF coefficient matrix C(k) recovers the " "Bloch crystalline orbital ψ_{n,k}(r); use the Python wrapper " "``vibeqc.evaluate_bloch_orbital`` for that one-liner."); // ----- XC functional (libxc wrapper) ---------------------------------- py::enum_(m, "XCKind") .value("LDA", vibeqc::XCKind::LDA) .value("GGA", vibeqc::XCKind::GGA) .value("MGGA", vibeqc::XCKind::MGGA); py::class_(m, "Functional") .def(py::init(), py::arg("name"), py::arg("spin") = 1, "Construct an XC functional by name ('LDA', 'PBE', 'BLYP', " "'B3LYP') or by comma-separated libxc IDs.") .def_property_readonly("name", &vibeqc::Functional::name) .def_property_readonly("kind", &vibeqc::Functional::kind) .def_property_readonly("is_hybrid", &vibeqc::Functional::is_hybrid) .def_property_readonly("hf_exchange_fraction", &vibeqc::Functional::hf_exchange_fraction) .def_property_readonly("is_range_separated", &vibeqc::Functional::is_range_separated, "True for a range-separated (CAM / RSH) " "hybrid — ωB97X, ωB97X-D, CAM-B3LYP, " "HSE06, …. The exact-exchange admixture " "is position-dependent: EXX(r) = " "cam_alpha + cam_beta·erf(rsh_omega·r).") .def_property_readonly("rsh_omega", &vibeqc::Functional::rsh_omega, "Range-separation parameter ω (bohr⁻¹) of " "a CAM/RSH hybrid; 0 for non-RSH.") .def_property_readonly("cam_alpha", &vibeqc::Functional::cam_alpha, "Full-range (always-on) HF-exchange " "fraction. Equals hf_exchange_fraction " "for a global hybrid; the constant CAM " "term α for an RSH hybrid.") .def_property_readonly("cam_beta", &vibeqc::Functional::cam_beta, "Long-range-only HF-exchange fraction β " "(switched on by erf(ω·r)). Zero for a " "global hybrid; for ωB97X-family " "long-range-corrected functionals " "cam_alpha + cam_beta = 1.") .def_property_readonly("is_double_hybrid", &vibeqc::Functional::is_double_hybrid, "True if this is a double-hybrid functional " "(B2PLYP, DSD-PBEP86, …). Double hybrids " "add a scaled-MP2 correlation correction " "on top of the SCF hybrid-DFT step; see " "``mp2_c_os`` / ``mp2_c_ss`` for the " "scaling coefficients and " "``vibeqc.run_b2plyp`` for the orchestrator.") .def_property_readonly("mp2_c_os", &vibeqc::Functional::mp2_c_os, "Opposite-spin MP2-correction coefficient " "for a double-hybrid functional. Zero for " "non-double-hybrid functionals. See " "``MP2Options.c_os``.") .def_property_readonly("mp2_c_ss", &vibeqc::Functional::mp2_c_ss, "Same-spin MP2-correction coefficient for " "a double-hybrid functional. Zero for " "non-double-hybrid functionals. See " "``MP2Options.c_ss``.") .def("eval_unpolarised", [](const vibeqc::Functional& f, const Eigen::VectorXd& rho, const Eigen::VectorXd& sigma) { Eigen::VectorXd exc, v_rho, v_sigma; f.eval_unpolarised(rho, sigma, exc, v_rho, v_sigma); return py::make_tuple(std::move(exc), std::move(v_rho), std::move(v_sigma)); }, py::arg("rho"), py::arg("sigma"), "Evaluate XC energy density and potentials on a grid. Returns " "(exc, v_rho, v_sigma). For LDA, sigma may be empty and v_sigma " "will be empty too.") .def("eval_unpolarised_fxc", [](const vibeqc::Functional& f, const Eigen::VectorXd& rho, const Eigen::VectorXd& sigma) { Eigen::VectorXd v2rho2, v2rhosigma, v2sigma2; f.eval_unpolarised_fxc(rho, sigma, v2rho2, v2rhosigma, v2sigma2); return py::make_tuple(std::move(v2rho2), std::move(v2rhosigma), std::move(v2sigma2)); }, py::arg("rho"), py::arg("sigma"), "Phase 17d — XC kernel for analytic KS Hessian / CPKS. Returns " "(v2rho2, v2rhosigma, v2sigma2). For LDA, only v2rho2 = " "∂²f/∂ρ² is populated; v2rhosigma and v2sigma2 are empty.") .def("eval_polarised_lda_fxc", [](const vibeqc::Functional& f, const Eigen::VectorXd& rho_a, const Eigen::VectorXd& rho_b) { Eigen::VectorXd v_aa, v_ab, v_bb; f.eval_polarised_lda_fxc(rho_a, rho_b, v_aa, v_ab, v_bb); return py::make_tuple(std::move(v_aa), std::move(v_ab), std::move(v_bb)); }, py::arg("rho_a"), py::arg("rho_b"), "LDA-only polarized XC kernel for UKS analytic Hessian. " "Returns (v2rho2_αα, v2rho2_αβ, v2rho2_ββ).") .def("eval_polarised_gga_fxc", [](const vibeqc::Functional& f, const Eigen::VectorXd& rho_a, const Eigen::VectorXd& rho_b, const Eigen::VectorXd& sigma_aa, const Eigen::VectorXd& sigma_ab, const Eigen::VectorXd& sigma_bb) { vibeqc::PolarisedGGAFxc out; f.eval_polarised_gga_fxc( rho_a, rho_b, sigma_aa, sigma_ab, sigma_bb, out); py::dict d; d["v2rho2_aa"] = out.v2rho2_aa; d["v2rho2_ab"] = out.v2rho2_ab; d["v2rho2_bb"] = out.v2rho2_bb; d["v2rhosigma_a_aa"] = out.v2rhosigma_a_aa; d["v2rhosigma_a_ab"] = out.v2rhosigma_a_ab; d["v2rhosigma_a_bb"] = out.v2rhosigma_a_bb; d["v2rhosigma_b_aa"] = out.v2rhosigma_b_aa; d["v2rhosigma_b_ab"] = out.v2rhosigma_b_ab; d["v2rhosigma_b_bb"] = out.v2rhosigma_b_bb; d["v2sigma2_aa_aa"] = out.v2sigma2_aa_aa; d["v2sigma2_aa_ab"] = out.v2sigma2_aa_ab; d["v2sigma2_aa_bb"] = out.v2sigma2_aa_bb; d["v2sigma2_ab_ab"] = out.v2sigma2_ab_ab; d["v2sigma2_ab_bb"] = out.v2sigma2_ab_bb; d["v2sigma2_bb_bb"] = out.v2sigma2_bb_bb; return d; }, py::arg("rho_a"), py::arg("rho_b"), py::arg("sigma_aa"), py::arg("sigma_ab"), py::arg("sigma_bb"), "Polarized LDA/GGA XC kernel for UKS analytic Hessian. " "Returns a dict with the 15 libxc second-derivative " "blocks. Raises on meta-GGA functionals.") .def("eval_polarised", [](const vibeqc::Functional& f, const Eigen::VectorXd& rho_a, const Eigen::VectorXd& rho_b, const Eigen::VectorXd& sigma_aa, const Eigen::VectorXd& sigma_ab, const Eigen::VectorXd& sigma_bb) { Eigen::VectorXd exc, v_rho_a, v_rho_b; Eigen::VectorXd v_sigma_aa, v_sigma_ab, v_sigma_bb; f.eval_polarised(rho_a, rho_b, sigma_aa, sigma_ab, sigma_bb, exc, v_rho_a, v_rho_b, v_sigma_aa, v_sigma_ab, v_sigma_bb); return py::make_tuple(std::move(exc), std::move(v_rho_a), std::move(v_rho_b), std::move(v_sigma_aa), std::move(v_sigma_ab), std::move(v_sigma_bb)); }, py::arg("rho_a"), py::arg("rho_b"), py::arg("sigma_aa"), py::arg("sigma_ab"), py::arg("sigma_bb"), "Polarized XC energy density + first derivatives on a grid. " "Returns (exc, v_rho_a, v_rho_b, v_sigma_aa, v_sigma_ab, " "v_sigma_bb). For LDA, sigma_* may be empty and the v_sigma_* " "outputs will be empty too. The Functional must be constructed " "with spin=2.") .def("eval_unpolarised_mgga", [](const vibeqc::Functional& f, const Eigen::VectorXd& rho, const Eigen::VectorXd& sigma, const Eigen::VectorXd& tau) { Eigen::VectorXd exc, v_rho, v_sigma, v_tau; f.eval_unpolarised_mgga(rho, sigma, tau, exc, v_rho, v_sigma, v_tau); return py::make_tuple(std::move(exc), std::move(v_rho), std::move(v_sigma), std::move(v_tau)); }, py::arg("rho"), py::arg("sigma"), py::arg("tau"), "Meta-GGA evaluation (τ-dependent). Returns " "(exc, v_rho, v_sigma, v_tau). τ = ½ Σ_i |∇ψ_i(r)|² is the " "total (both spins) for the unpolarised path.") .def("eval_polarised_mgga", [](const vibeqc::Functional& f, const Eigen::VectorXd& rho_a, const Eigen::VectorXd& rho_b, const Eigen::VectorXd& sigma_aa, const Eigen::VectorXd& sigma_ab, const Eigen::VectorXd& sigma_bb, const Eigen::VectorXd& tau_a, const Eigen::VectorXd& tau_b) { Eigen::VectorXd exc, v_rho_a, v_rho_b; Eigen::VectorXd v_sigma_aa, v_sigma_ab, v_sigma_bb; Eigen::VectorXd v_tau_a, v_tau_b; f.eval_polarised_mgga(rho_a, rho_b, sigma_aa, sigma_ab, sigma_bb, tau_a, tau_b, exc, v_rho_a, v_rho_b, v_sigma_aa, v_sigma_ab, v_sigma_bb, v_tau_a, v_tau_b); return py::make_tuple(std::move(exc), std::move(v_rho_a), std::move(v_rho_b), std::move(v_sigma_aa), std::move(v_sigma_ab), std::move(v_sigma_bb), std::move(v_tau_a), std::move(v_tau_b)); }, py::arg("rho_a"), py::arg("rho_b"), py::arg("sigma_aa"), py::arg("sigma_ab"), py::arg("sigma_bb"), py::arg("tau_a"), py::arg("tau_b"), "Polarized meta-GGA evaluation. Returns " "(exc, v_rho_a, v_rho_b, v_sigma_aa, v_sigma_ab, v_sigma_bb, " "v_tau_a, v_tau_b). τ_σ = ½ Σ_iσ |∇ψ_iσ(r)|² is per-spin."); // Phase D2c-KS — XC kernel matvec (W^XC[D_pert]) for the second-order // KS SCF path. Closed-shell (RKS) only in this commit; UKS path // exists as a factory stub that throws with a roadmap pointer. See // cpp/include/vibeqc/xc_kernel.hpp for the design + math. py::class_(m, "XCKernelBuilder", "Closed-shell unpolarised XC kernel builder. Pin the reference " "density at construction (via make_unpolarised_xc_kernel_builder); " "each ``apply(D_pert)`` returns W^XC[D_pert] in AO basis.") .def("apply", &vibeqc::XCKernelBuilder::apply, py::arg("D_pert"), "AO-basis W^XC[D_pert] matvec. Symmetric output."); m.def("make_unpolarised_xc_kernel_builder", [](const vibeqc::Functional& func, const vibeqc::Grid& grid, const Eigen::MatrixXd& chi, const Eigen::MatrixXd& dchi_x, const Eigen::MatrixXd& dchi_y, const Eigen::MatrixXd& dchi_z, const Eigen::MatrixXd& D_used) { std::array dchi = {dchi_x, dchi_y, dchi_z}; return vibeqc::make_unpolarised_xc_kernel_builder( func, grid, chi, dchi, D_used); }, py::arg("functional"), py::arg("grid"), py::arg("chi"), py::arg("dchi_x"), py::arg("dchi_y"), py::arg("dchi_z"), py::arg("D_used"), "Phase D2c-KS — build an unpolarised LDA / GGA XC kernel for " "the Hessian matvec at reference density D_used. The grid + " "chi + dchi triple comes from build_grid + " "evaluate_ao_with_gradient. Returns a pinned-state builder; " "call .apply(D_pert) for each Newton-CG matvec."); // Open-shell (polarised LDA) XC kernel. py::class_(m, "UHFXCKernelBuilder", "Open-shell polarised XC kernel builder. Per-spin matvec " "couples α and β via the αβ v2rho2 cross-term.") .def("apply", [](const vibeqc::UHFXCKernelBuilder& k, const Eigen::MatrixXd& D_pert_alpha, const Eigen::MatrixXd& D_pert_beta) { auto out = k.apply(D_pert_alpha, D_pert_beta); return py::make_tuple(std::move(out.alpha), std::move(out.beta)); }, py::arg("D_pert_alpha"), py::arg("D_pert_beta"), "Returns (W^XC_α, W^XC_β), each (n_bf, n_bf), symmetric."); m.def("make_polarised_lda_xc_kernel_builder", &vibeqc::make_polarised_lda_xc_kernel_builder, py::arg("functional"), py::arg("grid"), py::arg("chi"), py::arg("D_used_alpha"), py::arg("D_used_beta"), "Phase D2c-KS-UHF — build an open-shell polarised LDA XC " "kernel for the Hessian matvec at reference per-spin " "densities. Raises on GGA functionals (polarised GGA fxc " "isn't plumbed in Functional yet — Phase 17e)."); m.def("make_polarised_gga_xc_kernel_builder", [](const vibeqc::Functional& func, const vibeqc::Grid& grid, const Eigen::MatrixXd& chi, const Eigen::MatrixXd& dchi_x, const Eigen::MatrixXd& dchi_y, const Eigen::MatrixXd& dchi_z, const Eigen::MatrixXd& D_used_alpha, const Eigen::MatrixXd& D_used_beta) { std::array dchi = {dchi_x, dchi_y, dchi_z}; return vibeqc::make_polarised_gga_xc_kernel_builder( func, grid, chi, dchi, D_used_alpha, D_used_beta); }, py::arg("functional"), py::arg("grid"), py::arg("chi"), py::arg("dchi_x"), py::arg("dchi_y"), py::arg("dchi_z"), py::arg("D_used_alpha"), py::arg("D_used_beta"), "Phase D2c-KS-UHF — build an open-shell polarised GGA XC " "kernel for the Hessian matvec at reference per-spin " "densities."); m.def("make_polarised_xc_kernel_builder", [](const vibeqc::Functional& func, const vibeqc::Grid& grid, const Eigen::MatrixXd& chi, const Eigen::MatrixXd& dchi_x, const Eigen::MatrixXd& dchi_y, const Eigen::MatrixXd& dchi_z, const Eigen::MatrixXd& D_used_alpha, const Eigen::MatrixXd& D_used_beta) { std::array dchi = {dchi_x, dchi_y, dchi_z}; return vibeqc::make_polarised_xc_kernel_builder( func, grid, chi, dchi, D_used_alpha, D_used_beta); }, py::arg("functional"), py::arg("grid"), py::arg("chi"), py::arg("dchi_x"), py::arg("dchi_y"), py::arg("dchi_z"), py::arg("D_used_alpha"), py::arg("D_used_beta"), "Build an open-shell polarised XC kernel builder, dispatching " "to LDA or GGA based on the functional kind."); m.def("evaluate_ao_with_gradient", [](const vibeqc::BasisSet& basis, const Eigen::MatrixX3d& points) { auto r = vibeqc::evaluate_ao_with_gradient(basis, points); // Return values and the three gradient matrices as a plain // Python tuple. Construct outside the GIL-released region so // pybind11 can convert Eigen -> NumPy safely. py::gil_scoped_acquire gil; return py::make_tuple(r.values, r.gradients[0], r.gradients[1], r.gradients[2]); }, py::arg("basis"), py::arg("points"), "Evaluate AO values AND first spatial derivatives. Returns a tuple " "(values, grad_x, grad_y, grad_z), each an (n_points, n_basis) matrix."); // ----- UHF ----------------------------------------------------------- py::class_(m, "UHFOptions") .def(py::init<>()) .def_readwrite("max_iter", &vibeqc::UHFOptions::max_iter) .def_readwrite("conv_tol_energy", &vibeqc::UHFOptions::conv_tol_energy) .def_readwrite("conv_tol_grad", &vibeqc::UHFOptions::conv_tol_grad) .def_readwrite("damping", &vibeqc::UHFOptions::damping) .def_readwrite("dynamic_damping", &vibeqc::UHFOptions::dynamic_damping, "See RHFOptions::dynamic_damping. Default False.") .def_readwrite("dynamic_damping_min", &vibeqc::UHFOptions::dynamic_damping_min) .def_readwrite("dynamic_damping_max", &vibeqc::UHFOptions::dynamic_damping_max) .def_readwrite("fock_mixing", &vibeqc::UHFOptions::fock_mixing, "Per-spin Fock matrix mixing fraction. See " "RHFOptions.fock_mixing.") .def_readwrite("use_diis", &vibeqc::UHFOptions::use_diis) .def_readwrite("diis_start_iter", &vibeqc::UHFOptions::diis_start_iter) .def_readwrite("diis_subspace_size", &vibeqc::UHFOptions::diis_subspace_size) .def_readwrite("scf_accelerator", &vibeqc::UHFOptions::scf_accelerator, "SCF Fock-extrapolation accelerator. See " "RHFOptions.scf_accelerator. The EDIIS / " "EDIIS_DIIS variants apply a single coefficient " "set to both spins (the energy functional couples " "them through E_total).") .def_readwrite("ediis_diis_switch_threshold", &vibeqc::UHFOptions::ediis_diis_switch_threshold, "See RHFOptions.ediis_diis_switch_threshold. " "Default 1e-1.") .def_readwrite("initial_guess", &vibeqc::UHFOptions::initial_guess) .def_readwrite("linear_dep_threshold", &vibeqc::UHFOptions::linear_dep_threshold) .def_readwrite("ecp_centers", &vibeqc::UHFOptions::ecp_centers, "List[ECPCenter]. See RHFOptions.ecp_centers.") .def_readwrite("ecp_library", &vibeqc::UHFOptions::ecp_library, "ECP XML library name. See RHFOptions.ecp_library.") .def_readwrite("level_shift", &vibeqc::UHFOptions::level_shift, "Phase C1a-2 — per-spin Saunders-Hillier level " "shift (Hartree). Formula: " "F_σ + b·S − b·(S·D_σ·S). See " "RHFOptions.level_shift. Default 0.0.") .def_readwrite("quadratic_fallback_iter", &vibeqc::UHFOptions::quadratic_fallback_iter, "Phase C1c — per-spin second-order SCF " "fallback. See RHFOptions.quadratic_fallback_iter; " "the Newton step runs independently on the α " "and β MO sets. Default 0 (disabled).") .def_readwrite("quadratic_fallback_shift", &vibeqc::UHFOptions::quadratic_fallback_shift, "C1c orbital-rotation denominator shift. " "See RHFOptions.quadratic_fallback_shift.") .def_readwrite("quadratic_fallback_max_step", &vibeqc::UHFOptions::quadratic_fallback_max_step, "C1c per-step ‖κ‖ trust-region cap. " "See RHFOptions.quadratic_fallback_max_step.") .def_readwrite("newton_threshold", &vibeqc::UHFOptions::newton_threshold, "Phase D2c — Newton activation threshold " "(max-over-spins commutator norm). See " "RHFOptions.newton_threshold. The open-shell " "Newton step couples α and β rotations through " "the shared J in the orbital Hessian. " "Default 0.0 (disabled).") .def_readwrite("newton_opts", &vibeqc::UHFOptions::newton_opts, "Knobs for the open-shell Newton step (CG cap, " "tolerance, trust radius). Same NewtonOptions " "type as RHF; the trust radius applies to the " "combined ‖(κ_α, κ_β)‖_F.") .def_readwrite("soscf_threshold", &vibeqc::UHFOptions::soscf_threshold, "Phase D2d Neese SOSCF activation threshold " "for open-shell. Diagonal-dominant Hessian → " "per-spin AH eigsolves. Mutually exclusive " "with Newton. Default 0.0 (disabled).") .def_readwrite("soscf_opts", &vibeqc::UHFOptions::soscf_opts, "Knobs for the open-shell SOSCF step (trust " "radius). Same SOSCFOptions type as RHF.") .def_readwrite("trah_threshold", &vibeqc::UHFOptions::trah_threshold, "Phase D2e TRAH activation threshold for open-" "shell. Same coupled per-spin CG as Newton + " "adaptive trust radius. Priority: quadratic > " "Newton > TRAH > SOSCF. Default 0.0 (disabled).") .def_readwrite("trah_opts", &vibeqc::UHFOptions::trah_opts, "Knobs for the open-shell TRAH step (CG cap + " "trust-region schedule). Same TRAHOptions type " "as RHF; trust radius applies to combined " "‖(κ_α, κ_β)‖_F.") .def_readwrite("density_fit", &vibeqc::UHFOptions::density_fit, "Enable density fitting (RI) for the per-spin " "Fock build: J(D_total) and K(D_α/D_β) via " "the precomputed B-tensor instead of the four-" "index ERI. See RHFOptions.density_fit. " "Default False.") .def_readwrite("aux_basis", &vibeqc::UHFOptions::aux_basis, "Auxiliary basis name for density fitting. " "See RHFOptions.aux_basis.") .def_readwrite("scf_mode", &vibeqc::UHFOptions::scf_mode, "SCF Fock-build mode (SCFMode). See " "RHFOptions.scf_mode.") .def_readwrite("scf_mode_auto_threshold", &vibeqc::UHFOptions::scf_mode_auto_threshold, "AUTO cutoff (n_bf). See " "RHFOptions.scf_mode_auto_threshold.") .def_readwrite("schwarz_threshold", &vibeqc::UHFOptions::schwarz_threshold, "Per-quartet Schwarz skip bound for DIRECT mode. " "See RHFOptions.schwarz_threshold.") .def_readwrite("incremental_fock", &vibeqc::UHFOptions::incremental_fock, "Enable incremental ΔP Fock build for DIRECT. " "See RHFOptions.incremental_fock.") .def_readwrite("incremental_fock_reset_freq", &vibeqc::UHFOptions::incremental_fock_reset_freq, "Full-rebuild frequency for incremental_fock. " "See RHFOptions.incremental_fock_reset_freq.") .def_readwrite("schwarz_threshold_loose", &vibeqc::UHFOptions::schwarz_threshold_loose, "Two-phase Schwarz loose threshold. " "See RHFOptions.schwarz_threshold_loose.") .def_readwrite("schwarz_threshold_tighten_at", &vibeqc::UHFOptions::schwarz_threshold_tighten_at, "Grad-norm cutoff for the loose→tight transition. " "See RHFOptions.schwarz_threshold_tighten_at.") .def_readwrite("cosx", &vibeqc::UHFOptions::cosx, "Use chain-of-spheres exchange for K. See " "RHFOptions.cosx. Default False.") .def_readwrite("cosx_grid", &vibeqc::UHFOptions::cosx_grid, "GridOptions for the COSX grid. See " "RHFOptions.cosx_grid.") .def_readwrite("dft_plus_u_sites", &vibeqc::UHFOptions::dft_plus_u_sites, "Internal: list of _HubbardSiteCxx. Populated by " "the vibeqc.run_uhf Python wrapper.") .def_readwrite("dft_plus_u_ao_groups", &vibeqc::UHFOptions::dft_plus_u_ao_groups, "Internal: parallel AO-index lists per site."); py::class_(m, "UHFResult") .def_readonly("energy", &vibeqc::UHFResult::energy, "Total HF energy (Hartree).") .def_readonly("e_electronic", &vibeqc::UHFResult::e_electronic) .def_readonly("e_dft_plus_u", &vibeqc::UHFResult::e_dft_plus_u, "Dudarev DFT+U contribution to the total energy " "(Hartree). 0 unless UHFOptions.dft_plus_u_sites " "is non-empty. Sum of per-spin contributions.") .def_readonly("n_iter", &vibeqc::UHFResult::n_iter) .def_readonly("converged", &vibeqc::UHFResult::converged) .def_readonly("s_squared", &vibeqc::UHFResult::s_squared, " expectation value (spin-contamination diagnostic).") .def_readonly("s_squared_ideal", &vibeqc::UHFResult::s_squared_ideal, "S(S+1) for the requested multiplicity.") .def_readonly("mo_energies_alpha", &vibeqc::UHFResult::mo_energies_alpha) .def_readonly("mo_coeffs_alpha", &vibeqc::UHFResult::mo_coeffs_alpha) .def_readonly("density_alpha", &vibeqc::UHFResult::density_alpha) .def_readonly("fock_alpha", &vibeqc::UHFResult::fock_alpha) .def_readonly("mo_energies_beta", &vibeqc::UHFResult::mo_energies_beta) .def_readonly("mo_coeffs_beta", &vibeqc::UHFResult::mo_coeffs_beta) .def_readonly("density_beta", &vibeqc::UHFResult::density_beta) .def_readonly("fock_beta", &vibeqc::UHFResult::fock_beta) .def_readonly("scf_trace", &vibeqc::UHFResult::scf_trace) .def("__repr__", [](const vibeqc::UHFResult& r) { return std::string{"UHFResult(energy="} + std::to_string(r.energy) + ", n_iter=" + std::to_string(r.n_iter) + ", converged=" + (r.converged ? "True" : "False") + ", =" + std::to_string(r.s_squared) + ")"; }); // ----- RKS (closed-shell DFT) ---------------------------------------- py::class_(m, "RKSOptions") .def(py::init<>()) .def_readwrite("functional", &vibeqc::RKSOptions::functional) .def_readwrite("grid", &vibeqc::RKSOptions::grid) .def_readwrite("max_iter", &vibeqc::RKSOptions::max_iter) .def_readwrite("conv_tol_energy", &vibeqc::RKSOptions::conv_tol_energy) .def_readwrite("conv_tol_grad", &vibeqc::RKSOptions::conv_tol_grad) .def_readwrite("damping", &vibeqc::RKSOptions::damping) .def_readwrite("dynamic_damping", &vibeqc::RKSOptions::dynamic_damping, "See RHFOptions::dynamic_damping. Default False.") .def_readwrite("dynamic_damping_min", &vibeqc::RKSOptions::dynamic_damping_min) .def_readwrite("dynamic_damping_max", &vibeqc::RKSOptions::dynamic_damping_max) .def_readwrite("fock_mixing", &vibeqc::RKSOptions::fock_mixing, "Kohn-Sham matrix mixing fraction. See " "RHFOptions.fock_mixing.") .def_readwrite("use_diis", &vibeqc::RKSOptions::use_diis) .def_readwrite("diis_start_iter", &vibeqc::RKSOptions::diis_start_iter) .def_readwrite("diis_subspace_size", &vibeqc::RKSOptions::diis_subspace_size) .def_readwrite("scf_accelerator", &vibeqc::RKSOptions::scf_accelerator, "SCF Fock-extrapolation accelerator. See " "RHFOptions.scf_accelerator.") .def_readwrite("ediis_diis_switch_threshold", &vibeqc::RKSOptions::ediis_diis_switch_threshold, "See RHFOptions.ediis_diis_switch_threshold. " "Default 1e-1.") .def_readwrite("initial_guess", &vibeqc::RKSOptions::initial_guess) .def_readwrite("linear_dep_threshold", &vibeqc::RKSOptions::linear_dep_threshold) .def_readwrite("ecp_centers", &vibeqc::RKSOptions::ecp_centers, "List[ECPCenter]. See RHFOptions.ecp_centers.") .def_readwrite("ecp_library", &vibeqc::RKSOptions::ecp_library, "ECP XML library name. See RHFOptions.ecp_library.") .def_readwrite("level_shift", &vibeqc::RKSOptions::level_shift, "Phase C1a-2 — Saunders-Hillier level shift " "(Hartree). Same formula as RHF (closed-shell " "D = 2·C_occ·C_occ^T). See " "RHFOptions.level_shift. Default 0.0.") .def_readwrite("quadratic_fallback_iter", &vibeqc::RKSOptions::quadratic_fallback_iter, "Phase C1c — second-order SCF fallback. " "See RHFOptions.quadratic_fallback_iter. " "Default 0 (disabled).") .def_readwrite("quadratic_fallback_shift", &vibeqc::RKSOptions::quadratic_fallback_shift, "C1c orbital-rotation denominator shift. " "See RHFOptions.quadratic_fallback_shift.") .def_readwrite("quadratic_fallback_max_step", &vibeqc::RKSOptions::quadratic_fallback_max_step, "C1c per-step ‖κ‖ trust-region cap. " "See RHFOptions.quadratic_fallback_max_step.") .def_readwrite("newton_threshold", &vibeqc::RKSOptions::newton_threshold, "Phase D2c-KS Newton activation threshold. " "See RHFOptions.newton_threshold. Default 0.0 " "(disabled). XC kernel built via " "make_unpolarised_xc_kernel_builder; raises on " "meta-GGA functionals (Phase 17e+ classification).") .def_readwrite("newton_opts", &vibeqc::RKSOptions::newton_opts, "NewtonOptions (CG knobs + trust-region cap). " "See NewtonOptions.") .def_readwrite("soscf_threshold", &vibeqc::RKSOptions::soscf_threshold, "Phase D2d-KS Neese SOSCF activation threshold. " "Default 0.0. No XC kernel needed (F already " "carries V_xc).") .def_readwrite("soscf_opts", &vibeqc::RKSOptions::soscf_opts, "SOSCFOptions. See SOSCFOptions.") .def_readwrite("trah_threshold", &vibeqc::RKSOptions::trah_threshold, "Phase D2e-KS TRAH activation threshold. Default " "0.0. Same XC kernel as Newton. Mutual-exclusion " "priority: quadratic > Newton > TRAH > SOSCF.") .def_readwrite("trah_opts", &vibeqc::RKSOptions::trah_opts, "TRAHOptions (CG + Powell-ρ trust schedule).") .def_readwrite("density_fit", &vibeqc::RKSOptions::density_fit, "Enable density fitting (RI) for the J build " "(and K build for hybrid functionals). See " "RHFOptions.density_fit. Default False.") .def_readwrite("aux_basis", &vibeqc::RKSOptions::aux_basis, "Auxiliary basis name for density fitting. " "See RHFOptions.aux_basis.") .def_readwrite("scf_mode", &vibeqc::RKSOptions::scf_mode, "SCF Fock-build mode (SCFMode). See " "RHFOptions.scf_mode.") .def_readwrite("scf_mode_auto_threshold", &vibeqc::RKSOptions::scf_mode_auto_threshold, "AUTO cutoff (n_bf). See " "RHFOptions.scf_mode_auto_threshold.") .def_readwrite("schwarz_threshold", &vibeqc::RKSOptions::schwarz_threshold, "Per-quartet Schwarz skip bound for DIRECT mode. " "See RHFOptions.schwarz_threshold.") .def_readwrite("incremental_fock", &vibeqc::RKSOptions::incremental_fock, "Enable incremental ΔP Fock build for DIRECT. " "See RHFOptions.incremental_fock.") .def_readwrite("incremental_fock_reset_freq", &vibeqc::RKSOptions::incremental_fock_reset_freq, "Full-rebuild frequency for incremental_fock. " "See RHFOptions.incremental_fock_reset_freq.") .def_readwrite("schwarz_threshold_loose", &vibeqc::RKSOptions::schwarz_threshold_loose, "Two-phase Schwarz loose threshold. " "See RHFOptions.schwarz_threshold_loose.") .def_readwrite("schwarz_threshold_tighten_at", &vibeqc::RKSOptions::schwarz_threshold_tighten_at, "Grad-norm cutoff for the loose→tight transition. " "See RHFOptions.schwarz_threshold_tighten_at.") .def_readwrite("cosx", &vibeqc::RKSOptions::cosx, "Use chain-of-spheres exchange (COSX) for the K " "build. Pair with density_fit=True for the " "standard RIJCOSX hybrid-DFT acceleration. " "No-op for pure DFT (α_HF = 0). Default False.") .def_readwrite("cosx_grid", &vibeqc::RKSOptions::cosx_grid, "GridOptions for the COSX integration grid (only " "used when cosx=True). Defaults to a sparser " "tier than the XC grid; see " "vibeqc.cosx.default_cosx_grid_options.") .def_readwrite("dft_plus_u_sites", &vibeqc::RKSOptions::dft_plus_u_sites, "Internal: list of _HubbardSiteCxx (Hartree). " "Populated by the vibeqc.run_rks Python wrapper " "from the user-facing dft_plus_u=[HubbardSite(...)] " "kwarg.") .def_readwrite("dft_plus_u_ao_groups", &vibeqc::RKSOptions::dft_plus_u_ao_groups, "Internal: parallel array of AO-index lists, one " "per Hubbard site. Precomputed by the Python " "wrapper from ao_group_indices(basis)."); py::class_(m, "RKSResult") .def_readonly("energy", &vibeqc::RKSResult::energy) .def_readonly("e_electronic", &vibeqc::RKSResult::e_electronic) .def_readonly("e_coulomb", &vibeqc::RKSResult::e_coulomb) .def_readonly("e_hf_exchange", &vibeqc::RKSResult::e_hf_exchange) .def_readonly("e_xc", &vibeqc::RKSResult::e_xc) .def_readonly("e_nuclear", &vibeqc::RKSResult::e_nuclear) .def_readonly("e_dft_plus_u", &vibeqc::RKSResult::e_dft_plus_u, "Dudarev DFT+U contribution to the total energy " "(Hartree). 0 when RKSOptions.dft_plus_u_sites is " "empty.") .def_readonly("n_iter", &vibeqc::RKSResult::n_iter) .def_readonly("converged", &vibeqc::RKSResult::converged) .def_readonly("mo_energies", &vibeqc::RKSResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::RKSResult::mo_coeffs) .def_readonly("density", &vibeqc::RKSResult::density) .def_readonly("fock", &vibeqc::RKSResult::fock) .def_readonly("scf_trace", &vibeqc::RKSResult::scf_trace) .def_readonly("functional", &vibeqc::RKSResult::functional) .def("__repr__", [](const vibeqc::RKSResult& r) { return std::string{"RKSResult(energy="} + std::to_string(r.energy) + ", functional=" + r.functional + ", converged=" + (r.converged ? "True" : "False") + ")"; }); m.def("run_rks", &vibeqc::run_rks, py::arg("molecule"), py::arg("basis"), py::arg("options") = vibeqc::RKSOptions{}, py::call_guard(), "Run restricted Kohn-Sham DFT on a closed-shell molecule."); m.def("rhf_result_from_rks", &vibeqc::rhf_result_from_rks, py::arg("rks_result"), "Adapter that returns an RHFResult populated from an " "RKSResult — the same orbitals, energy, density, and Fock " "matrix, just typed as an RHF reference so it can be passed " "to ``run_mp2``. The intended consumer is the double-hybrid " "dispatch (``vibeqc.run_b2plyp`` etc.), which runs RKS with " "the hybrid SCF piece of a double-hybrid functional and then " "feeds the KS orbitals into MP2 with the functional's " "``mp2_c_os`` / ``mp2_c_ss`` scaling."); m.def("run_rks_scf_with_jk", [](const vibeqc::BasisSet& basis, int n_electrons, const Eigen::MatrixXd& S, const Eigen::MatrixXd& Hcore, double E_nuc, const std::shared_ptr& jk, const vibeqc::Grid& xc_grid, const vibeqc::RKSOptions& options, const Eigen::MatrixXd& initial_density) { if (!jk) throw std::invalid_argument( "run_rks_scf_with_jk: jk_builder must not be None"); return vibeqc::run_rks_scf_with_jk( basis, n_electrons, S, Hcore, E_nuc, *jk, xc_grid, options, initial_density); }, py::arg("basis"), py::arg("n_electrons"), py::arg("S"), py::arg("Hcore"), py::arg("E_nuc"), py::arg("jk_builder"), py::arg("xc_grid"), py::arg("options") = vibeqc::RKSOptions{}, py::arg("initial_density") = Eigen::MatrixXd{}, py::call_guard(), "Lower-level closed-shell KS-DFT SCF entry point. Mirrors " "run_rhf_scf_with_jk but adds a pre-built ``xc_grid`` for " "the numerical XC integration. See run_rhf_scf_with_jk and " "RKSOptions.functional for the rest."); // ----- UKS (open-shell DFT) ------------------------------------------ py::class_(m, "UKSOptions") .def(py::init<>()) .def_readwrite("functional", &vibeqc::UKSOptions::functional) .def_readwrite("grid", &vibeqc::UKSOptions::grid) .def_readwrite("max_iter", &vibeqc::UKSOptions::max_iter) .def_readwrite("conv_tol_energy", &vibeqc::UKSOptions::conv_tol_energy) .def_readwrite("conv_tol_grad", &vibeqc::UKSOptions::conv_tol_grad) .def_readwrite("damping", &vibeqc::UKSOptions::damping) .def_readwrite("dynamic_damping", &vibeqc::UKSOptions::dynamic_damping, "See RHFOptions::dynamic_damping. Default False.") .def_readwrite("dynamic_damping_min", &vibeqc::UKSOptions::dynamic_damping_min) .def_readwrite("dynamic_damping_max", &vibeqc::UKSOptions::dynamic_damping_max) .def_readwrite("fock_mixing", &vibeqc::UKSOptions::fock_mixing, "Per-spin Kohn-Sham matrix mixing fraction. See " "RHFOptions.fock_mixing.") .def_readwrite("use_diis", &vibeqc::UKSOptions::use_diis) .def_readwrite("diis_start_iter", &vibeqc::UKSOptions::diis_start_iter) .def_readwrite("diis_subspace_size", &vibeqc::UKSOptions::diis_subspace_size) .def_readwrite("scf_accelerator", &vibeqc::UKSOptions::scf_accelerator, "SCF Fock-extrapolation accelerator. See " "RHFOptions.scf_accelerator and " "UHFOptions.scf_accelerator.") .def_readwrite("ediis_diis_switch_threshold", &vibeqc::UKSOptions::ediis_diis_switch_threshold, "See RHFOptions.ediis_diis_switch_threshold. " "Default 1e-1.") .def_readwrite("initial_guess", &vibeqc::UKSOptions::initial_guess) .def_readwrite("linear_dep_threshold", &vibeqc::UKSOptions::linear_dep_threshold) .def_readwrite("ecp_centers", &vibeqc::UKSOptions::ecp_centers, "List[ECPCenter]. See RHFOptions.ecp_centers.") .def_readwrite("ecp_library", &vibeqc::UKSOptions::ecp_library, "ECP XML library name. See RHFOptions.ecp_library.") .def_readwrite("level_shift", &vibeqc::UKSOptions::level_shift, "Phase C1a-2 — per-spin Saunders-Hillier level " "shift (Hartree). Same formula as UHF. See " "UHFOptions.level_shift. Default 0.0.") .def_readwrite("quadratic_fallback_iter", &vibeqc::UKSOptions::quadratic_fallback_iter, "Phase C1c — per-spin second-order SCF " "fallback. See RHFOptions.quadratic_fallback_iter. " "Default 0 (disabled).") .def_readwrite("quadratic_fallback_shift", &vibeqc::UKSOptions::quadratic_fallback_shift, "C1c orbital-rotation denominator shift. " "See RHFOptions.quadratic_fallback_shift.") .def_readwrite("quadratic_fallback_max_step", &vibeqc::UKSOptions::quadratic_fallback_max_step, "C1c per-step ‖κ‖ trust-region cap. " "See RHFOptions.quadratic_fallback_max_step.") .def_readwrite("newton_threshold", &vibeqc::UKSOptions::newton_threshold, "Phase D2c-KS-UHF Newton activation threshold. " "Default 0.0. XC kernel via polarised LDA fxc; " "raises on GGA functionals (Phase 17e pending).") .def_readwrite("newton_opts", &vibeqc::UKSOptions::newton_opts, "NewtonOptions. See NewtonOptions.") .def_readwrite("soscf_threshold", &vibeqc::UKSOptions::soscf_threshold, "Phase D2d-KS-UHF Neese SOSCF threshold. " "Default 0.0. No XC kernel needed.") .def_readwrite("soscf_opts", &vibeqc::UKSOptions::soscf_opts, "SOSCFOptions. See SOSCFOptions.") .def_readwrite("trah_threshold", &vibeqc::UKSOptions::trah_threshold, "Phase D2e-KS-UHF TRAH threshold. Default 0.0. " "Same XC kernel as Newton.") .def_readwrite("trah_opts", &vibeqc::UKSOptions::trah_opts, "TRAHOptions. See TRAHOptions.") .def_readwrite("density_fit", &vibeqc::UKSOptions::density_fit, "Enable density fitting (RI) for the per-spin " "Fock build. See RHFOptions.density_fit. " "Default False.") .def_readwrite("aux_basis", &vibeqc::UKSOptions::aux_basis, "Auxiliary basis name for density fitting. " "See RHFOptions.aux_basis.") .def_readwrite("scf_mode", &vibeqc::UKSOptions::scf_mode, "SCF Fock-build mode (SCFMode). See " "RHFOptions.scf_mode.") .def_readwrite("scf_mode_auto_threshold", &vibeqc::UKSOptions::scf_mode_auto_threshold, "AUTO cutoff (n_bf). See " "RHFOptions.scf_mode_auto_threshold.") .def_readwrite("schwarz_threshold", &vibeqc::UKSOptions::schwarz_threshold, "Per-quartet Schwarz skip bound for DIRECT mode. " "See RHFOptions.schwarz_threshold.") .def_readwrite("incremental_fock", &vibeqc::UKSOptions::incremental_fock, "Enable incremental ΔP Fock build for DIRECT. " "See RHFOptions.incremental_fock.") .def_readwrite("incremental_fock_reset_freq", &vibeqc::UKSOptions::incremental_fock_reset_freq, "Full-rebuild frequency for incremental_fock. " "See RHFOptions.incremental_fock_reset_freq.") .def_readwrite("schwarz_threshold_loose", &vibeqc::UKSOptions::schwarz_threshold_loose, "Two-phase Schwarz loose threshold. " "See RHFOptions.schwarz_threshold_loose.") .def_readwrite("schwarz_threshold_tighten_at", &vibeqc::UKSOptions::schwarz_threshold_tighten_at, "Grad-norm cutoff for the loose→tight transition. " "See RHFOptions.schwarz_threshold_tighten_at.") .def_readwrite("cosx", &vibeqc::UKSOptions::cosx, "Use chain-of-spheres exchange for K. See " "RKSOptions.cosx. Default False.") .def_readwrite("cosx_grid", &vibeqc::UKSOptions::cosx_grid, "GridOptions for the COSX grid. See " "RKSOptions.cosx_grid.") .def_readwrite("dft_plus_u_sites", &vibeqc::UKSOptions::dft_plus_u_sites, "Internal: list of _HubbardSiteCxx. Populated by " "the vibeqc.run_uks Python wrapper.") .def_readwrite("dft_plus_u_ao_groups", &vibeqc::UKSOptions::dft_plus_u_ao_groups, "Internal: parallel AO-index lists per site."); py::class_(m, "UKSResult") .def_readonly("energy", &vibeqc::UKSResult::energy) .def_readonly("e_electronic", &vibeqc::UKSResult::e_electronic) .def_readonly("e_coulomb", &vibeqc::UKSResult::e_coulomb) .def_readonly("e_hf_exchange",&vibeqc::UKSResult::e_hf_exchange) .def_readonly("e_dft_plus_u", &vibeqc::UKSResult::e_dft_plus_u, "Dudarev DFT+U contribution (Hartree). Sum of " "per-spin contributions; 0 when +U is not active.") .def_readonly("e_xc", &vibeqc::UKSResult::e_xc) .def_readonly("e_nuclear", &vibeqc::UKSResult::e_nuclear) .def_readonly("n_iter", &vibeqc::UKSResult::n_iter) .def_readonly("converged", &vibeqc::UKSResult::converged) .def_readonly("s_squared", &vibeqc::UKSResult::s_squared) .def_readonly("s_squared_ideal", &vibeqc::UKSResult::s_squared_ideal) .def_readonly("mo_energies_alpha", &vibeqc::UKSResult::mo_energies_alpha) .def_readonly("mo_coeffs_alpha", &vibeqc::UKSResult::mo_coeffs_alpha) .def_readonly("density_alpha",&vibeqc::UKSResult::density_alpha) .def_readonly("fock_alpha", &vibeqc::UKSResult::fock_alpha) .def_readonly("mo_energies_beta", &vibeqc::UKSResult::mo_energies_beta) .def_readonly("mo_coeffs_beta", &vibeqc::UKSResult::mo_coeffs_beta) .def_readonly("density_beta", &vibeqc::UKSResult::density_beta) .def_readonly("fock_beta", &vibeqc::UKSResult::fock_beta) .def_readonly("scf_trace", &vibeqc::UKSResult::scf_trace) .def_readonly("functional", &vibeqc::UKSResult::functional) .def("__repr__", [](const vibeqc::UKSResult& r) { return std::string{"UKSResult(energy="} + std::to_string(r.energy) + ", functional=" + r.functional + ", converged=" + (r.converged ? "True" : "False") + ", =" + std::to_string(r.s_squared) + ")"; }); m.def("run_uks", &vibeqc::run_uks, py::arg("molecule"), py::arg("basis"), py::arg("options") = vibeqc::UKSOptions{}, py::call_guard(), "Run unrestricted Kohn-Sham DFT on an open-shell molecule. " "Alpha/beta occupations follow from multiplicity."); m.def("run_uhf", &vibeqc::run_uhf, py::arg("molecule"), py::arg("basis"), py::arg("options") = vibeqc::UHFOptions{}, py::call_guard(), "Run unrestricted Hartree-Fock SCF on an open- or closed-shell " "molecule. Alpha/beta occupations follow from molecule.multiplicity. " "Returns a UHFResult."); m.def("run_uhf_scf_with_jk", [](const vibeqc::BasisSet& basis, int n_alpha, int n_beta, const Eigen::MatrixXd& S, const Eigen::MatrixXd& Hcore, double E_nuc, const std::shared_ptr& jk, const vibeqc::UHFOptions& options, const Eigen::MatrixXd& init_alpha, const Eigen::MatrixXd& init_beta) { if (!jk) throw std::invalid_argument( "run_uhf_scf_with_jk: jk_builder must not be None"); return vibeqc::run_uhf_scf_with_jk( basis, n_alpha, n_beta, S, Hcore, E_nuc, *jk, options, init_alpha, init_beta); }, py::arg("basis"), py::arg("n_alpha"), py::arg("n_beta"), py::arg("S"), py::arg("Hcore"), py::arg("E_nuc"), py::arg("jk_builder"), py::arg("options") = vibeqc::UHFOptions{}, py::arg("init_alpha") = Eigen::MatrixXd{}, py::arg("init_beta") = Eigen::MatrixXd{}, py::call_guard(), "Lower-level UHF SCF entry point. Mirrors run_rhf_scf_with_jk " "but with per-spin electron counts. ``init_alpha`` / " "``init_beta`` are optional density-matrix seeds (both supplied " "or both empty); empty falls back to a Hcore-diag guess."); m.def("run_uks_scf_with_jk", [](const vibeqc::BasisSet& basis, int n_alpha, int n_beta, const Eigen::MatrixXd& S, const Eigen::MatrixXd& Hcore, double E_nuc, const std::shared_ptr& jk, const vibeqc::Grid& xc_grid, const vibeqc::UKSOptions& options, const Eigen::MatrixXd& init_alpha, const Eigen::MatrixXd& init_beta) { if (!jk) throw std::invalid_argument( "run_uks_scf_with_jk: jk_builder must not be None"); return vibeqc::run_uks_scf_with_jk( basis, n_alpha, n_beta, S, Hcore, E_nuc, *jk, xc_grid, options, init_alpha, init_beta); }, py::arg("basis"), py::arg("n_alpha"), py::arg("n_beta"), py::arg("S"), py::arg("Hcore"), py::arg("E_nuc"), py::arg("jk_builder"), py::arg("xc_grid"), py::arg("options") = vibeqc::UKSOptions{}, py::arg("init_alpha") = Eigen::MatrixXd{}, py::arg("init_beta") = Eigen::MatrixXd{}, py::call_guard(), "Lower-level UKS SCF entry point. Mirrors run_rks_scf_with_jk " "but with per-spin electron counts and per-spin initial " "densities. The ``xc_grid`` is the numerical-integration grid " "for the XC piece — pass a periodic grid for periodic-Γ " "open-shell DFT."); // ----- MP2 (second-order Moller-Plesset) ----------------------------- py::class_(m, "MP2Options") .def(py::init<>()) .def_readwrite("density_fit", &vibeqc::MP2Options::density_fit, "Enable density fitting for the AO→MO ERI " "transformation. Replaces the four-index ERI " "and O(n^5) ao_to_mo_ovov transform with the " "DF factorisation (ia|jb) ≈ Σ_P B^P_ia B^P_jb " "(Vahtras-Almlöf-Feyereisen 1993). Default False.") .def_readwrite("aux_basis", &vibeqc::MP2Options::aux_basis, "Auxiliary basis name for density fitting " "(libint-recognised, e.g. \"def2-tzvp-rifit\", " "\"cc-pvtz-ri\"). Use " "vibeqc.default_aux_basis_for(orbital_basis_name, " "kind=\"ri\") for autodetection. Empty + " "density_fit=True raises.") .def_readwrite("c_os", &vibeqc::MP2Options::c_os, "Opposite-spin (singlet pair) scaling coefficient " "in the Grimme JCP 118, 9095 (2003) decomposition. " "Default 1.0 (canonical MP2). 6/5 for SCS-MP2, " "1.3 for SOS-MP2, 0.27 for the B2PLYP MP2 " "correction.") .def_readwrite("c_ss", &vibeqc::MP2Options::c_ss, "Same-spin (triplet pair) scaling coefficient. " "Default 1.0 (canonical MP2). 1/3 for SCS-MP2, " "0.0 for SOS-MP2, 0.27 for the B2PLYP MP2 " "correction.") .def_readwrite("use_float_intermediates", &vibeqc::MP2Options::use_float_intermediates, "When density_fit=True, store the half-transformed " "B-mo intermediate as single-precision float " "instead of double. Halves the memory for the " "(n_aux × n_occ·n_vir) tensor. Default False.") .def_readwrite("report_ri_residual", &vibeqc::MP2Options::report_ri_residual, "Diagnostic: when density_fit=True, also build " "the canonical four-index (ia|jb) tensor and " "report (e_os_RI − e_os_canonical) / " "(e_ss_RI − e_ss_canonical) on the MP2Result. " "Doubles the runtime cost — opt-in, intended " "for verifying aux-basis adequacy. Default False."); py::class_(m, "MP2Result") .def_readonly("e_hf", &vibeqc::MP2Result::e_hf) .def_readonly("e_correlation", &vibeqc::MP2Result::e_correlation, "Spin-component-scaled correlation energy: " "c_os · e_os + c_ss · e_ss. At c_os = c_ss = 1 " "this is canonical RMP2.") .def_readonly("e_total", &vibeqc::MP2Result::e_total) .def_readonly("e_os", &vibeqc::MP2Result::e_os, "Unscaled opposite-spin (αβ singlet pair) energy " "Σ (ia|jb)²/Δ (Grimme JCP 118, 9095 (2003)). " "Multiply by ``MP2Options.c_os`` for the scaled " "contribution to e_correlation.") .def_readonly("e_ss", &vibeqc::MP2Result::e_ss, "Unscaled same-spin (αα+ββ triplet pair) energy " "Σ [(ia|jb)² − (ia|jb)(ib|ja)]/Δ. Multiply by " "``MP2Options.c_ss`` for the scaled contribution " "to e_correlation.") .def_readonly("e_os_ri_residual", &vibeqc::MP2Result::e_os_ri_residual, "RI fit residual e_os_RI − e_os_canonical, " "populated only when ``MP2Options.density_fit && " "report_ri_residual``. 0.0 otherwise; check " "``ri_residual_reported`` to distinguish.") .def_readonly("e_ss_ri_residual", &vibeqc::MP2Result::e_ss_ri_residual, "RI fit residual e_ss_RI − e_ss_canonical, " "populated only when ``MP2Options.density_fit && " "report_ri_residual``.") .def_readonly("ri_residual_reported", &vibeqc::MP2Result::ri_residual_reported, "True iff ``report_ri_residual`` was set on the " "MP2Options for this run AND density_fit was on. " "Use this to distinguish 'residual not requested' " "from 'residual requested and happens to be 0'.") .def("__repr__", [](const vibeqc::MP2Result& r) { return std::string{"MP2Result(e_total="} + std::to_string(r.e_total) + ", e_corr=" + std::to_string(r.e_correlation) + ")"; }); m.def("run_mp2", &vibeqc::run_mp2, py::arg("molecule"), py::arg("basis"), py::arg("rhf_result"), py::arg("options") = vibeqc::MP2Options{}, py::call_guard(), "Run closed-shell RMP2 correlation on a converged RHF reference. " "Pass options.density_fit=True with options.aux_basis set to " "use density fitting (RI-MP2)."); // ----- UMP2 (open-shell MP2) ----------------------------------------- py::class_(m, "UMP2Options") .def(py::init<>()) .def_readwrite("density_fit", &vibeqc::UMP2Options::density_fit, "Enable density fitting for the AO→MO ERI " "transform. The αα / ββ / αβ OVOV tensors all " "route through the DF factorisation; α and β " "MO B-tensors are built once each and the three " "blocks form via single GEMMs. See " "MP2Options.density_fit. Default False.") .def_readwrite("aux_basis", &vibeqc::UMP2Options::aux_basis, "Auxiliary basis name for density fitting. " "See MP2Options.aux_basis.") .def_readwrite("c_os", &vibeqc::UMP2Options::c_os, "Opposite-spin (αβ channel) scaling coefficient. " "Default 1.0 (canonical UMP2). See " "``MP2Options.c_os`` for the named SCS / SOS / " "B2PLYP recipes.") .def_readwrite("c_ss", &vibeqc::UMP2Options::c_ss, "Same-spin (αα+ββ channels) scaling coefficient. " "Default 1.0 (canonical UMP2). See " "``MP2Options.c_ss``.") .def_readwrite("report_ri_residual", &vibeqc::UMP2Options::report_ri_residual, "Diagnostic: when density_fit=True, also build " "the canonical four-index OVOV tensors for all " "three (αα, ββ, αβ) channels and report the " "per-channel residual on the UMP2Result. Doubles " "the runtime cost — opt-in, intended for " "verifying aux-basis adequacy. See " "MP2Options.report_ri_residual. Default False."); py::class_(m, "UMP2Result") .def_readonly("e_hf", &vibeqc::UMP2Result::e_hf) .def_readonly("e_correlation", &vibeqc::UMP2Result::e_correlation) .def_readonly("e_total", &vibeqc::UMP2Result::e_total) .def_readonly("e_aa", &vibeqc::UMP2Result::e_aa, "αα same-spin channel") .def_readonly("e_bb", &vibeqc::UMP2Result::e_bb, "ββ same-spin channel") .def_readonly("e_ab", &vibeqc::UMP2Result::e_ab, "αβ opposite-spin channel") .def_readonly("e_aa_ri_residual", &vibeqc::UMP2Result::e_aa_ri_residual, "RI fit residual e_aa_RI − e_aa_canonical, " "populated only when ``UMP2Options.density_fit && " "report_ri_residual``. 0.0 otherwise; check " "``ri_residual_reported`` to distinguish.") .def_readonly("e_bb_ri_residual", &vibeqc::UMP2Result::e_bb_ri_residual, "RI fit residual e_bb_RI − e_bb_canonical, " "populated only when ``UMP2Options.density_fit && " "report_ri_residual``.") .def_readonly("e_ab_ri_residual", &vibeqc::UMP2Result::e_ab_ri_residual, "RI fit residual e_ab_RI − e_ab_canonical, " "populated only when ``UMP2Options.density_fit && " "report_ri_residual``.") .def_readonly("ri_residual_reported", &vibeqc::UMP2Result::ri_residual_reported, "True iff ``report_ri_residual`` was set on the " "UMP2Options for this run AND density_fit was on.") .def("__repr__", [](const vibeqc::UMP2Result& r) { return std::string{"UMP2Result(e_total="} + std::to_string(r.e_total) + ", e_corr=" + std::to_string(r.e_correlation) + ")"; }); m.def("run_ump2", &vibeqc::run_ump2, py::arg("molecule"), py::arg("basis"), py::arg("uhf_result"), py::arg("options") = vibeqc::UMP2Options{}, py::call_guard(), "Run open-shell UMP2 correlation on a converged UHF reference. " "Pass options.density_fit=True with options.aux_basis set to " "use density fitting (RI-UMP2)."); // ----- CCSD / CCSD(T) (coupled-cluster) ----------------------------- py::class_(m, "CCSDOptions") .def(py::init<>()) .def_readwrite("density_fit", &vibeqc::CCSDOptions::density_fit, "Enable density fitting for the two-electron integrals. " "Default True — DF is the standard approach for CC.") .def_readwrite("aux_basis", &vibeqc::CCSDOptions::aux_basis, "Auxiliary basis name for density fitting " "(libint-recognised, e.g. \"def2-tzvp-rifit\"). " "Empty + density_fit=True raises.") .def_readwrite("max_iter", &vibeqc::CCSDOptions::max_iter, "Maximum CCSD iterations. Default 100.") .def_readwrite("conv_tol_energy", &vibeqc::CCSDOptions::conv_tol_energy, "Convergence threshold on energy change (Ha). Default 1e-8.") .def_readwrite("conv_tol_residual", &vibeqc::CCSDOptions::conv_tol_residual, "Convergence threshold on ||R1|| + ||R2||. Default 1e-7.") .def_readwrite("diis_subspace_size", &vibeqc::CCSDOptions::diis_subspace_size, "DIIS subspace size for amplitude extrapolation. " "Set to 0 to disable DIIS. Default 6.") .def_readwrite("n_frozen_core", &vibeqc::CCSDOptions::n_frozen_core, "Number of lowest-energy occupied MOs to freeze. " "Default 0 (all-electron correlation).") .def_readwrite("compute_triples", &vibeqc::CCSDOptions::compute_triples, "Compute the perturbative (T) correction after CCSD " "converges. Default True.") .def_readwrite("triples_memory_mode", &vibeqc::CCSDOptions::triples_memory_mode, "Triples memory strategy: \"fast\" (default) or " "\"low\" (aggressive batching, lower memory)."); py::class_(m, "CCSDIteration", "One row of the CCSD iteration trace.") .def_readonly("iter", &vibeqc::CCSDIteration::iter) .def_readonly("energy", &vibeqc::CCSDIteration::energy) .def_readonly("delta_e", &vibeqc::CCSDIteration::delta_e) .def_readonly("r1_norm", &vibeqc::CCSDIteration::r1_norm) .def_readonly("r2_norm", &vibeqc::CCSDIteration::r2_norm) .def_readonly("diis_subspace", &vibeqc::CCSDIteration::diis_subspace); py::class_(m, "CCSDResult") .def_readonly("e_hf", &vibeqc::CCSDResult::e_hf, "Underlying RHF total energy (Hartree).") .def_readonly("e_ccsd_correlation", &vibeqc::CCSDResult::e_ccsd_correlation, "CCSD correlation energy contribution.") .def_readonly("e_ccsd", &vibeqc::CCSDResult::e_ccsd, "e_hf + e_ccsd_correlation.") .def_readonly("e_t", &vibeqc::CCSDResult::e_t, "Perturbative (T) correction (0 if not computed).") .def_readonly("e_ccsd_t", &vibeqc::CCSDResult::e_ccsd_t, "e_ccsd + e_t.") .def_readonly("e_total", &vibeqc::CCSDResult::e_total, "Total CCSD(T) energy (alias for e_ccsd_t).") .def_readonly("n_iter", &vibeqc::CCSDResult::n_iter) .def_readonly("converged", &vibeqc::CCSDResult::converged) .def_readonly("t1_norm", &vibeqc::CCSDResult::t1_norm, "Frobenius norm of converged T1 amplitudes.") .def_readonly("t2_norm", &vibeqc::CCSDResult::t2_norm, "Frobenius norm of converged T2 amplitudes.") .def_readonly("cc_trace", &vibeqc::CCSDResult::cc_trace, "Per-iteration CCSD trace.") .def("__repr__", [](const vibeqc::CCSDResult& r) { return std::string{"CCSDResult(e_total="} + std::to_string(r.e_total) + ", e_corr=" + std::to_string(r.e_ccsd_correlation) + ", e_t=" + std::to_string(r.e_t) + ")"; }); m.def("run_ccsd", &vibeqc::run_ccsd, py::arg("molecule"), py::arg("basis"), py::arg("rhf_result"), py::arg("options") = vibeqc::CCSDOptions{}, py::call_guard(), "Run closed-shell DF-CCSD (and optionally DF-CCSD(T)) on a " "converged RHF reference. Requires density_fit=True with " "an auxiliary basis set."); // ----- Crystal / space-group analysis (spglib) ----------------------- py::class_(m, "Crystal") .def(py::init<>()) .def(py::init([](const Eigen::Matrix3d& lattice, const Eigen::Matrix3Xd& fractional_coords, const std::vector& species) { vibeqc::Crystal c; c.lattice = lattice; c.fractional_coords = fractional_coords; c.species = species; if (static_cast(fractional_coords.cols()) != static_cast(species.size())) { throw std::runtime_error( "Crystal: fractional_coords must have one column " "per species entry"); } return c; }), py::arg("lattice"), py::arg("fractional_coords"), py::arg("species"), "Construct a Crystal. ``lattice`` columns are the Cartesian " "lattice vectors in bohr; ``fractional_coords`` is a 3×N matrix " "with one atom per column.") .def_readwrite("lattice", &vibeqc::Crystal::lattice) .def_readwrite("fractional_coords", &vibeqc::Crystal::fractional_coords) .def_readwrite("species", &vibeqc::Crystal::species) .def_property_readonly("n_atoms", &vibeqc::Crystal::n_atoms) .def("__repr__", [](const vibeqc::Crystal& c) { return std::string{"Crystal(n_atoms="} + std::to_string(c.n_atoms()) + ")"; }); py::class_(m, "SymmetryOp") .def_readonly("rotation", &vibeqc::SymmetryOp::rotation, "Integer rotation matrix in the fractional basis.") .def_readonly("translation", &vibeqc::SymmetryOp::translation, "Translation in fractional coordinates.") .def("__repr__", [](const vibeqc::SymmetryOp&) { return std::string{"SymmetryOp(...)"}; }); py::class_(m, "SpaceGroup") .def_readonly("number", &vibeqc::SpaceGroup::number) .def_readonly("international_symbol", &vibeqc::SpaceGroup::international_symbol) .def_readonly("hall_number", &vibeqc::SpaceGroup::hall_number) .def_readonly("point_group", &vibeqc::SpaceGroup::point_group) .def_readonly("operations", &vibeqc::SpaceGroup::operations) .def_readonly("equivalent_atoms", &vibeqc::SpaceGroup::equivalent_atoms) .def_property_readonly("order", [](const vibeqc::SpaceGroup& sg) { return static_cast(sg.operations.size()); }) .def("__repr__", [](const vibeqc::SpaceGroup& sg) { return std::string{"SpaceGroup(number="} + std::to_string(sg.number) + ", symbol='" + sg.international_symbol + "', order=" + std::to_string(sg.operations.size()) + ")"; }); m.def("analyze_symmetry", &vibeqc::analyze, py::arg("crystal"), py::arg("symprec") = 1.0e-5, py::call_guard(), "Run spglib space-group analysis on a Crystal."); m.def("to_primitive", &vibeqc::to_primitive, py::arg("crystal"), py::arg("symprec") = 1.0e-5, py::call_guard(), "Return the primitive cell of an input Crystal."); m.def("spglib_version", &vibeqc::spglib_version, "Version of the linked spglib library."); py::class_(m, "IrreducibleKMesh") .def_readonly("fractional_kpoints", &vibeqc::IrreducibleKMesh::fractional_kpoints) .def_readonly("weights", &vibeqc::IrreducibleKMesh::weights) .def_readonly("ir_mapping", &vibeqc::IrreducibleKMesh::ir_mapping); m.def("irreducible_kpoints", &vibeqc::irreducible_kpoints, py::arg("crystal"), py::arg("mesh"), py::arg("is_shift") = std::array{0, 0, 0}, py::arg("symprec") = 1.0e-5, py::call_guard(), "Irreducible-Brillouin-zone k-point list for a Monkhorst–Pack mesh " "(fractional coordinates in the reciprocal lattice)."); // ----- Periodic systems + lattice sums ------------------------------- py::class_(m, "PeriodicSystem") .def(py::init<>()) .def(py::init([](int dim, const Eigen::Matrix3d& lattice, const std::vector& unit_cell, int charge, int multiplicity) { vibeqc::PeriodicSystem s; s.dim = dim; s.lattice = lattice; s.unit_cell = unit_cell; s.charge = charge; s.multiplicity = multiplicity; if (dim < 1 || dim > 3) { throw std::runtime_error( "PeriodicSystem: dim must be 1, 2, or 3"); } return s; }), py::arg("dim"), py::arg("lattice"), py::arg("unit_cell"), py::arg("charge") = 0, py::arg("multiplicity") = 1, "Construct a PeriodicSystem. ``lattice`` columns are Cartesian " "lattice vectors in bohr (always 3×3 regardless of ``dim``; " "columns beyond ``dim`` are implicit vacuum directions).") .def_readwrite("dim", &vibeqc::PeriodicSystem::dim) .def_readwrite("lattice", &vibeqc::PeriodicSystem::lattice) .def_readwrite("unit_cell", &vibeqc::PeriodicSystem::unit_cell) .def_readwrite("charge", &vibeqc::PeriodicSystem::charge) .def_readwrite("multiplicity", &vibeqc::PeriodicSystem::multiplicity) .def_readwrite("symmetry", &vibeqc::PeriodicSystem::symmetry) .def("reciprocal_lattice", &vibeqc::PeriodicSystem::reciprocal_lattice) .def("n_electrons", &vibeqc::PeriodicSystem::n_electrons) .def("unit_cell_molecule", &vibeqc::PeriodicSystem::unit_cell_molecule) .def("__repr__", [](const vibeqc::PeriodicSystem& s) { return std::string{"PeriodicSystem(dim="} + std::to_string(s.dim) + ", n_atoms=" + std::to_string(s.unit_cell.size()) + ")"; }); m.def("attach_symmetry", &vibeqc::attach_symmetry, py::arg("system"), py::arg("symprec") = 1.0e-5, py::call_guard(), "Populate system.symmetry with a spglib analysis of the unit cell."); py::class_(m, "LatticeCell") .def_readonly("index", &vibeqc::LatticeCell::index) .def_readonly("r_cart", &vibeqc::LatticeCell::r_cart); py::enum_(m, "CoulombMethod") .value("DIRECT_TRUNCATED", vibeqc::CoulombMethod::DIRECT_TRUNCATED) .value("EWALD_3D", vibeqc::CoulombMethod::EWALD_3D) .value("SLAB_EWALD_2D", vibeqc::CoulombMethod::SLAB_EWALD_2D) .value("NEUTRALIZED_1D", vibeqc::CoulombMethod::NEUTRALIZED_1D); py::class_(m, "LatticeSumOptions") .def(py::init<>()) .def_readwrite("cutoff_bohr", &vibeqc::LatticeSumOptions::cutoff_bohr) .def_readwrite("nuclear_cutoff_bohr", &vibeqc::LatticeSumOptions::nuclear_cutoff_bohr) .def_readwrite("coulomb_method", &vibeqc::LatticeSumOptions::coulomb_method) .def_readwrite("screening_overlap_threshold", &vibeqc::LatticeSumOptions::screening_overlap_threshold) .def_readwrite("screening_exchange_threshold", &vibeqc::LatticeSumOptions::screening_exchange_threshold) .def_readwrite("schwarz_threshold", &vibeqc::LatticeSumOptions::schwarz_threshold, "Schwarz screening threshold for the periodic 2-e " "Fock build. Default 1e-10 Ha. Set to 0.0 to disable " "(unscreened, O(n_c^3 * n_shells^4) — very slow).") .def_readwrite("schwarz_threshold_forces", &vibeqc::LatticeSumOptions::schwarz_threshold_forces, "Schwarz screening threshold for the periodic 2-e " "Fock derivative pass (gradient ERIs). Default " "1e-14 — 100× tighter than schwarz_threshold " "because plain Schwarz on derivatives is non-rigorous " "(integral bound is not a derivative bound). Set to " "0.0 to disable."); py::class_(m, "LatticeMatrixSet") .def(py::init<>()) .def_readonly("nbf", &vibeqc::LatticeMatrixSet::nbf) .def_readonly("cells", &vibeqc::LatticeMatrixSet::cells) .def_readonly("blocks", &vibeqc::LatticeMatrixSet::blocks, "Per-cell ``(n_bf, n_bf)`` block matrices. " "Note: ``.blocks`` returns a fresh Python list " "each access (pybind11's default for " "std::vector); element assignment to " "the returned list does NOT propagate to C++. " "Use ``set_block(i, M)`` to mutate the underlying " "C++ vector in place.") .def("set_block", [](vibeqc::LatticeMatrixSet& s, int i, const Eigen::MatrixXd& M) { if (i < 0 || static_cast(i) >= s.blocks.size()) { throw py::index_error( "LatticeMatrixSet.set_block: index " + std::to_string(i) + " out of range (size " + std::to_string(s.blocks.size()) + ")"); } if (static_cast(M.rows()) != s.nbf || static_cast(M.cols()) != s.nbf) { throw py::value_error( "LatticeMatrixSet.set_block: matrix shape " "must equal (nbf, nbf)"); } s.blocks[static_cast(i)] = M; }, py::arg("index"), py::arg("matrix"), "Replace the block at the given cell index with ``matrix``. " "Mutates the underlying C++ storage in place — unlike " "``self.blocks[i] = matrix`` which writes to a transient " "Python-list copy and silently no-ops at the C++ level.") .def("__len__", [](const vibeqc::LatticeMatrixSet& s) { return s.size(); }); m.def("direct_lattice_cells", &vibeqc::direct_lattice_cells, py::arg("system"), py::arg("cutoff_bohr"), py::call_guard(), "Enumerate integer lattice cells within a Cartesian cutoff."); m.def("compute_overlap_lattice", &vibeqc::compute_overlap_lattice, py::arg("basis"), py::arg("system"), py::arg("options"), py::call_guard(), "Lattice-summed overlap integrals S_μν(g)."); m.def("compute_kinetic_lattice", &vibeqc::compute_kinetic_lattice, py::arg("basis"), py::arg("system"), py::arg("options"), py::call_guard(), "Lattice-summed kinetic-energy integrals T_μν(g)."); m.def("compute_nuclear_lattice", &vibeqc::compute_nuclear_lattice, py::arg("basis"), py::arg("system"), py::arg("options"), py::call_guard(), "Lattice-summed nuclear-attraction integrals V_μν(g)."); m.def("compute_overlap_lattice_explicit", &vibeqc::compute_overlap_lattice_explicit, py::arg("basis"), py::arg("system"), py::arg("cells"), py::call_guard(), "Overlap S(g) on caller-supplied cell list (Phase SYM3b)."); m.def("compute_kinetic_lattice_explicit", &vibeqc::compute_kinetic_lattice_explicit, py::arg("basis"), py::arg("system"), py::arg("cells"), py::call_guard(), "Kinetic T(g) on caller-supplied cell list (Phase SYM3b)."); m.def("compute_nuclear_lattice_explicit", &vibeqc::compute_nuclear_lattice_explicit, py::arg("basis"), py::arg("system"), py::arg("options"), py::arg("cells"), py::call_guard(), "Nuclear V(g) on caller-supplied cell list (Phase SYM3b)."); m.def("compute_nuclear_erfc_lattice", &vibeqc::compute_nuclear_erfc_lattice, py::arg("basis"), py::arg("system"), py::arg("omega"), py::arg("options"), py::call_guard(), "Short-range (erfc-screened) component of the nuclear-attraction " "lattice sum. Ewald building block: exponentially convergent real-" "space sum for any ω > 0. Combine with the reciprocal-space long-" "range part (Phase 12e-c) for the full 3D Ewald treatment of V(g)."); // ----- Lattice multipole moments (BIPOLE Phase 1) --------------------- py::class_(m, "LatticeMultipoleSet") .def_readonly("nbf", &vibeqc::LatticeMultipoleSet::nbf) .def_readonly("L_max", &vibeqc::LatticeMultipoleSet::L_max) .def_readonly("cells", &vibeqc::LatticeMultipoleSet::cells) .def_readonly("origin", &vibeqc::LatticeMultipoleSet::origin) .def_readonly("blocks", &vibeqc::LatticeMultipoleSet::blocks, "Per-cell, per-component nbf×nbf moment matrix. " "Layout: blocks[c][component] is the (nbf, nbf) " "matrix of ⟨μ_0 | x^i y^j z^k | ν_c⟩ where (i, j, k) " "encodes the Cartesian multipole order. Components " "follow libint emultipole{1,2,3} ordering: " "[S, μ_x, μ_y, μ_z, Q_xx, Q_xy, Q_xz, Q_yy, Q_yz, Q_zz, ...].") .def("__len__", [](const vibeqc::LatticeMultipoleSet& s) { return s.cells.size(); }); m.def("cartesian_multipole_n_components", &vibeqc::cartesian_multipole_n_components, py::arg("L_max"), "Number of Cartesian-multipole components for moments up " "through L_max (libint convention): emultipole1→4, " "emultipole2→10, emultipole3→20."); m.def("compute_multipole_moments_lattice", &vibeqc::compute_multipole_moments_lattice, py::arg("basis"), py::arg("system"), py::arg("options"), py::arg("L_max") = 2, py::arg("origin") = std::array{0.0, 0.0, 0.0}, py::call_guard(), "Compute shell-pair Cartesian multipole moments at every " "lattice cell out to options.cutoff_bohr. Foundation for the " "BIPOLE multipole-far-pair branch. L_max ∈ {1, 2, 3} maps to " "libint emultipole{1, 2, 3}. The single expansion origin " "(typically (0,0,0)) is shared across all shell pairs; per-pair " "Gaussian-product centers are constructed via polynomial-shift " "in Python (see vibeqc.bipole_multipole)."); // ----- Spheropole kernel (EXT EL-SPHEROPOLE) ------------------------- py::class_(m, "SpheropolePairData") .def_readonly("a_prim", &vibeqc::SpheropolePairData::a_prim) .def_readonly("b_prim", &vibeqc::SpheropolePairData::b_prim) .def_readonly("l_a", &vibeqc::SpheropolePairData::l_a) .def_readonly("l_b", &vibeqc::SpheropolePairData::l_b) .def_readonly("gamma", &vibeqc::SpheropolePairData::gamma) .def_readonly("S0", &vibeqc::SpheropolePairData::S0) .def_readonly("P", &vibeqc::SpheropolePairData::P) .def_readonly("M1", &vibeqc::SpheropolePairData::M1) .def_readonly("M2", &vibeqc::SpheropolePairData::M2); m.def("compute_spheropole_bare_integrals", &vibeqc::compute_spheropole_bare_integrals, py::arg("basis"), py::arg("system"), py::arg("options"), py::call_guard(), "Compute bare Gaussian-product integrals for every primitive " "shell pair at every lattice cell. Returns a list-of-lists: " "one vector of SpheropolePairData per lattice cell."); m.def("compute_spheropole_kernel_lattice", &vibeqc::compute_spheropole_kernel_lattice, py::arg("prim_basis"), py::arg("system"), py::arg("options"), py::arg("orig_nbf"), py::arg("c_map"), py::call_guard(), "Compute the full EXT EL-SPHEROPOLE kernel K(g) for every " "lattice cell, ready for tracing with P_real. Uses analytic " "Gaussian-product formulas with the bond-symmetrised second-" "moment operator and CJAT33 per-(l_a,l_b) decomposition. " "The global TOTCHR/(2V·π^{3/2}) factor is NOT included — " "multiply after contraction with P."); // ----- Bloch sums & k-mesh ------------------------------------------- py::class_(m, "BlochKMesh") .def_readonly("mesh", &vibeqc::BlochKMesh::mesh) .def_readonly("is_shift", &vibeqc::BlochKMesh::is_shift) .def_readonly("kpoints", &vibeqc::BlochKMesh::kpoints) .def_readonly("weights", &vibeqc::BlochKMesh::weights) .def_readonly("ir_mapping", &vibeqc::BlochKMesh::ir_mapping) .def("__len__", [](const vibeqc::BlochKMesh& m) { return m.size(); }); m.def("monkhorst_pack", &vibeqc::monkhorst_pack, py::arg("system"), py::arg("mesh"), py::arg("is_shift") = std::array{0, 0, 0}, py::arg("use_symmetry") = false, py::call_guard(), "Monkhorst–Pack k-point mesh, optionally reduced to the IBZ."); // K4 — factory for building a BlochKMesh from an explicit Cartesian // k-point list + matching weights. Used by KPoints.to_bloch_kmesh() // when kind ∈ {"explicit", "band-path"} so user-supplied k-lists // and HPKOT band paths can feed any periodic SCF driver that // expects a BlochKMesh. m.def("bloch_kmesh_from_lists", [](const std::vector& kpoints_cart, const std::vector& weights) { if (kpoints_cart.size() != weights.size()) { throw std::runtime_error( "bloch_kmesh_from_lists: kpoints_cart and weights " "must have the same length."); } vibeqc::BlochKMesh bm; bm.mesh = {1, 1, 1}; bm.is_shift = {0, 0, 0}; bm.kpoints = kpoints_cart; bm.weights = weights; bm.ir_mapping.clear(); return bm; }, py::arg("kpoints_cart"), py::arg("weights"), "Build a BlochKMesh from explicit (Cartesian k-vector, weight) " "lists. mesh / is_shift / ir_mapping are left as 1×1×1 / 0 / " "empty since those are MP-specific concepts that don't apply " "to explicit lists or band paths."); m.def("bloch_sum", &vibeqc::bloch_sum, py::arg("real_space"), py::arg("k_cart"), py::call_guard(), "Bloch sum M(k) = Σ_g exp(i k·g) M(g)."); py::class_(m, "BandDiag") .def_readonly("energies", &vibeqc::BandDiag::energies) .def_readonly("coefficients", &vibeqc::BandDiag::coefficients); m.def("diagonalize_bloch", &vibeqc::diagonalize_bloch, py::arg("F_k"), py::arg("S_k"), py::arg("lindep_threshold") = 1.0e-9, py::call_guard(), "Solve F·C = S·C·diag(ε) at a single k-point by symmetric " "orthogonalisation."); // ----- JKBuilder (J / K Fock-build strategy) ------------------------- // Polymorphic interface used by the molecular SCF drivers and (via // the periodic factory below) the periodic Γ-only adapter. Concrete // factories return std::unique_ptr; pybind11 holds them // as ``std::shared_ptr`` for Python-side ownership convenience, and // we expose only build_J / build_K / build_g_rhf — the abstract // class doesn't need a Python __init__ since it's never directly // constructible. py::class_>( m, "JKBuilder") .def("build_J", &vibeqc::JKBuilder::build_J, py::arg("density"), py::call_guard(), "Coulomb matrix J(D)_{μν} = Σ_{λρ} D_{λρ} (μν|λρ).") .def("build_K", &vibeqc::JKBuilder::build_K, py::arg("density"), py::call_guard(), "Exchange matrix K(D)_{μν} = Σ_{λρ} D_{λρ} (μλ|νρ).") .def("build_g_rhf", &vibeqc::JKBuilder::build_g_rhf, py::arg("density"), py::arg("alpha_hf") = 1.0, py::call_guard(), "Closed-shell Fock 2-electron piece G(D) = J − ½·α_HF·K. " "Some implementations (DF, COSX-RIJCOSX, periodic-Γ) fuse " "the J and K kernels in one pass and override this; the " "default falls back to build_J + build_K."); m.def("make_four_index_jk_builder", [](const vibeqc::BasisSet& basis) { return std::shared_ptr( vibeqc::make_four_index_jk_builder(basis)); }, py::arg("basis"), "Direct four-index ERI JKBuilder. Materialises the full " "(μν|λρ) tensor at construction; reuses it across SCF " "iterations."); m.def("make_direct_jk_builder", [](const vibeqc::BasisSet& basis, double schwarz_threshold, bool incremental, int reset_freq) { return std::shared_ptr( vibeqc::make_direct_jk_builder( basis, schwarz_threshold, incremental, reset_freq)); }, py::arg("basis"), py::arg("schwarz_threshold") = 1e-10, py::arg("incremental") = false, py::arg("reset_freq") = 8, "Direct (integral-driven) JKBuilder. Caches the Schwarz Q " "matrix at construction; per SCF iter, walks an 8-fold-" "symmetric shell-quartet loop and evaluates only the " "surviving quartets via libint. O(n_shells² + n_bf²) " "memory instead of O(n_bf⁴) — closes the ~50-atom / " "def2-SVP OOM wall."); m.def("make_df_jk_builder", [](const vibeqc::BasisSet& basis, const vibeqc::BasisSet& aux) { return std::shared_ptr( vibeqc::make_df_jk_builder(basis, aux)); }, py::arg("basis"), py::arg("aux"), "Density-fitting (RIJK) JKBuilder. Owns a DensityFitting " "instance with the V-metric Cholesky and B-tensor."); m.def("make_cosx_jk_builder", [](const vibeqc::BasisSet& basis, const vibeqc::BasisSet& aux, const vibeqc::Grid& cosx_grid) { return std::shared_ptr( vibeqc::make_cosx_jk_builder(basis, aux, cosx_grid)); }, py::arg("basis"), py::arg("aux"), py::arg("cosx_grid"), "RIJCOSX JKBuilder. RI-J for the Coulomb piece, " "seminumerical chain-of-spheres for the K piece on the " "supplied grid."); m.def("make_periodic_gamma_jk_builder", [](const vibeqc::BasisSet& basis, const vibeqc::PeriodicSystem& system, const vibeqc::LatticeSumOptions& opts, double omega) { return std::shared_ptr( vibeqc::make_periodic_gamma_jk_builder( basis, system, opts, omega)); }, py::arg("basis"), py::arg("system"), py::arg("options"), py::arg("omega") = 0.0, "Γ-only molecular-limit periodic JKBuilder. Wraps " "build_jk_gamma_molecular_limit so the molecular SCF " "loop body can drive a periodic-Γ calculation through " "the same JKBuilder interface."); // ----- Γ-only periodic RHF ------------------------------------------- py::class_(m, "JKMatrices") .def_readonly("J", &vibeqc::JKMatrices::J) .def_readonly("K", &vibeqc::JKMatrices::K); m.def("build_jk_gamma_molecular_limit", &vibeqc::build_jk_gamma_molecular_limit, py::arg("basis"), py::arg("system"), py::arg("options"), py::arg("density"), py::arg("omega") = 0.0, py::call_guard(), "Direct-SCF J and K matrices at \u0393 for a \u0393-only molecular-limit " "density. Set ``omega > 0`` to use the erfc-screened Coulomb " "kernel (short-range Ewald split); default ``omega=0`` uses " "the full 1/r_12 kernel."); m.def("build_jk_gamma_molecular_limit_explicit", &vibeqc::build_jk_gamma_molecular_limit_explicit, py::arg("basis"), py::arg("system"), py::arg("cells"), py::arg("options"), py::arg("density"), py::arg("omega") = 0.0, py::call_guard(), "Same as build_jk_gamma_molecular_limit, but with caller-supplied " "cell list (Phase SYM3b — symmetry-reduced Fock build)."); // ----- Phase M3b — per-cell-pair J/K contributions ------------------- py::class_(m, "PairJKContribution") .def_readonly("c_g", &vibeqc::PairJKContribution::c_g) .def_readonly("c_p", &vibeqc::PairJKContribution::c_p) .def_readonly("J_contrib", &vibeqc::PairJKContribution::J_contrib) .def_readonly("K_contrib", &vibeqc::PairJKContribution::K_contrib); m.def("build_jk_pair_contributions", &vibeqc::build_jk_pair_contributions, py::arg("basis"), py::arg("system"), py::arg("cells"), py::arg("pairs"), py::arg("options"), py::arg("density"), py::arg("omega") = 0.0, py::call_guard(), "Per-pair J/K contributions for symmetry-reduced Fock build (M3b). " "Runs the same Γ-only molecular-limit kernel as " "build_jk_gamma_molecular_limit_explicit but stores each (c_g, c_p) " "cell pair's J/K block separately instead of summing. Summing over " "the full n_c × n_c pair set reproduces the explicit kernel's J/K."); m.def("nuclear_repulsion_per_cell", &vibeqc::nuclear_repulsion_per_cell, py::arg("system"), py::arg("options"), py::call_guard(), "Nuclear-nuclear repulsion per unit cell. Dispatches on " "``options.coulomb_method``: DIRECT_TRUNCATED (half-summed " "direct truncation — conditionally convergent in 3D) or " "EWALD_3D (Ewald summation — absolutely convergent)."); // ----- Ewald summation (Phase 12e-a) --------------------------------- py::class_(m, "EwaldOptions") .def(py::init<>()) .def_readwrite("alpha", &vibeqc::EwaldOptions::alpha, "Gaussian screening parameter; ≤ 0 auto-selects.") .def_readwrite("real_cutoff_bohr", &vibeqc::EwaldOptions::real_cutoff_bohr) .def_readwrite("recip_cutoff_bohr_inv", &vibeqc::EwaldOptions::recip_cutoff_bohr_inv, "Reciprocal-space cutoff in bohr⁻¹; ≤ 0 auto-selects.") .def_readwrite("tolerance", &vibeqc::EwaldOptions::tolerance); m.def("ewald_point_charge_energy", &vibeqc::ewald_point_charge_energy, py::arg("lattice"), py::arg("positions_cart"), py::arg("charges"), py::arg("options") = vibeqc::EwaldOptions{}, py::call_guard(), "Ewald-summed Madelung energy for an arbitrary 3D lattice of " "point charges. Handles non-neutral cells via the jellium " "(uniform compensating background) convention."); m.def("ewald_nuclear_repulsion", &vibeqc::ewald_nuclear_repulsion, py::arg("system"), py::arg("options") = vibeqc::EwaldOptions{}, py::call_guard(), "Nuclear Madelung energy for a PeriodicSystem via Ewald summation " "(dim == 3 only)."); m.def("ewald_point_charge_potential", &vibeqc::ewald_point_charge_potential, py::arg("lattice"), py::arg("charge_positions_cart"), py::arg("charges"), py::arg("eval_points_cart"), py::arg("options") = vibeqc::EwaldOptions{}, py::arg("include_short_range") = true, py::call_guard(), "Ewald-summed Coulomb potential of a periodic lattice of point " "charges, evaluated at a list of Cartesian points. Set " "``include_short_range=False`` to omit the erfc 1/r spikes at " "each charge — used when the short-range part will be added " "analytically via libint's erfc_nuclear."); m.def("ewald_nuclear_potential", &vibeqc::ewald_nuclear_potential, py::arg("system"), py::arg("eval_points_cart"), py::arg("options") = vibeqc::EwaldOptions{}, py::arg("include_short_range") = true, py::call_guard(), "Electronic nuclear-attraction potential V_nuc(r) = -Σ Z_A/|r-R_A| " "summed over a 3D-periodic lattice via Ewald."); m.def("compute_nuclear_lattice_ewald", &vibeqc::compute_nuclear_lattice_ewald, py::arg("basis"), py::arg("system"), py::arg("grid"), py::arg("options"), py::arg("ewald_options") = vibeqc::EwaldOptions{}, py::call_guard(), "Full 3D Ewald-summed nuclear-attraction lattice sum V_μν(g) via " "numerical integration on the molecular grid. Replacement for " "compute_nuclear_lattice in 3D bulk where the point-charge sum " "is conditionally convergent."); py::class_(m, "PeriodicRHFOptions") .def(py::init<>()) .def_readwrite("max_iter", &vibeqc::PeriodicRHFOptions::max_iter) .def_readwrite("conv_tol_energy", &vibeqc::PeriodicRHFOptions::conv_tol_energy) .def_readwrite("conv_tol_grad", &vibeqc::PeriodicRHFOptions::conv_tol_grad) .def_readwrite("damping", &vibeqc::PeriodicRHFOptions::damping) .def_readwrite("dynamic_damping", &vibeqc::PeriodicRHFOptions::dynamic_damping, "See RHFOptions.dynamic_damping. Default False.") .def_readwrite("dynamic_damping_min", &vibeqc::PeriodicRHFOptions::dynamic_damping_min) .def_readwrite("dynamic_damping_max", &vibeqc::PeriodicRHFOptions::dynamic_damping_max) .def_readwrite("fock_mixing", &vibeqc::PeriodicRHFOptions::fock_mixing, "Fock matrix mixing fraction. CRYSTAL FMIXING " "30 corresponds to 0.30. Default 0.0 (off).") .def_readwrite("use_diis", &vibeqc::PeriodicRHFOptions::use_diis) .def_readwrite("diis_start_iter", &vibeqc::PeriodicRHFOptions::diis_start_iter) .def_readwrite("diis_subspace_size", &vibeqc::PeriodicRHFOptions::diis_subspace_size) .def_readwrite("scf_accelerator", &vibeqc::PeriodicRHFOptions::scf_accelerator, "SCFAccelerator family member (DIIS / KDIIS / " "EDIIS / EDIIS_DIIS / ADIIS). Default DIIS. " "Wired for the Γ-only RHF driver; multi-k " "dispatch lands in a follow-up commit.") .def_readwrite("ediis_diis_switch_threshold", &vibeqc::PeriodicRHFOptions::ediis_diis_switch_threshold) .def_readwrite("initial_guess", &vibeqc::PeriodicRHFOptions::initial_guess) .def_readwrite("level_shift", &vibeqc::PeriodicRHFOptions::level_shift, "Saunders-Hillier level shift (Hartree). Adds " "b · (S − ½ S D S) to F before diagonalisation, " "raising virtual orbital eigenvalues by b. " "Default 0.0 (no shift). Useful when DIIS " "oscillates on small-HOMO–LUMO-gap insulators.") .def_readwrite("level_shift_warmup_cycles", &vibeqc::PeriodicRHFOptions::level_shift_warmup_cycles, "CRYSTAL-style level-shift warm-up cycles. " "-1 auto, 0 persistent shift, positive values " "request an explicit shifted startup length.") .def_readwrite("smearing_temperature", &vibeqc::PeriodicRHFOptions::smearing_temperature, "Fermi-Dirac smearing temperature (Hartree). " "T > 0 replaces hard Aufbau with smooth " "occupations and reports the free energy " "A = E − T·S. Default 0.0 (no smearing).") .def_readwrite("quadratic_fallback_iter", &vibeqc::PeriodicRHFOptions::quadratic_fallback_iter, "Phase C1c — second-order SCF fallback " "activation iteration. When > 0, after this " "iteration the Ewald drivers switch from " "'diagonalize F' to a Newton step in MO space " "with diagonal-Hessian preconditioning. Use " "for small-gap systems where DIIS + level " "shift fail to converge. Default 0 (disabled).") .def_readwrite("quadratic_fallback_shift", &vibeqc::PeriodicRHFOptions::quadratic_fallback_shift, "Damping shift λ in the orbital-rotation " "denominator (ε_a − ε_i + λ) of the C1c " "Newton step. Default 0.1 Hartree. Larger " "values produce more conservative steps.") .def_readwrite("quadratic_fallback_max_step", &vibeqc::PeriodicRHFOptions::quadratic_fallback_max_step, "Trust-region cap on ‖κ‖ per C1c Newton " "step. Default 0.1. Larger values produce " "more aggressive rotations; smaller values " "are safer.") .def_readwrite("lattice_opts", &vibeqc::PeriodicRHFOptions::lattice_opts) .def_readwrite("dft_plus_u_sites", &vibeqc::PeriodicRHFOptions::dft_plus_u_sites, "Internal: list of _HubbardSiteCxx. Populated by " "the Python run_rhf_periodic_gamma wrapper.") .def_readwrite("dft_plus_u_ao_groups", &vibeqc::PeriodicRHFOptions::dft_plus_u_ao_groups, "Internal: parallel AO-index lists per site."); py::class_(m, "PeriodicRHFResult") .def_readonly("energy", &vibeqc::PeriodicRHFResult::energy) .def_readonly("e_electronic", &vibeqc::PeriodicRHFResult::e_electronic) .def_readonly("e_nuclear", &vibeqc::PeriodicRHFResult::e_nuclear) .def_readonly("e_dft_plus_u", &vibeqc::PeriodicRHFResult::e_dft_plus_u, "Dudarev DFT+U contribution per unit cell (Hartree). " "0 unless PeriodicRHFOptions.dft_plus_u_sites is " "non-empty.") .def_readonly("n_iter", &vibeqc::PeriodicRHFResult::n_iter) .def_readonly("converged", &vibeqc::PeriodicRHFResult::converged) .def_readonly("mo_energies", &vibeqc::PeriodicRHFResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::PeriodicRHFResult::mo_coeffs) .def_readonly("density", &vibeqc::PeriodicRHFResult::density) .def_readonly("fock", &vibeqc::PeriodicRHFResult::fock) .def_readonly("overlap", &vibeqc::PeriodicRHFResult::overlap) .def_readonly("scf_trace", &vibeqc::PeriodicRHFResult::scf_trace) .def("__repr__", [](const vibeqc::PeriodicRHFResult& r) { return std::string{"PeriodicRHFResult(energy="} + std::to_string(r.energy) + ", converged=" + (r.converged ? "True" : "False") + ", n_iter=" + std::to_string(r.n_iter) + ")"; }); m.def("run_rhf_periodic_gamma", &vibeqc::run_rhf_periodic_gamma, py::arg("system"), py::arg("basis"), py::arg("options") = vibeqc::PeriodicRHFOptions{}, py::call_guard(), "Γ-only closed-shell RHF for a periodic system (molecular-limit " "regime)."); // ----- Multi-k periodic RHF (Phase 12c) ------------------------------ m.def("build_fock_2e_real_space", &vibeqc::build_fock_2e_real_space, py::arg("basis"), py::arg("system"), py::arg("options"), py::arg("density"), py::arg("exchange_scale") = 1.0, py::arg("omega") = 0.0, py::call_guard(), "General multi-k two-electron Fock F^{2e}(g) = J(g) − " "½·exchange_scale·K(g) from a real-space density matrix. " "exchange_scale=0 skips the K computation (pure DFT). " "omega > 0 swaps the 1/r_12 kernel for erfc(ω·r_12)/r_12 " "(short-range Ewald split)."); py::class_(m, "JKLatticeMatrixSets") .def_readonly("J", &vibeqc::JKLatticeMatrixSets::J) .def_readonly("K", &vibeqc::JKLatticeMatrixSets::K); m.def("build_jk_2e_real_space", &vibeqc::build_jk_2e_real_space, py::arg("basis"), py::arg("system"), py::arg("options"), py::arg("density"), py::arg("omega") = 0.0, py::call_guard(), "General multi-k two-electron component build returning " "separate real-space J(g) and full K(g) blocks. K has no " "-1/2 prefactor; callers form J - 1/2*exchange_scale*K."); m.def("build_jk_2e_real_space_explicit", &vibeqc::build_jk_2e_real_space_explicit, py::arg("basis"), py::arg("system"), py::arg("options"), py::arg("density"), py::arg("cells"), py::arg("omega") = 0.0, py::call_guard(), "Same as build_jk_2e_real_space, but with caller-supplied " "cell list (Phase SYM3b — symmetry-reduced Fock build)."); m.def("real_space_density_from_kpoints", &vibeqc::real_space_density_from_kpoints, py::arg("C_per_k"), py::arg("n_occ_per_k"), py::arg("kmesh"), py::arg("cells"), py::call_guard(), "Fold per-k MO coefficients into a real-space density matrix P(g)."); m.def("real_space_density_from_kpoints_fractional", &vibeqc::real_space_density_from_kpoints_fractional, py::arg("C_per_k"), py::arg("occ_per_k"), py::arg("kmesh"), py::arg("cells"), py::call_guard(), "Fractional-occupation variant of real_space_density_from_kpoints. " "``occ_per_k`` is a list of double arrays — one occupation per MO " "per k-point, typically in [0, 2] for closed-shell RHF. Used by " "Fermi-Dirac-smearing-driven SCF (Phase C1b)."); py::class_(m, "PeriodicSCFOptions") .def(py::init<>()) .def_readwrite("max_iter", &vibeqc::PeriodicSCFOptions::max_iter) .def_readwrite("conv_tol_energy", &vibeqc::PeriodicSCFOptions::conv_tol_energy) .def_readwrite("conv_tol_grad", &vibeqc::PeriodicSCFOptions::conv_tol_grad) .def_readwrite("damping", &vibeqc::PeriodicSCFOptions::damping) .def_readwrite("dynamic_damping", &vibeqc::PeriodicSCFOptions::dynamic_damping, "See RHFOptions.dynamic_damping. Default False.") .def_readwrite("dynamic_damping_min", &vibeqc::PeriodicSCFOptions::dynamic_damping_min) .def_readwrite("dynamic_damping_max", &vibeqc::PeriodicSCFOptions::dynamic_damping_max) .def_readwrite("fock_mixing", &vibeqc::PeriodicSCFOptions::fock_mixing, "Fock matrix mixing fraction. See " "RHFOptions.fock_mixing.") .def_readwrite("use_diis", &vibeqc::PeriodicSCFOptions::use_diis) .def_readwrite("diis_start_iter", &vibeqc::PeriodicSCFOptions::diis_start_iter) .def_readwrite("diis_subspace_size", &vibeqc::PeriodicSCFOptions::diis_subspace_size) .def_readwrite("scf_accelerator", &vibeqc::PeriodicSCFOptions::scf_accelerator, "SCFAccelerator family. Multi-k dispatch lands in " "a follow-up commit; field is exposed now for " "API uniformity.") .def_readwrite("ediis_diis_switch_threshold", &vibeqc::PeriodicSCFOptions::ediis_diis_switch_threshold) .def_readwrite("initial_guess", &vibeqc::PeriodicSCFOptions::initial_guess) .def_readwrite("level_shift", &vibeqc::PeriodicSCFOptions::level_shift, "Saunders-Hillier level shift (Hartree). Default 0.0.") .def_readwrite("level_shift_warmup_cycles", &vibeqc::PeriodicSCFOptions::level_shift_warmup_cycles, "CRYSTAL-style level-shift warm-up cycles. " "-1 auto, 0 persistent shift, positive values " "request an explicit shifted startup length.") .def_readwrite("smearing_temperature", &vibeqc::PeriodicSCFOptions::smearing_temperature, "Fermi-Dirac smearing temperature (Hartree). Default 0.0.") .def_readwrite("quadratic_fallback_iter", &vibeqc::PeriodicSCFOptions::quadratic_fallback_iter, "C1c second-order SCF fallback activation iter. " "Default 0 (disabled).") .def_readwrite("quadratic_fallback_shift", &vibeqc::PeriodicSCFOptions::quadratic_fallback_shift, "C1c shift λ in (ε_a − ε_i + λ). Default 0.1.") .def_readwrite("quadratic_fallback_max_step", &vibeqc::PeriodicSCFOptions::quadratic_fallback_max_step, "C1c trust-region cap on ‖κ‖ per Newton step. " "Default 0.1.") .def_readwrite("lattice_opts", &vibeqc::PeriodicSCFOptions::lattice_opts); m.def("run_rhf_periodic", &vibeqc::run_rhf_periodic, py::arg("system"), py::arg("basis"), py::arg("kmesh"), py::arg("options") = vibeqc::PeriodicSCFOptions{}, py::call_guard(), "Multi-k closed-shell RHF for a periodic system — real-space " "density matrix, general three-lattice-index Fock build."); // ----- Multi-k periodic KS-DFT (Phase 12d) --------------------------- py::class_(m, "PeriodicXCContribution") .def_readonly("V_xc", &vibeqc::PeriodicXCContribution::V_xc) .def_readonly("e_xc", &vibeqc::PeriodicXCContribution::e_xc); m.def("build_xc_periodic", &vibeqc::build_xc_periodic, py::arg("basis"), py::arg("system"), py::arg("grid"), py::arg("functional"), py::arg("density"), py::arg("options"), py::call_guard(), "Periodic XC Fock contribution V_xc(g) and E_xc per unit cell."); // Open-shell periodic XC — companion of build_xc_periodic for UKS. py::class_(m, "PeriodicUKSXCContribution") .def_readonly("V_alpha", &vibeqc::PeriodicUKSXCContribution::V_alpha, "V_xc^α(g) as a LatticeMatrixSet — per-cell α-spin " "XC potential matrices.") .def_readonly("V_beta", &vibeqc::PeriodicUKSXCContribution::V_beta, "V_xc^β(g) as a LatticeMatrixSet — per-cell β-spin " "XC potential matrices.") .def_readonly("e_xc", &vibeqc::PeriodicUKSXCContribution::e_xc, "Total E_xc per unit cell (Hartree)."); m.def("build_xc_periodic_uks", &vibeqc::build_xc_periodic_uks, py::arg("basis"), py::arg("system"), py::arg("grid"), py::arg("functional"), py::arg("density_alpha"), py::arg("density_beta"), py::arg("options"), py::call_guard(), "Open-shell (UKS) periodic XC: takes per-spin LatticeMatrixSet " "densities P_α(g), P_β(g) and returns per-spin V_xc matrices + " "total E_xc per unit cell. The Functional must be constructed " "with spin=2 (eval_polarised path); spin=1 will throw."); py::class_(m, "PeriodicKSOptions") .def(py::init<>()) .def_readwrite("functional", &vibeqc::PeriodicKSOptions::functional) .def_readwrite("grid", &vibeqc::PeriodicKSOptions::grid) .def_readwrite("max_iter", &vibeqc::PeriodicKSOptions::max_iter) .def_readwrite("conv_tol_energy", &vibeqc::PeriodicKSOptions::conv_tol_energy) .def_readwrite("conv_tol_grad", &vibeqc::PeriodicKSOptions::conv_tol_grad) .def_readwrite("damping", &vibeqc::PeriodicKSOptions::damping) .def_readwrite("dynamic_damping", &vibeqc::PeriodicKSOptions::dynamic_damping, "See RHFOptions.dynamic_damping. Default False.") .def_readwrite("dynamic_damping_min", &vibeqc::PeriodicKSOptions::dynamic_damping_min) .def_readwrite("dynamic_damping_max", &vibeqc::PeriodicKSOptions::dynamic_damping_max) .def_readwrite("fock_mixing", &vibeqc::PeriodicKSOptions::fock_mixing, "Kohn-Sham matrix mixing fraction. See " "RHFOptions.fock_mixing.") .def_readwrite("use_diis", &vibeqc::PeriodicKSOptions::use_diis) .def_readwrite("diis_start_iter", &vibeqc::PeriodicKSOptions::diis_start_iter) .def_readwrite("scf_accelerator", &vibeqc::PeriodicKSOptions::scf_accelerator, "SCFAccelerator family. Multi-k dispatch lands in " "a follow-up commit.") .def_readwrite("ediis_diis_switch_threshold", &vibeqc::PeriodicKSOptions::ediis_diis_switch_threshold) .def_readwrite("diis_subspace_size", &vibeqc::PeriodicKSOptions::diis_subspace_size) .def_readwrite("initial_guess", &vibeqc::PeriodicKSOptions::initial_guess) .def_readwrite("level_shift", &vibeqc::PeriodicKSOptions::level_shift, "Saunders-Hillier level shift (Hartree). Default 0.0.") .def_readwrite("level_shift_warmup_cycles", &vibeqc::PeriodicKSOptions::level_shift_warmup_cycles, "CRYSTAL-style level-shift warm-up cycles. " "-1 auto, 0 persistent shift, positive values " "request an explicit shifted startup length.") .def_readwrite("smearing_temperature", &vibeqc::PeriodicKSOptions::smearing_temperature, "Fermi-Dirac smearing temperature (Hartree). Default 0.0.") .def_readwrite("quadratic_fallback_iter", &vibeqc::PeriodicKSOptions::quadratic_fallback_iter, "C1c second-order SCF fallback activation iter. " "Default 0 (disabled).") .def_readwrite("quadratic_fallback_shift", &vibeqc::PeriodicKSOptions::quadratic_fallback_shift, "C1c shift λ in (ε_a − ε_i + λ). Default 0.1.") .def_readwrite("quadratic_fallback_max_step", &vibeqc::PeriodicKSOptions::quadratic_fallback_max_step, "C1c trust-region cap on ‖κ‖ per Newton step. " "Default 0.1.") .def_readwrite("lattice_opts", &vibeqc::PeriodicKSOptions::lattice_opts) .def_readwrite("use_periodic_becke", &vibeqc::PeriodicKSOptions::use_periodic_becke, "When true (default, since v0.9.x), the DFT " "integration grid uses the periodic Becke " "partition (extends the partition denominator " "over image atoms within becke_image_radius_bohr). " "Reduces exactly to the molecular partition in the " "molecular-limit regime, so it is safe to leave on " "for any periodic calculation. Set false only to " "reproduce v0.8.x numerics — silently wrong on " "tight cells.") .def_readwrite("becke_image_radius_bohr", &vibeqc::PeriodicKSOptions::becke_image_radius_bohr, "Radius (bohr) within which image atoms are " "included in the periodic Becke partition. Only " "consulted when use_periodic_becke = true. " "Default 10 bohr.") .def_readwrite("dft_plus_u_sites", &vibeqc::PeriodicKSOptions::dft_plus_u_sites, "Internal: list of _HubbardSiteCxx for the Γ-only " "periodic RKS +U path (Increment 4b). Populated by " "the Python run_rks_periodic wrapper; multi-k +U " "support is queued for Increment 4c.") .def_readwrite("dft_plus_u_ao_groups", &vibeqc::PeriodicKSOptions::dft_plus_u_ao_groups, "Internal: parallel AO-index lists per site."); py::class_(m, "PeriodicKSResult") .def_readonly("energy", &vibeqc::PeriodicKSResult::energy) .def_readonly("e_electronic", &vibeqc::PeriodicKSResult::e_electronic) .def_readonly("e_coulomb", &vibeqc::PeriodicKSResult::e_coulomb) .def_readonly("e_hf_exchange", &vibeqc::PeriodicKSResult::e_hf_exchange) .def_readonly("e_xc", &vibeqc::PeriodicKSResult::e_xc) .def_readonly("e_nuclear", &vibeqc::PeriodicKSResult::e_nuclear) .def_readonly("e_dft_plus_u", &vibeqc::PeriodicKSResult::e_dft_plus_u, "Dudarev DFT+U contribution per unit cell " "(Hartree). 0 unless dft_plus_u_sites is set.") .def_readonly("n_iter", &vibeqc::PeriodicKSResult::n_iter) .def_readonly("converged", &vibeqc::PeriodicKSResult::converged) .def_readonly("mo_energies", &vibeqc::PeriodicKSResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::PeriodicKSResult::mo_coeffs) .def_readonly("density", &vibeqc::PeriodicKSResult::density) .def_readonly("fock", &vibeqc::PeriodicKSResult::fock) .def_readonly("overlap", &vibeqc::PeriodicKSResult::overlap) .def_readonly("scf_trace", &vibeqc::PeriodicKSResult::scf_trace) .def_readonly("functional", &vibeqc::PeriodicKSResult::functional) .def("__repr__", [](const vibeqc::PeriodicKSResult& r) { return std::string{"PeriodicKSResult(energy="} + std::to_string(r.energy) + ", functional=" + r.functional + ", converged=" + (r.converged ? "True" : "False") + ")"; }); m.def("run_rks_periodic", &vibeqc::run_rks_periodic, py::arg("system"), py::arg("basis"), py::arg("kmesh"), py::arg("options") = vibeqc::PeriodicKSOptions{}, py::call_guard(), "Multi-k closed-shell Kohn-Sham DFT for a periodic system."); // ---- DFT+U (Dudarev) — C++ kernel surface -------------------------- // // ``HubbardSiteCxx`` is the internal POD used by the C++ SCF Fock // builder. Users construct ``vibeqc.HubbardSite`` (Python frozen // dataclass, eV inputs); the Python wrapper converts to this struct // (eV → Hartree, ao_groups precomputed) before calling into C++. py::class_(m, "_HubbardSiteCxx", "Internal POD for the C++ DFT+U kernel. Mirrors " "vibeqc.HubbardSite but stores U_eff in Hartree (the unit " "conversion happens at the Python boundary). Not intended for " "direct user construction — use vibeqc.HubbardSite.") .def(py::init<>()) .def(py::init([](int atom_index, int l, double U_eff_au) { return vibeqc::HubbardSiteCxx{atom_index, l, U_eff_au}; }), py::arg("atom_index"), py::arg("l"), py::arg("U_eff_au")) .def_readwrite("atom_index", &vibeqc::HubbardSiteCxx::atom_index) .def_readwrite("l", &vibeqc::HubbardSiteCxx::l) .def_readwrite("U_eff_au", &vibeqc::HubbardSiteCxx::U_eff_au); m.def("_compute_dft_plus_u_cxx", [](const std::vector& sites, const std::vector>& ao_groups, const Eigen::MatrixXd& P, const Eigen::MatrixXd& S) { auto r = vibeqc::compute_dft_plus_u(sites, ao_groups, P, S); return py::make_tuple(r.energy, std::move(r.V)); }, py::arg("sites"), py::arg("ao_groups"), py::arg("P"), py::arg("S"), "Internal: per-spin Dudarev +U energy + AO Fock contribution. " "Returns (E_U_sigma, V_U). Spin convention is the caller's " "responsibility (pass P_sigma for the per-spin formula)."); // ---- Dispersion (D3-BJ) -------------------------------------------- py::class_(m, "D3BJParams", "Four damping parameters of the D3 Becke-Johnson correction.") .def(py::init<>()) .def(py::init([](double s6, double s8, double a1, double a2) { return vibeqc::D3BJParams{s6, s8, a1, a2}; }), py::arg("s6") = 1.0, py::arg("s8") = 0.0, py::arg("a1") = 0.0, py::arg("a2") = 0.0) .def_readwrite("s6", &vibeqc::D3BJParams::s6) .def_readwrite("s8", &vibeqc::D3BJParams::s8) .def_readwrite("a1", &vibeqc::D3BJParams::a1) .def_readwrite("a2", &vibeqc::D3BJParams::a2) .def("__repr__", [](const vibeqc::D3BJParams& p) { return "D3BJParams(s6=" + std::to_string(p.s6) + ", s8=" + std::to_string(p.s8) + ", a1=" + std::to_string(p.a1) + ", a2=" + std::to_string(p.a2) + ")"; }); py::class_(m, "DispersionResult", "Dispersion energy and (optional) atomic gradient.") .def_readonly("energy", &vibeqc::DispersionResult::energy) .def_readonly("gradient", &vibeqc::DispersionResult::gradient); m.def("d3bj_params_for", [](const std::string& functional) -> py::object { auto p = vibeqc::d3bj_params_for(functional); if (!p) return py::none(); return py::cast(*p); }, py::arg("functional"), "Look up D3-BJ damping parameters for a DFT functional " "(case-insensitive). Returns None if the functional is not " "in the shipped parameter table."); m.def("compute_d3bj", &vibeqc::compute_d3bj, py::arg("molecule"), py::arg("params"), py::arg("with_gradient") = false, py::call_guard(), "Pairwise D3(BJ) dispersion energy (and optionally gradient) " "for a molecule given the functional's damping parameters."); m.def("d3_coordination_numbers", [](const vibeqc::Molecule& mol) { return vibeqc::coordination_numbers(mol).cn; }, py::arg("molecule"), py::call_guard(), "Smooth D3 coordination numbers (Grimme counting function) for " "every atom in the molecule."); m.def("d3_r2r4", &vibeqc::r2r4, py::arg("Z"), "r2r4 atomic ratio sqrt(/) used in the D3 C8 scaling."); m.def("d3_rcov", &vibeqc::rcov, py::arg("Z"), "D3 covalent radius in bohr (Pyykko single-bond radius scaled " "by k2 = 4/3)."); // ---- Chai-Head-Gordon dispersion — the ωB97X-D "-D" term ---------- m.def("compute_chg_dispersion", &vibeqc::compute_chg_dispersion, py::arg("molecule"), py::arg("with_gradient") = false, py::call_guard(), "Chai-Head-Gordon ('CHG' / DFT-D2-type) empirical dispersion " "energy — the intrinsic '-D' correction of ωB97X-D. " "E = -s6 Σ_{A(m, "EEQResult", "Result of an EEQ charge calculation. ``charges`` is an " "ndarray of n atomic partial charges (electrons); " "``chemical_potential`` is the Lagrange multiplier λ of the " "charge-conservation constraint (== the system chemical " "potential).") .def_readonly("charges", &vibeqc::EEQResult::charges) .def_readonly("chemical_potential", &vibeqc::EEQResult::chemical_potential) .def("__repr__", [](const vibeqc::EEQResult& r) { return std::string{"EEQResult(n_atoms="} + std::to_string(r.charges.size()) + ", chemical_potential=" + std::to_string(r.chemical_potential) + ")"; }); py::class_(m, "EEQOptions", "EEQ-CN counting + solver knobs. Defaults match " "multicharge::new_eeq2019_model (cn_cutoff = 25.0 bohr, " "cn_exp = 7.5, cn_max = 8.0).") .def(py::init<>()) .def_readwrite("cn_cutoff", &vibeqc::EEQOptions::cn_cutoff) .def_readwrite("cn_exp", &vibeqc::EEQOptions::cn_exp) .def_readwrite("cn_max", &vibeqc::EEQOptions::cn_max); m.def("eeq_charges", &vibeqc::eeq_charges, py::arg("molecule"), py::arg("total_charge") = 0.0, py::arg("options") = vibeqc::EEQOptions{}, py::call_guard(), "Compute electronegativity-equilibration (EEQ, Caldeweyher " "2019) atomic partial charges for a molecule. Returns an " "EEQResult. Apache-2.0 port of Grimme group's " "multicharge::eeq2019 model. The algorithmically " "distinctive piece of D4 dispersion over D3 — full D4 " "dispersion energy in vibe-qc's native backend will land " "in a follow-on (Phase D4b); for production D4 numbers, " "use vibeqc.compute_d4 (which wraps the optional dftd4 " "Python package)."); m.def("eeq_coordination_numbers", &vibeqc::eeq_coordination_numbers, py::arg("molecule"), py::arg("options") = vibeqc::EEQOptions{}, py::call_guard(), "EEQ-specific coordination numbers (erf counting function, " "kcn = 7.5, bare Pyykko covalent radii). Different from " "vibeqc.d3_coordination_numbers, which uses the D3 " "counting function + 4/3-scaled radii."); // ---- Reciprocal-space periodic Poisson (Phase 12e-c-3a) ----------- // // Accept + return a flat numpy array rather than a custom // ScalarField3D pybind11 class; the C++-side struct is just a // shape + std::vector pair, not worth exposing. Users pass a // (nx, ny, nz) numpy array in and get a same-shaped one out. auto numpy_to_field = [](py::array_t arr) -> vibeqc::ScalarField3D { if (arr.ndim() != 3) { throw std::invalid_argument( "solve_poisson: rho must be a 3D numpy array"); } vibeqc::ScalarField3D out; out.resize(static_cast(arr.shape(0)), static_cast(arr.shape(1)), static_cast(arr.shape(2))); const double* src = arr.data(); std::copy(src, src + out.data.size(), out.data.begin()); return out; }; auto field_to_numpy = [](const vibeqc::ScalarField3D& f) -> py::array_t { py::array_t out({f.nx, f.ny, f.nz}); auto buf = out.mutable_unchecked<3>(); for (std::size_t i = 0; i < f.nx; ++i) for (std::size_t j = 0; j < f.ny; ++j) for (std::size_t k = 0; k < f.nz; ++k) buf(i, j, k) = f(i, j, k); return out; }; m.def("solve_poisson_erf_screened", [numpy_to_field, field_to_numpy]( py::array_t rho, const Eigen::Matrix3d& lattice, double omega) { auto rho_f = numpy_to_field(rho); auto V = vibeqc::solve_poisson_erf_screened(rho_f, lattice, omega); return field_to_numpy(V); }, py::arg("rho"), py::arg("lattice"), py::arg("omega"), "Periodic Poisson solver with the erf-screened Coulomb kernel " "K(r) = erf(omega r) / r. Uses FFTW3 to apply Ṽ(G) = (4π/G²)" " exp(-G²/4ω²) ρ̃(G). Returned V has zero mean (G=0 component " "dropped)."); m.def("solve_poisson_coulomb", [numpy_to_field, field_to_numpy]( py::array_t rho, const Eigen::Matrix3d& lattice) { auto rho_f = numpy_to_field(rho); auto V = vibeqc::solve_poisson_coulomb(rho_f, lattice); return field_to_numpy(V); }, py::arg("rho"), py::arg("lattice"), "Periodic Poisson solver with the unscreened Coulomb kernel " "K(r) = 1/r. Ṽ(G) = (4π/G²) ρ̃(G). Returned V has zero mean."); m.def("hartree_energy_on_grid", [numpy_to_field]( py::array_t rho, py::array_t V, double cell_volume_bohr3) { auto rho_f = numpy_to_field(rho); auto V_f = numpy_to_field(V); return vibeqc::hartree_energy(rho_f, V_f, cell_volume_bohr3); }, py::arg("rho"), py::arg("V"), py::arg("cell_volume_bohr3"), "E_H = (1/2) integral rho(r) V(r) dr computed on the grid. " "Test utility — most callers compute this themselves."); // ================================================================ // Semiempirical DFTB0 (Stage 1) // ================================================================ py::module_ m_semi = m.def_submodule("semiempirical", "Semiempirical tight-binding methods (DFTB0, SCC-DFTB, ...)"); // --- SemiempiricalBasis --- py::class_( m_semi, "SemiempiricalBasis", "Minimal valence STO-NG basis builder.") .def_static("build", py::overload_cast( &vibeqc::semiempirical::SemiempiricalBasis::build), py::arg("mol"), py::arg("params"), py::arg("n_primitives") = 6, "Build basis from CoreParameterSet."); py::class_( m_semi, "SemiempiricalParameters", "DFTB0 / SCC-DFTB parameter set: on-site energies, STO exponents, " "repulsive pair potentials. Use SemiempiricalParameters.dftb0_default() " "for the built-in H/C/N/O table.") .def(py::init<>()) .def_static("dftb0_production",&vibeqc::semiempirical::SemiempiricalParameters::dftb0_production,"Production DFTB parameter set with spline-fitted repulsives.").def_static("dftb0_default", &vibeqc::semiempirical::SemiempiricalParameters::dftb0_default, "Return the built-in DFTB0 parameter set for H, C, N, O.") .def("add_element", [](vibeqc::semiempirical::SemiempiricalParameters& p, int Z, const std::vector& on_site, const std::vector& zeta, double U, int nval) { vibeqc::semiempirical::ElementData e; e.Z = Z; e.on_site = on_site; e.zeta = zeta; e.hubbard_u = U; e.valence_electrons = nval; p.add_element(e); }, py::arg("Z"), py::arg("on_site"), py::arg("zeta"), py::arg("hubbard_u"), py::arg("valence_electrons"), "Add an element to the parameter set.") .def("has_element", &vibeqc::semiempirical::SemiempiricalParameters::has_element, py::arg("Z"), "True if element Z is in the parameter set.") .def("on_site_energy", &vibeqc::semiempirical::SemiempiricalParameters::on_site_energy, py::arg("Z"), py::arg("l"), "On-site energy for element Z, ang.mom. l (Hartree).") .def("sto_exponent", &vibeqc::semiempirical::SemiempiricalParameters::sto_exponent, py::arg("Z"), py::arg("l"), "STO exponent for element Z, ang.mom. l (a0^-1).") .def("average_on_site", &vibeqc::semiempirical::SemiempiricalParameters::average_on_site, py::arg("Z"), "Average on-site energy for element Z (Hartree).") .def("hubbard_u", &vibeqc::semiempirical::SemiempiricalParameters::hubbard_u, py::arg("Z"), "Hubbard U parameter for element Z (Hartree).") .def("valence_electrons", &vibeqc::semiempirical::SemiempiricalParameters::valence_electrons, py::arg("Z"), "Neutral-atom valence electron count for element Z.") .def("gamma_onsite", &vibeqc::semiempirical::SemiempiricalParameters::gamma_onsite, py::arg("Z"), "On-site gamma = U_A (Hartree).") .def("set_repulsive_pair_analytic", [](vibeqc::semiempirical::SemiempiricalParameters& p, int Z1, int Z2, double A, double B) { vibeqc::semiempirical::RepulsivePairV2 rp; rp.A = A; rp.B = B; p.set_repulsive_pair(Z1, Z2, rp); }, py::arg("Z1"), py::arg("Z2"), py::arg("A"), py::arg("B") = 0.0, "Set analytic repulsive pair: A/R^12 if B=0, else A*exp(-B*R).") .def("set_repulsive_pair_spline", [](vibeqc::semiempirical::SemiempiricalParameters& p, int Z1, int Z2, const std::vector& R, const std::vector& V) { vibeqc::semiempirical::RepulsivePairV2 rp; rp.spline.build(R, V); p.set_repulsive_pair(Z1, Z2, rp); }, py::arg("Z1"), py::arg("Z2"), py::arg("R"), py::arg("V"), "Set spline-fitted repulsive pair from knot arrays (R, V). " "R and V must be the same length, R strictly increasing.") .def("has_repulsive_pair", [](const vibeqc::semiempirical::SemiempiricalParameters& p, int Z1, int Z2) -> bool { return p.repulsive_pair(Z1, Z2) != nullptr; }, py::arg("Z1"), py::arg("Z2"), "True if an explicit repulsive pair is set for (Z1,Z2).") .def("repulsive_energy", &vibeqc::semiempirical::SemiempiricalParameters::repulsive_energy, py::arg("Z1"), py::arg("Z2"), py::arg("R"), "Evaluate V_rep(R) for pair (Z1, Z2) at distance R (bohr).") .def_property("kappa", &vibeqc::semiempirical::SemiempiricalParameters::kappa, &vibeqc::semiempirical::SemiempiricalParameters::set_kappa, "Wolfsberg-Helmholtz scaling constant (default 1.75)."); py::class_( m_semi, "DFTB0Result", "Result of a non-self-consistent DFTB0 energy calculation.") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::DFTB0Result::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::DFTB0Result::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::DFTB0Result::e_repulsive) .def_readonly("mo_energies", &vibeqc::semiempirical::DFTB0Result::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::DFTB0Result::mo_coeffs) .def_readonly("density", &vibeqc::semiempirical::DFTB0Result::density) .def_readonly("overlap", &vibeqc::semiempirical::DFTB0Result::overlap) .def_readonly("hamiltonian", &vibeqc::semiempirical::DFTB0Result::hamiltonian) .def_readonly("n_basis", &vibeqc::semiempirical::DFTB0Result::n_basis) .def_readonly("n_occ", &vibeqc::semiempirical::DFTB0Result::n_occ); m_semi.def("run_dftb0", &vibeqc::semiempirical::run_dftb0, py::arg("mol"), py::arg("params"), py::call_guard(), "Run a non-self-consistent DFTB0 energy calculation."); m_semi.def("compute_dftb0_gradient", &vibeqc::semiempirical::compute_dftb0_gradient, py::arg("mol"), py::arg("result"), py::arg("params"), py::call_guard(), "Compute the analytic nuclear gradient for a DFTB0 result. " "Returns (n_atoms, 3) matrix in Hartree/bohr."); m_semi.def("compute_scc_dftb_gradient_response", &vibeqc::semiempirical::compute_scc_dftb_gradient_response, py::arg("mol"), py::arg("result"), py::arg("params"), py::call_guard(), "SCC-DFTB gradient with CP-SCC charge-response correction."); m_semi.def("compute_scc_dftb_gradient", &vibeqc::semiempirical::compute_scc_dftb_gradient, py::arg("mol"), py::arg("result"), py::arg("params"), py::call_guard(), "Compute the analytic nuclear gradient for an SCC-DFTB result."); m_semi.def("compute_uscc_dftb_gradient", &vibeqc::semiempirical::compute_uscc_dftb_gradient, py::arg("mol"), py::arg("result"), py::arg("params"), py::call_guard(), "Compute the analytic nuclear gradient for an unrestricted SCC-DFTB result."); // --- Periodic gradients --- m_semi.def("compute_periodic_dftb0_gradient", &vibeqc::semiempirical::compute_periodic_dftb0_gradient, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic DFTB0 atomic gradient at Gamma-point."); m_semi.def("compute_periodic_scc_dftb_gradient", &vibeqc::semiempirical::compute_periodic_scc_dftb_gradient, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic SCC-DFTB atomic gradient at Gamma-point."); m_semi.def("compute_periodic_udftb0_gradient", &vibeqc::semiempirical::compute_periodic_udftb0_gradient, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic UDFTB0 atomic gradient at Gamma-point."); m_semi.def("compute_periodic_uscc_dftb_gradient", &vibeqc::semiempirical::compute_periodic_uscc_dftb_gradient, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic USCC-DFTB atomic gradient at Gamma-point."); // --- Periodic stress --- m_semi.def("compute_periodic_dftb0_stress", &vibeqc::semiempirical::compute_periodic_dftb0_stress, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic DFTB0 analytic stress tensor at Gamma-point (3x3, Ha/bohr^3)."); m_semi.def("compute_periodic_scc_dftb_stress", &vibeqc::semiempirical::compute_periodic_scc_dftb_stress, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic SCC-DFTB analytic stress tensor at Gamma-point (3x3, Ha/bohr^3)."); m_semi.def("dftb0_repulsive_gradient", &vibeqc::semiempirical::dftb0_repulsive_gradient, py::arg("mol"), py::arg("params"), "Repulsive-energy gradient only (for testing). " "Returns (n_atoms, 3) matrix in Hartree/bohr."); // --- SCC-DFTB (Stage 3) --- py::class_( m_semi, "SCCOptions", "SCC-DFTB SCF options.") .def(py::init<>()) .def_readwrite("max_iter", &vibeqc::semiempirical::SCCOptions::max_iter) .def_readwrite("conv_tol_charge", &vibeqc::semiempirical::SCCOptions::conv_tol_charge) .def_readwrite("charge_mixing", &vibeqc::semiempirical::SCCOptions::charge_mixing); py::class_( m_semi, "SCCDFTBResult", "Result of an SCC-DFTB calculation.") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::SCCDFTBResult::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::SCCDFTBResult::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::SCCDFTBResult::e_repulsive) .def_readonly("e_scc", &vibeqc::semiempirical::SCCDFTBResult::e_scc) .def_readonly("mo_energies", &vibeqc::semiempirical::SCCDFTBResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::SCCDFTBResult::mo_coeffs) .def_readonly("density", &vibeqc::semiempirical::SCCDFTBResult::density) .def_readonly("overlap", &vibeqc::semiempirical::SCCDFTBResult::overlap) .def_readonly("hamiltonian", &vibeqc::semiempirical::SCCDFTBResult::hamiltonian) .def_readonly("charges", &vibeqc::semiempirical::SCCDFTBResult::charges) .def_readonly("n_basis", &vibeqc::semiempirical::SCCDFTBResult::n_basis) .def_readonly("n_occ", &vibeqc::semiempirical::SCCDFTBResult::n_occ) .def_readonly("n_iter", &vibeqc::semiempirical::SCCDFTBResult::n_iter) .def_readonly("converged", &vibeqc::semiempirical::SCCDFTBResult::converged); m_semi.def("run_scc_dftb", &vibeqc::semiempirical::run_scc_dftb, py::arg("mol"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::SCCOptions{}, py::call_guard(), "Run an SCC-DFTB calculation with charge self-consistency."); // Add Hubbard U accessor to SemiempiricalParameters (already defined above, // just need to add the binding). We patch it in via a follow-up .def call. // --- Periodic DFTB0 (Stage 4) --- py::class_( m_semi, "PeriodicDFTB0Options", "Options for periodic Gamma-point DFTB0.") .def(py::init<>()) .def_readwrite("cutoff_bohr", &vibeqc::semiempirical::PeriodicDFTB0Options::cutoff_bohr) .def_readwrite("gamma_only_0", &vibeqc::semiempirical::PeriodicDFTB0Options::gamma_only_0); py::class_( m_semi, "PeriodicDFTB0Result", "Result of a periodic Gamma-point DFTB0 calculation.") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::PeriodicDFTB0Result::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::PeriodicDFTB0Result::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::PeriodicDFTB0Result::e_repulsive) .def_readonly("mo_energies", &vibeqc::semiempirical::PeriodicDFTB0Result::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::PeriodicDFTB0Result::mo_coeffs) .def_readonly("density", &vibeqc::semiempirical::PeriodicDFTB0Result::density) .def_readonly("overlap_gamma", &vibeqc::semiempirical::PeriodicDFTB0Result::overlap_gamma) .def_readonly("hamiltonian_gamma", &vibeqc::semiempirical::PeriodicDFTB0Result::hamiltonian_gamma) .def_readonly("n_basis", &vibeqc::semiempirical::PeriodicDFTB0Result::n_basis) .def_readonly("n_occ", &vibeqc::semiempirical::PeriodicDFTB0Result::n_occ) .def_readonly("n_cells", &vibeqc::semiempirical::PeriodicDFTB0Result::n_cells); m_semi.def("run_dftb0_gamma", &vibeqc::semiempirical::run_dftb0_gamma, py::arg("system"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::PeriodicDFTB0Options{}, py::call_guard(), "Run periodic Gamma-point DFTB0 energy calculation."); // --- Periodic UDFTB0 --- py::class_( m_semi, "PeriodicUDFTB0Result", "Result of an unrestricted periodic Gamma-point DFTB0 calculation.") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::PeriodicUDFTB0Result::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::PeriodicUDFTB0Result::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::PeriodicUDFTB0Result::e_repulsive) .def_readonly("mo_energies", &vibeqc::semiempirical::PeriodicUDFTB0Result::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::PeriodicUDFTB0Result::mo_coeffs) .def_readonly("density_alpha", &vibeqc::semiempirical::PeriodicUDFTB0Result::density_alpha) .def_readonly("density_beta", &vibeqc::semiempirical::PeriodicUDFTB0Result::density_beta) .def_readonly("overlap_gamma", &vibeqc::semiempirical::PeriodicUDFTB0Result::overlap_gamma) .def_readonly("hamiltonian_gamma", &vibeqc::semiempirical::PeriodicUDFTB0Result::hamiltonian_gamma) .def_readonly("n_basis", &vibeqc::semiempirical::PeriodicUDFTB0Result::n_basis) .def_readonly("n_alpha", &vibeqc::semiempirical::PeriodicUDFTB0Result::n_alpha) .def_readonly("n_beta", &vibeqc::semiempirical::PeriodicUDFTB0Result::n_beta) .def_readonly("n_cells", &vibeqc::semiempirical::PeriodicUDFTB0Result::n_cells); m_semi.def("run_udftb0_gamma", &vibeqc::semiempirical::run_udftb0_gamma, py::arg("system"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::PeriodicDFTB0Options{}, py::call_guard(), "Run unrestricted periodic Gamma-point DFTB0 calculation."); // --- Periodic SCC-DFTB (Stage 5) --- py::class_( m_semi, "PeriodicSCCOptions") .def(py::init<>()) .def_readwrite("cutoff_bohr", &vibeqc::semiempirical::PeriodicSCCOptions::cutoff_bohr) .def_readwrite("max_iter", &vibeqc::semiempirical::PeriodicSCCOptions::max_iter) .def_readwrite("conv_tol_charge", &vibeqc::semiempirical::PeriodicSCCOptions::conv_tol_charge) .def_readwrite("charge_mixing", &vibeqc::semiempirical::PeriodicSCCOptions::charge_mixing); py::class_( m_semi, "PeriodicSCCDFTBResult") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::PeriodicSCCDFTBResult::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::PeriodicSCCDFTBResult::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::PeriodicSCCDFTBResult::e_repulsive) .def_readonly("e_scc", &vibeqc::semiempirical::PeriodicSCCDFTBResult::e_scc) .def_readonly("mo_energies", &vibeqc::semiempirical::PeriodicSCCDFTBResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::PeriodicSCCDFTBResult::mo_coeffs) .def_readonly("density", &vibeqc::semiempirical::PeriodicSCCDFTBResult::density) .def_readonly("overlap_gamma", &vibeqc::semiempirical::PeriodicSCCDFTBResult::overlap_gamma) .def_readonly("hamiltonian_gamma", &vibeqc::semiempirical::PeriodicSCCDFTBResult::hamiltonian_gamma) .def_readonly("charges", &vibeqc::semiempirical::PeriodicSCCDFTBResult::charges) .def_readonly("n_basis", &vibeqc::semiempirical::PeriodicSCCDFTBResult::n_basis) .def_readonly("n_occ", &vibeqc::semiempirical::PeriodicSCCDFTBResult::n_occ) .def_readonly("n_cells", &vibeqc::semiempirical::PeriodicSCCDFTBResult::n_cells) .def_readonly("n_iter", &vibeqc::semiempirical::PeriodicSCCDFTBResult::n_iter) .def_readonly("converged", &vibeqc::semiempirical::PeriodicSCCDFTBResult::converged); m_semi.def("run_scc_dftb_gamma", &vibeqc::semiempirical::run_scc_dftb_gamma, py::arg("system"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::PeriodicSCCOptions{}, py::call_guard(), "Run periodic SCC-DFTB Gamma-point calculation."); // --- Periodic USCC-DFTB --- py::class_( m_semi, "PeriodicUSCCDFTBResult") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::e_repulsive) .def_readonly("e_scc", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::e_scc) .def_readonly("mo_energies", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::mo_coeffs) .def_readonly("density_alpha", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::density_alpha) .def_readonly("density_beta", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::density_beta) .def_readonly("overlap_gamma", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::overlap_gamma) .def_readonly("hamiltonian_gamma", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::hamiltonian_gamma) .def_readonly("charges", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::charges) .def_readonly("n_basis", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::n_basis) .def_readonly("n_alpha", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::n_alpha) .def_readonly("n_beta", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::n_beta) .def_readonly("n_cells", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::n_cells) .def_readonly("n_iter", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::n_iter) .def_readonly("converged", &vibeqc::semiempirical::PeriodicUSCCDFTBResult::converged); m_semi.def("run_uscc_dftb_gamma", &vibeqc::semiempirical::run_uscc_dftb_gamma, py::arg("system"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::PeriodicSCCOptions{}, py::call_guard(), "Run unrestricted periodic SCC-DFTB Gamma-point calculation."); // --- Unrestricted DFTB0 --- py::class_( m_semi, "UDFTB0Result", "Result of an unrestricted DFTB0 calculation.") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::UDFTB0Result::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::UDFTB0Result::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::UDFTB0Result::e_repulsive) .def_readonly("mo_energies", &vibeqc::semiempirical::UDFTB0Result::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::UDFTB0Result::mo_coeffs) .def_readonly("density_alpha", &vibeqc::semiempirical::UDFTB0Result::density_alpha) .def_readonly("density_beta", &vibeqc::semiempirical::UDFTB0Result::density_beta) .def_readonly("overlap", &vibeqc::semiempirical::UDFTB0Result::overlap) .def_readonly("hamiltonian", &vibeqc::semiempirical::UDFTB0Result::hamiltonian) .def_readonly("n_basis", &vibeqc::semiempirical::UDFTB0Result::n_basis) .def_readonly("n_alpha", &vibeqc::semiempirical::UDFTB0Result::n_alpha) .def_readonly("n_beta", &vibeqc::semiempirical::UDFTB0Result::n_beta); m_semi.def("run_udftb0", &vibeqc::semiempirical::run_udftb0, py::arg("mol"), py::arg("params"), py::call_guard(), "Run unrestricted DFTB0 for open-shell molecules."); py::class_( m_semi, "USCCDFTBResult") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::USCCDFTBResult::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::USCCDFTBResult::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::USCCDFTBResult::e_repulsive) .def_readonly("e_scc", &vibeqc::semiempirical::USCCDFTBResult::e_scc) .def_readonly("mo_energies", &vibeqc::semiempirical::USCCDFTBResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::USCCDFTBResult::mo_coeffs) .def_readonly("density_alpha", &vibeqc::semiempirical::USCCDFTBResult::density_alpha) .def_readonly("density_beta", &vibeqc::semiempirical::USCCDFTBResult::density_beta) .def_readonly("overlap", &vibeqc::semiempirical::USCCDFTBResult::overlap) .def_readonly("hamiltonian", &vibeqc::semiempirical::USCCDFTBResult::hamiltonian) .def_readonly("charges", &vibeqc::semiempirical::USCCDFTBResult::charges) .def_readonly("n_basis", &vibeqc::semiempirical::USCCDFTBResult::n_basis) .def_readonly("n_alpha", &vibeqc::semiempirical::USCCDFTBResult::n_alpha) .def_readonly("n_beta", &vibeqc::semiempirical::USCCDFTBResult::n_beta) .def_readonly("n_iter", &vibeqc::semiempirical::USCCDFTBResult::n_iter) .def_readonly("converged", &vibeqc::semiempirical::USCCDFTBResult::converged); m_semi.def("run_uscc_dftb", &vibeqc::semiempirical::run_uscc_dftb, py::arg("mol"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::SCCOptions{}, py::call_guard(), "Run unrestricted SCC-DFTB for open-shell molecules."); m_semi.def("compute_udftb0_gradient", &vibeqc::semiempirical::compute_udftb0_gradient, py::arg("mol"), py::arg("result"), py::arg("params"), py::call_guard(), "Analytic gradient for unrestricted DFTB0."); // --- k-point DFTB0 (Stage 10) --- py::class_( m_semi, "KPointDFTB0Result") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::KPointDFTB0Result::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::KPointDFTB0Result::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::KPointDFTB0Result::e_repulsive) .def_readonly("band_energies", &vibeqc::semiempirical::KPointDFTB0Result::band_energies) .def_readonly("eps_per_k", &vibeqc::semiempirical::KPointDFTB0Result::eps_per_k) .def_readonly("n_basis", &vibeqc::semiempirical::KPointDFTB0Result::n_basis) .def_readonly("n_occ", &vibeqc::semiempirical::KPointDFTB0Result::n_occ) .def_readonly("n_kpoints", &vibeqc::semiempirical::KPointDFTB0Result::n_kpoints); m_semi.def("run_dftb0_kpoints", &vibeqc::semiempirical::run_dftb0_kpoints, py::arg("system"), py::arg("params"), py::arg("kmesh"), py::arg("cutoff_bohr") = 15.0, "Run periodic DFTB0 on a k-point mesh."); m_semi.def("run_dftb0_bandpath", &vibeqc::semiempirical::run_dftb0_bandpath, py::arg("system"), py::arg("params"), py::arg("kpath"), py::arg("cutoff_bohr") = 15.0, "Run DFTB0 along a band-structure k-path."); py::class_(m_semi, "KPointSCCDFTBResult") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::KPointSCCDFTBResult::energy) .def_readonly("e_scc", &vibeqc::semiempirical::KPointSCCDFTBResult::e_scc) .def_readonly("charges", &vibeqc::semiempirical::KPointSCCDFTBResult::charges) .def_readonly("band_energies", &vibeqc::semiempirical::KPointSCCDFTBResult::band_energies) .def_readonly("eps_per_k", &vibeqc::semiempirical::KPointSCCDFTBResult::eps_per_k) .def_readonly("n_kpoints", &vibeqc::semiempirical::KPointSCCDFTBResult::n_kpoints) .def_readonly("n_iter", &vibeqc::semiempirical::KPointSCCDFTBResult::n_iter) .def_readonly("converged", &vibeqc::semiempirical::KPointSCCDFTBResult::converged); m_semi.def("run_scc_dftb_kpoints", &vibeqc::semiempirical::run_scc_dftb_kpoints, py::arg("system"), py::arg("params"), py::arg("kmesh"), py::arg("scc_opts") = vibeqc::semiempirical::SCCOptions{}, py::arg("cutoff_bohr") = 15.0, "Run periodic SCC-DFTB on a k-point mesh."); // --- Method registry --- py::enum_(m_semi, "MethodFamily") .value("DFTB", vibeqc::semiempirical::MethodFamily::DFTB) .value("XTB", vibeqc::semiempirical::MethodFamily::XTB) .value("NDDO", vibeqc::semiempirical::MethodFamily::NDDO) .value("ML", vibeqc::semiempirical::MethodFamily::ML) .value("LEGACY", vibeqc::semiempirical::MethodFamily::LEGACY); py::enum_(m_semi, "PeriodicTier") .value("NATIVE", vibeqc::semiempirical::PeriodicTier::Native) .value("GENERALIZED", vibeqc::semiempirical::PeriodicTier::Generalized) .value("EXPERIMENTAL", vibeqc::semiempirical::PeriodicTier::Experimental); py::class_( m_semi, "SemiempiricalMethodConfig") .def(py::init<>()) .def_readwrite("name", &vibeqc::semiempirical::SemiempiricalMethodConfig::name) .def_readwrite("display_name", &vibeqc::semiempirical::SemiempiricalMethodConfig::display_name) .def_readwrite("family", &vibeqc::semiempirical::SemiempiricalMethodConfig::family) .def_readwrite("description", &vibeqc::semiempirical::SemiempiricalMethodConfig::description) .def_readwrite("supports_open_shell", &vibeqc::semiempirical::SemiempiricalMethodConfig::supports_open_shell) .def_readwrite("supports_periodic", &vibeqc::semiempirical::SemiempiricalMethodConfig::supports_periodic) .def_readwrite("periodic_tier", &vibeqc::semiempirical::SemiempiricalMethodConfig::periodic_tier) .def_readwrite("gradient_quality", &vibeqc::semiempirical::SemiempiricalMethodConfig::gradient_quality) .def_readwrite("supports_stress", &vibeqc::semiempirical::SemiempiricalMethodConfig::supports_stress) .def_readwrite("supports_kpoints", &vibeqc::semiempirical::SemiempiricalMethodConfig::supports_kpoints) .def_readwrite("has_dispersion", &vibeqc::semiempirical::SemiempiricalMethodConfig::has_dispersion) .def_readwrite("parameter_version", &vibeqc::semiempirical::SemiempiricalMethodConfig::parameter_version); py::class_( m_semi, "CoreElementData") .def(py::init<>()) .def_readwrite("Z", &vibeqc::semiempirical::CoreElementData::Z) .def_readwrite("valence_electrons", &vibeqc::semiempirical::CoreElementData::valence_electrons) .def_readwrite("hubbard_u", &vibeqc::semiempirical::CoreElementData::hubbard_u) .def_readonly("on_site", &vibeqc::semiempirical::CoreElementData::on_site) .def_readonly("zeta", &vibeqc::semiempirical::CoreElementData::zeta); py::class_>( m_semi, "CoreParameterSet", "Base class for semiempirical parameter sets.") .def("has_element", &vibeqc::semiempirical::CoreParameterSet::has_element) .def("element", &vibeqc::semiempirical::CoreParameterSet::element, py::return_value_policy::reference_internal) .def("repulsive_energy", &vibeqc::semiempirical::CoreParameterSet::repulsive_energy) .def("repulsive_derivative", &vibeqc::semiempirical::CoreParameterSet::repulsive_derivative) .def("method_name", &vibeqc::semiempirical::CoreParameterSet::method_name) .def("parameter_version", &vibeqc::semiempirical::CoreParameterSet::parameter_version) .def("n_elements", &vibeqc::semiempirical::CoreParameterSet::n_elements); py::class_( m_semi, "ParameterSetMetadata") .def(py::init<>()) .def_readwrite("method_name", &vibeqc::semiempirical::ParameterSetMetadata::method_name) .def_readwrite("version", &vibeqc::semiempirical::ParameterSetMetadata::version) .def_readwrite("origin", &vibeqc::semiempirical::ParameterSetMetadata::origin) .def_readwrite("license", &vibeqc::semiempirical::ParameterSetMetadata::license) .def_readwrite("n_elements", &vibeqc::semiempirical::ParameterSetMetadata::n_elements) .def_readwrite("doi_or_url", &vibeqc::semiempirical::ParameterSetMetadata::doi_or_url) .def_readwrite("parameter_hash", &vibeqc::semiempirical::ParameterSetMetadata::parameter_hash); py::class_( m_semi, "SemiempiricalMethodRegistry") .def_static("register_dftb", []() { auto& reg = vibeqc::semiempirical::SemiempiricalMethodRegistry::instance(); // DFTB0 { vibeqc::semiempirical::SemiempiricalMethodPlugin p; p.config.name = "dftb0"; p.config.display_name = "DFTB0"; p.config.family = vibeqc::semiempirical::MethodFamily::DFTB; p.config.description = "Non-self-consistent tight-binding"; p.config.supports_open_shell = true; p.config.supports_periodic = true; p.config.periodic_tier = vibeqc::semiempirical::PeriodicTier::Native; p.config.gradient_quality = vibeqc::semiempirical::GradientQuality::Exact; p.config.supports_stress = true; p.config.supports_kpoints = true; p.config.has_dispersion = true; p.config.parameter_version = "in-house-85-2026"; reg.register_method(p); } // PM6 { vibeqc::semiempirical::SemiempiricalMethodPlugin p; p.config.name = "pm6"; p.config.display_name = "PM6"; p.config.family = vibeqc::semiempirical::MethodFamily::NDDO; p.config.description = "Parameterized Model 6 (Stewart 2007)"; p.config.supports_open_shell = true; p.config.gradient_quality = vibeqc::semiempirical::GradientQuality::FiniteDifference; p.config.periodic_tier = vibeqc::semiempirical::PeriodicTier::Experimental; p.config.parameter_version = "pm6-2007"; reg.register_method(p); } // OM1 { vibeqc::semiempirical::SemiempiricalMethodPlugin p; p.config.name = "om1"; p.config.display_name = "OM1"; p.config.family = vibeqc::semiempirical::MethodFamily::NDDO; p.config.description = "Orthogonalization Model 1 (Kolb/Thiel 1993)"; p.config.supports_open_shell = true; p.config.gradient_quality = vibeqc::semiempirical::GradientQuality::FiniteDifference; p.config.parameter_version = "om1-1993"; reg.register_method(p); } // OM2 { vibeqc::semiempirical::SemiempiricalMethodPlugin p; p.config.name = "om2"; p.config.display_name = "OM2"; p.config.family = vibeqc::semiempirical::MethodFamily::NDDO; p.config.description = "Orthogonalization Model 2 (Weber/Thiel 2000)"; p.config.supports_open_shell = true; p.config.gradient_quality = vibeqc::semiempirical::GradientQuality::FiniteDifference; p.config.parameter_version = "om2-2000"; reg.register_method(p); } // OM3 { vibeqc::semiempirical::SemiempiricalMethodPlugin p; p.config.name = "om3"; p.config.display_name = "OM3"; p.config.family = vibeqc::semiempirical::MethodFamily::NDDO; p.config.description = "Orthogonalization Model 3 (Scholten 2003)"; p.config.supports_open_shell = true; p.config.gradient_quality = vibeqc::semiempirical::GradientQuality::FiniteDifference; p.config.parameter_version = "om3-2003"; reg.register_method(p); } // GFN2-xTB { vibeqc::semiempirical::SemiempiricalMethodPlugin p; p.config.name = "gfn2-xtb"; p.config.display_name = "GFN2-xTB"; p.config.family = vibeqc::semiempirical::MethodFamily::XTB; p.config.description = "Extended tight-binding (Grimme 2019)"; p.config.supports_open_shell = true; p.config.supports_periodic = true; p.config.periodic_tier = vibeqc::semiempirical::PeriodicTier::Generalized; p.config.gradient_quality = vibeqc::semiempirical::GradientQuality::Approximate; p.config.supports_stress = false; p.config.supports_kpoints = false; p.config.has_dispersion = true; p.config.parameter_version = "gfn2-xtb-2019"; reg.register_method(p); } // SCC-DFTB { vibeqc::semiempirical::SemiempiricalMethodPlugin p; p.config.name = "scc-dftb"; p.config.display_name = "SCC-DFTB"; p.config.family = vibeqc::semiempirical::MethodFamily::DFTB; p.config.description = "Self-consistent-charge tight-binding"; p.config.supports_open_shell = true; p.config.supports_periodic = true; p.config.periodic_tier = vibeqc::semiempirical::PeriodicTier::Native; p.config.gradient_quality = vibeqc::semiempirical::GradientQuality::Approximate; p.config.supports_stress = true; p.config.supports_kpoints = true; p.config.has_dispersion = true; p.config.parameter_version = "in-house-85-2026"; reg.register_method(p); } }) .def_static("instance", &vibeqc::semiempirical::SemiempiricalMethodRegistry::instance, py::return_value_policy::reference) .def("find_config", [](const vibeqc::semiempirical::SemiempiricalMethodRegistry& reg, const std::string& name) -> vibeqc::semiempirical::SemiempiricalMethodConfig { const auto* p = reg.find(name); if (p) return p->config; throw std::runtime_error("Method not found: " + name); }) .def("method_names", &vibeqc::semiempirical::SemiempiricalMethodRegistry::method_names) .def("family_methods", &vibeqc::semiempirical::SemiempiricalMethodRegistry::family_methods); // --- GFN2-xTB parameters --- py::module_ m_xtb = m_semi.def_submodule("xtb", "GFN-xTB method family"); py::class_( m_xtb, "GFN2ShellData", "Per-shell data for GFN2-xTB.") .def(py::init<>()) .def_readwrite("n", &vibeqc::semiempirical::xtb::GFN2ElementData::ShellData::n) .def_readwrite("l", &vibeqc::semiempirical::xtb::GFN2ElementData::ShellData::l) .def_readwrite("en", &vibeqc::semiempirical::xtb::GFN2ElementData::ShellData::en) .def_readwrite("zeta", &vibeqc::semiempirical::xtb::GFN2ElementData::ShellData::zeta) .def_readwrite("k_en", &vibeqc::semiempirical::xtb::GFN2ElementData::ShellData::k_en); py::class_( m_xtb, "GFN2ElementData", "Per-element data for GFN2-xTB.") .def(py::init<>()) .def_readwrite("Z", &vibeqc::semiempirical::xtb::GFN2ElementData::Z) .def_property("shells", [](vibeqc::semiempirical::xtb::GFN2ElementData& ed) -> std::vector& { return ed.shells; }, [](vibeqc::semiempirical::xtb::GFN2ElementData& ed, const std::vector& v) { ed.shells = v; }) .def_readwrite("gam", &vibeqc::semiempirical::xtb::GFN2ElementData::gam) .def_readwrite("alpha", &vibeqc::semiempirical::xtb::GFN2ElementData::alpha) .def("add_shell", [](vibeqc::semiempirical::xtb::GFN2ElementData& ed, int l, double en, double zeta, double k_en) { vibeqc::semiempirical::xtb::GFN2ElementData::ShellData sd; sd.l = l; sd.en = en; sd.zeta = zeta; sd.k_en = k_en; ed.shells.push_back(sd); }); py::class_( m_xtb, "GFN2RepulsivePair", "Repulsive pair data for GFN2-xTB.") .def(py::init<>()) .def_readwrite("alpha", &vibeqc::semiempirical::xtb::GFN2RepulsivePair::alpha) .def_readwrite("k_ab", &vibeqc::semiempirical::xtb::GFN2RepulsivePair::k_ab); py::class_( m_xtb, "GFN2ParameterSet", "GFN2-xTB parameter set (Grimme group, 2019).") .def(py::init<>()) .def("has_element", &vibeqc::semiempirical::xtb::GFN2ParameterSet::has_element) .def("n_elements", &vibeqc::semiempirical::xtb::GFN2ParameterSet::n_elements) .def("repulsive_energy", &vibeqc::semiempirical::xtb::GFN2ParameterSet::repulsive_energy) .def("add_element", &vibeqc::semiempirical::xtb::GFN2ParameterSet::add_element) .def("set_repulsive_pair", &vibeqc::semiempirical::xtb::GFN2ParameterSet::set_repulsive_pair) .def("metadata", &vibeqc::semiempirical::xtb::GFN2ParameterSet::metadata); // --- GFN2-xTB driver --- py::class_( m_xtb, "GFN2Result", "Result of a GFN2-xTB calculation.") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::xtb::GFN2Result::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::xtb::GFN2Result::e_electronic) .def_readonly("e_repulsive", &vibeqc::semiempirical::xtb::GFN2Result::e_repulsive) .def_readonly("e_scc", &vibeqc::semiempirical::xtb::GFN2Result::e_scc) .def_readonly("charges", &vibeqc::semiempirical::xtb::GFN2Result::charges) .def_readonly("converged", &vibeqc::semiempirical::xtb::GFN2Result::converged) .def_readonly("n_iter", &vibeqc::semiempirical::xtb::GFN2Result::n_iter); py::class_( m_xtb, "XTBSccOptions") .def(py::init<>()) .def_readwrite("max_iter", &vibeqc::semiempirical::xtb::XTBSccOptions::max_iter) .def_readwrite("conv_tol_charge", &vibeqc::semiempirical::xtb::XTBSccOptions::conv_tol_charge) .def_readwrite("charge_mixing", &vibeqc::semiempirical::xtb::XTBSccOptions::charge_mixing); m_xtb.def("run_gfn2_xtb", &vibeqc::semiempirical::xtb::run_gfn2_xtb, py::arg("mol"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::xtb::XTBSccOptions{}, py::call_guard(), "Run GFN2-xTB SCC energy calculation."); // GFN2 gradient m_semi.def("compute_gfn2_gradient", &vibeqc::semiempirical::compute_gfn2_gradient, py::arg("mol"), py::arg("result"), py::arg("params"), py::call_guard(), "GFN2-xTB analytic gradient (fixed-charge approximation)."); // Periodic GFN2-xTB py::class_( m_xtb, "PeriodicGFN2Result") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::xtb::PeriodicGFN2Result::energy) .def_readonly("e_repulsive", &vibeqc::semiempirical::xtb::PeriodicGFN2Result::e_repulsive) .def_readonly("converged", &vibeqc::semiempirical::xtb::PeriodicGFN2Result::converged) .def_readonly("n_iter", &vibeqc::semiempirical::xtb::PeriodicGFN2Result::n_iter) .def_readonly("n_cells", &vibeqc::semiempirical::xtb::PeriodicGFN2Result::n_cells); m_xtb.def("run_gfn2_xtb_gamma", &vibeqc::semiempirical::xtb::run_gfn2_xtb_gamma, py::arg("system"), py::arg("params"), py::arg("scc_opts") = vibeqc::semiempirical::xtb::XTBSccOptions{}, py::arg("cutoff_bohr") = 15.0, py::call_guard(), "Run periodic GFN2-xTB SCC calculation at Gamma-point."); m_semi.def("compute_periodic_gfn2_gradient", &vibeqc::semiempirical::compute_periodic_gfn2_gradient, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic GFN2-xTB atomic gradient at Gamma-point (fixed-charge approx)."); m_semi.def("compute_periodic_gfn2_stress", &vibeqc::semiempirical::compute_periodic_gfn2_stress, py::arg("system"), py::arg("result"), py::arg("params"), py::call_guard(), "Periodic GFN2-xTB stress tensor at Gamma-point (fixed-charge approx)."); // Unrestricted GFN2-xTB py::class_( m_xtb, "UGFN2Result") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::xtb::UGFN2Result::energy) .def_readonly("n_alpha", &vibeqc::semiempirical::xtb::UGFN2Result::n_alpha) .def_readonly("n_beta", &vibeqc::semiempirical::xtb::UGFN2Result::n_beta) .def_readonly("converged", &vibeqc::semiempirical::xtb::UGFN2Result::converged); m_xtb.def("run_ugfn2_xtb", &vibeqc::semiempirical::xtb::run_ugfn2_xtb, py::arg("mol"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::xtb::XTBSccOptions{}, py::call_guard(), "Run unrestricted GFN2-xTB SCC calculation."); // --- NDDO methods --- py::module_ m_nddo = m_semi.def_submodule("nddo", "NDDO-family methods (MNDO/AM1/PMx/OMx)"); py::class_( m_nddo, "NDDOElementData") .def(py::init<>()) .def_readwrite("Z", &vibeqc::semiempirical::nddo::NDDOElementData::Z) .def_readwrite("uss", &vibeqc::semiempirical::nddo::NDDOElementData::uss) .def_readwrite("upp", &vibeqc::semiempirical::nddo::NDDOElementData::upp) .def_readwrite("betas", &vibeqc::semiempirical::nddo::NDDOElementData::betas) .def_readwrite("betap", &vibeqc::semiempirical::nddo::NDDOElementData::betap) .def_readwrite("zs", &vibeqc::semiempirical::nddo::NDDOElementData::zs) .def_readwrite("zp", &vibeqc::semiempirical::nddo::NDDOElementData::zp) .def_readwrite("gss", &vibeqc::semiempirical::nddo::NDDOElementData::gss) .def_readwrite("gpp", &vibeqc::semiempirical::nddo::NDDOElementData::gpp) .def_readwrite("gsp", &vibeqc::semiempirical::nddo::NDDOElementData::gsp) .def_readwrite("gp2", &vibeqc::semiempirical::nddo::NDDOElementData::gp2) .def_readwrite("hsp", &vibeqc::semiempirical::nddo::NDDOElementData::hsp) .def_readwrite("alpha", &vibeqc::semiempirical::nddo::NDDOElementData::alpha) .def_readwrite("polvo", &vibeqc::semiempirical::nddo::NDDOElementData::polvo) .def("add_gamma_term", [](vibeqc::semiempirical::nddo::NDDOElementData& ed, double coeff, double exponent, double factor) { ed.gamma_terms.push_back({coeff, exponent, factor}); }) .def("add_gaussian_term", &vibeqc::semiempirical::nddo::NDDOElementData::add_gaussian_term, py::arg("coeff"), py::arg("exponent"), py::arg("factor")); py::class_( m_nddo, "PM6ParameterSet") .def(py::init<>()) .def("has_element", &vibeqc::semiempirical::nddo::PM6ParameterSet::has_element) .def("n_elements", &vibeqc::semiempirical::nddo::PM6ParameterSet::n_elements) .def("add_element", &vibeqc::semiempirical::nddo::PM6ParameterSet::add_element) .def("add_diatomic", [](vibeqc::semiempirical::nddo::PM6ParameterSet& p, int Z1, int Z2, double alpb, double xfac) { vibeqc::semiempirical::nddo::NDDODiatomicParams dp; dp.d1 = alpb; dp.d2 = xfac; p.set_diatomic(Z1, Z2, dp); }, py::arg("Z1"), py::arg("Z2"), py::arg("alpb"), py::arg("xfac")) .def("metadata", &vibeqc::semiempirical::nddo::PM6ParameterSet::metadata); // --- OMx parameter set --- py::class_( m_nddo, "OMxElementData") .def(py::init<>()) .def_readwrite("Z", &vibeqc::semiempirical::nddo::OMxElementData::Z) .def_readwrite("uss", &vibeqc::semiempirical::nddo::OMxElementData::uss) .def_readwrite("upp", &vibeqc::semiempirical::nddo::OMxElementData::upp) .def_readwrite("gss", &vibeqc::semiempirical::nddo::OMxElementData::gss) .def_readwrite("gpp", &vibeqc::semiempirical::nddo::OMxElementData::gpp) .def_readwrite("gsp", &vibeqc::semiempirical::nddo::OMxElementData::gsp) .def_readwrite("gp2", &vibeqc::semiempirical::nddo::OMxElementData::gp2) .def_readwrite("hsp", &vibeqc::semiempirical::nddo::OMxElementData::hsp) .def_readwrite("beta_s", &vibeqc::semiempirical::nddo::OMxElementData::beta_s) .def_readwrite("beta_p", &vibeqc::semiempirical::nddo::OMxElementData::beta_p) .def_readwrite("beta_pi", &vibeqc::semiempirical::nddo::OMxElementData::beta_pi) .def_readwrite("alpha_s", &vibeqc::semiempirical::nddo::OMxElementData::alpha_s) .def_readwrite("alpha_p", &vibeqc::semiempirical::nddo::OMxElementData::alpha_p) .def_readwrite("alpha_pi", &vibeqc::semiempirical::nddo::OMxElementData::alpha_pi) .def_readwrite("F1", &vibeqc::semiempirical::nddo::OMxElementData::F1) .def_readwrite("G1", &vibeqc::semiempirical::nddo::OMxElementData::G1) .def_readwrite("F2", &vibeqc::semiempirical::nddo::OMxElementData::F2) .def_readwrite("G2", &vibeqc::semiempirical::nddo::OMxElementData::G2) .def_readwrite("zeta", &vibeqc::semiempirical::nddo::OMxElementData::zeta) .def_readwrite("zs", &vibeqc::semiempirical::nddo::OMxElementData::zs) .def_readwrite("zp", &vibeqc::semiempirical::nddo::OMxElementData::zp) .def_readwrite("alpha", &vibeqc::semiempirical::nddo::OMxElementData::alpha) .def_readwrite("n_orbitals", &vibeqc::semiempirical::nddo::OMxElementData::n_orbitals) .def_readwrite("beta_s_xh", &vibeqc::semiempirical::nddo::OMxElementData::beta_s_xh) .def_readwrite("beta_p_xh", &vibeqc::semiempirical::nddo::OMxElementData::beta_p_xh) .def_readwrite("alpha_s_xh", &vibeqc::semiempirical::nddo::OMxElementData::alpha_s_xh) .def_readwrite("alpha_p_xh", &vibeqc::semiempirical::nddo::OMxElementData::alpha_p_xh) .def_readwrite("zeta_alpha", &vibeqc::semiempirical::nddo::OMxElementData::zeta_alpha) .def_readwrite("F_alpha_alpha", &vibeqc::semiempirical::nddo::OMxElementData::F_alpha_alpha) .def_readwrite("beta_alpha", &vibeqc::semiempirical::nddo::OMxElementData::beta_alpha) .def_readwrite("alpha_alpha", &vibeqc::semiempirical::nddo::OMxElementData::alpha_alpha); py::class_( m_nddo, "OMxParameterSet") .def(py::init(), py::arg("variant")) .def("variant", &vibeqc::semiempirical::nddo::OMxParameterSet::variant) .def("method_name", &vibeqc::semiempirical::nddo::OMxParameterSet::method_name) .def("add_omx_element", &vibeqc::semiempirical::nddo::OMxParameterSet::add_omx_element); // --- PM6 driver --- py::class_( m_nddo, "PM6Result") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::nddo::PM6Result::energy) .def_readonly("e_core", &vibeqc::semiempirical::nddo::PM6Result::e_core) .def_readonly("converged", &vibeqc::semiempirical::nddo::PM6Result::converged) .def_readonly("n_iter", &vibeqc::semiempirical::nddo::PM6Result::n_iter); m_nddo.def("compute_pm6_gradient_fd", &vibeqc::semiempirical::nddo::compute_pm6_gradient_fd, py::arg("mol"), py::arg("params"), py::arg("h") = 0.001, py::call_guard(), "Finite-difference PM6 gradient (Hartree/bohr)."); m_nddo.def("run_pm6", &vibeqc::semiempirical::nddo::run_pm6, py::arg("mol"), py::arg("params"), py::arg("max_iter") = 100, py::arg("conv_tol") = 1e-7, py::call_guard(), "Run PM6 SCF energy calculation."); // --- OMx driver --- m_nddo.def("run_omx", &vibeqc::semiempirical::nddo::run_omx, py::arg("mol"), py::arg("params"), py::arg("max_iter") = 100, py::arg("conv_tol") = 1e-7, py::call_guard(), "Run OMx SCF energy calculation (Loewdin orthogonalization)."); m_nddo.def("run_omx_v2", &vibeqc::semiempirical::nddo::run_omx_v2, py::arg("mol"), py::arg("params"), py::arg("max_iter") = 100, py::arg("conv_tol") = 1e-7, py::call_guard(), "Run OMx SCF with the proper OM1/OM2/OM3 Hamiltonian " "(exponential resonance integrals + VORT corrections)."); m_nddo.def("compute_omx_gradient_fd", &vibeqc::semiempirical::nddo::compute_omx_gradient_fd, py::arg("mol"), py::arg("params"), py::arg("h") = 0.001, py::call_guard(), "Finite-difference OMx gradient (v1, Hartree/bohr)."); m_nddo.def("compute_omx_v2_gradient_fd", &vibeqc::semiempirical::nddo::compute_omx_v2_gradient_fd, py::arg("mol"), py::arg("params"), py::arg("h") = 0.001, py::call_guard(), "Finite-difference OMx v2 gradient (Hartree/bohr)."); // OMx py::enum_(m_nddo, "OMxVariant") .value("OM1", vibeqc::semiempirical::nddo::OMxVariant::OM1) .value("OM2", vibeqc::semiempirical::nddo::OMxVariant::OM2) .value("OM3", vibeqc::semiempirical::nddo::OMxVariant::OM3); // --- Periodic PM6 --- py::class_( m_nddo, "PeriodicPM6Options") .def(py::init<>()) .def_readwrite("cutoff_bohr", &vibeqc::semiempirical::nddo::PeriodicPM6Options::cutoff_bohr) .def_readwrite("conv_tol", &vibeqc::semiempirical::nddo::PeriodicPM6Options::conv_tol) .def_readwrite("max_iter", &vibeqc::semiempirical::nddo::PeriodicPM6Options::max_iter) .def_readwrite("gamma_only_0", &vibeqc::semiempirical::nddo::PeriodicPM6Options::gamma_only_0); py::class_( m_nddo, "PeriodicPM6Result") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::nddo::PeriodicPM6Result::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::nddo::PeriodicPM6Result::e_electronic) .def_readonly("e_core", &vibeqc::semiempirical::nddo::PeriodicPM6Result::e_core) .def_readonly("mo_energies", &vibeqc::semiempirical::nddo::PeriodicPM6Result::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::nddo::PeriodicPM6Result::mo_coeffs) .def_readonly("density", &vibeqc::semiempirical::nddo::PeriodicPM6Result::density) .def_readonly("fock_gamma", &vibeqc::semiempirical::nddo::PeriodicPM6Result::fock_gamma) .def_readonly("n_basis", &vibeqc::semiempirical::nddo::PeriodicPM6Result::n_basis) .def_readonly("n_occ", &vibeqc::semiempirical::nddo::PeriodicPM6Result::n_occ) .def_readonly("n_cells", &vibeqc::semiempirical::nddo::PeriodicPM6Result::n_cells) .def_readonly("n_iter", &vibeqc::semiempirical::nddo::PeriodicPM6Result::n_iter) .def_readonly("converged", &vibeqc::semiempirical::nddo::PeriodicPM6Result::converged); m_nddo.def("run_pm6_gamma", &vibeqc::semiempirical::nddo::run_pm6_gamma, py::arg("system"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::nddo::PeriodicPM6Options{}, py::call_guard(), "Run periodic PM6 SCF calculation at Gamma-point."); // --- Periodic OMx --- py::class_( m_nddo, "PeriodicOMxOptions") .def(py::init<>()) .def_readwrite("cutoff_bohr", &vibeqc::semiempirical::nddo::PeriodicOMxOptions::cutoff_bohr) .def_readwrite("conv_tol", &vibeqc::semiempirical::nddo::PeriodicOMxOptions::conv_tol) .def_readwrite("max_iter", &vibeqc::semiempirical::nddo::PeriodicOMxOptions::max_iter) .def_readwrite("gamma_only_0", &vibeqc::semiempirical::nddo::PeriodicOMxOptions::gamma_only_0) .def_readwrite("warmup_iters", &vibeqc::semiempirical::nddo::PeriodicOMxOptions::warmup_iters) .def_readwrite("density_mixing", &vibeqc::semiempirical::nddo::PeriodicOMxOptions::density_mixing) .def_readwrite("eval_floor", &vibeqc::semiempirical::nddo::PeriodicOMxOptions::eval_floor); py::class_( m_nddo, "PeriodicOMxResult") .def(py::init<>()) .def_readonly("energy", &vibeqc::semiempirical::nddo::PeriodicOMxResult::energy) .def_readonly("e_electronic", &vibeqc::semiempirical::nddo::PeriodicOMxResult::e_electronic) .def_readonly("e_core", &vibeqc::semiempirical::nddo::PeriodicOMxResult::e_core) .def_readonly("mo_energies", &vibeqc::semiempirical::nddo::PeriodicOMxResult::mo_energies) .def_readonly("mo_coeffs", &vibeqc::semiempirical::nddo::PeriodicOMxResult::mo_coeffs) .def_readonly("density", &vibeqc::semiempirical::nddo::PeriodicOMxResult::density) .def_readonly("fock_gamma", &vibeqc::semiempirical::nddo::PeriodicOMxResult::fock_gamma) .def_readonly("overlap_gamma", &vibeqc::semiempirical::nddo::PeriodicOMxResult::overlap_gamma) .def_readonly("n_basis", &vibeqc::semiempirical::nddo::PeriodicOMxResult::n_basis) .def_readonly("n_occ", &vibeqc::semiempirical::nddo::PeriodicOMxResult::n_occ) .def_readonly("n_cells", &vibeqc::semiempirical::nddo::PeriodicOMxResult::n_cells) .def_readonly("n_iter", &vibeqc::semiempirical::nddo::PeriodicOMxResult::n_iter) .def_readonly("converged", &vibeqc::semiempirical::nddo::PeriodicOMxResult::converged); m_nddo.def("run_omx_gamma", &vibeqc::semiempirical::nddo::run_omx_gamma, py::arg("system"), py::arg("params"), py::arg("opts") = vibeqc::semiempirical::nddo::PeriodicOMxOptions{}, py::call_guard(), "Run periodic OMx SCF calculation at Gamma-point."); // --- Hamiltonian builders --- m_semi.def("build_gfn2_hamiltonian_zero", &vibeqc::semiempirical::build_gfn2_hamiltonian_zero, py::arg("basis"), py::arg("S"), py::arg("mol"), py::arg("params"), "Build GFN2-xTB H⁰ from overlap and parameters."); m_semi.def("build_gfn2_gamma", &vibeqc::semiempirical::build_gfn2_gamma, py::arg("mol"), py::arg("params"), "Build GFN2-xTB gamma (Klopman-Ohno) matrix."); // --- Repulsive spline --- py::class_(m_semi, "RepulsiveSpline") .def(py::init<>()) .def("build", &vibeqc::semiempirical::RepulsiveSpline::build) .def("evaluate", &vibeqc::semiempirical::RepulsiveSpline::evaluate) .def("derivative", &vibeqc::semiempirical::RepulsiveSpline::derivative) .def("n_knots", &vibeqc::semiempirical::RepulsiveSpline::n_knots) .def("cutoff", &vibeqc::semiempirical::RepulsiveSpline::cutoff) .def("empty", &vibeqc::semiempirical::RepulsiveSpline::empty); }