diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py index bbd763e0..b19f6618 100644 --- a/src/nomad_simulations/schema_packages/basis_set.py +++ b/src/nomad_simulations/schema_packages/basis_set.py @@ -13,7 +13,7 @@ from nomad import utils from nomad.datamodel.data import ArchiveSection from nomad.datamodel.metainfo.annotations import ELNAnnotation -from nomad.metainfo import MEnum, Quantity, SubSection +from nomad.metainfo import JSON, MEnum, Quantity, SubSection from nomad.units import ureg from nomad_simulations.schema_packages.atoms_state import AtomsState @@ -190,9 +190,119 @@ class AtomCenteredFunction(ArchiveSection): Specifies a single function (term) in an atom-centered basis set. """ - pass + harmonic_type = Quantity( + type=MEnum( + 'spherical', + 'cartesian', + ), + default='spherical', + description=""" + Specifies whether the basis functions are spherical-harmonic or cartesian functions. + """, + ) + + function_type = Quantity( + type=MEnum( + 's', + 'p', + 'd', + 'f', + 'g', + 'h', + 'i', + 'j', + 'k', + 'l', + 'sp', + 'spd', + 'spdf', + ), + description=""" + L=a+b+c + The angular momentum of GTO to be added. + """, + ) + + n_primitive = Quantity( + type=np.int32, + description=""" + Number of primitives. + Linear combinations of the primitive Gaussians are formed to approximate the radial extent of an STO. + """, + ) + + exponents = Quantity( + type=np.float32, + shape=['n_primitive'], + description=""" + List of exponents for the basis function. + """, + ) + + contraction_coefficients = Quantity( + type=np.float32, + shape=['*'], # Flexible shape to handle combined types (e.g. SP, SPD..) + description=""" + List of contraction coefficients corresponding to the exponents. + """, + ) + + point_charge = Quantity( + type=np.float32, + description=""" + the value of the point charge. + """, + ) - # TODO: design system for writing basis functions like gaussian or slater orbitals + def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: + """ + Validates the input data + and resolves combined types like SP, SPD, SPDF, etc. + + Raises ValueError: If the data is inconsistent (e.g., mismatch in exponents and coefficients). + """ + super().normalize(archive, logger) + + # Validate number of primitives + if self.n_primitive is not None: + if self.exponents is not None and len(self.exponents) != self.n_primitive: + raise ValueError( + f'Mismatch in number of exponents: expected {self.n_primitive}, ' + f'found {len(self.exponents)}.' + ) + + # Resolve combined types + if self.function_type and len(self.function_type) > 1: + num_types = len(self.function_type) # For SP: 2, SPD: 3, etc. + if self.contraction_coefficients is not None: + expected_coeffs = num_types * self.n_primitive + if len(self.contraction_coefficients) != expected_coeffs: + raise ValueError( + f'Mismatch in contraction coefficients for {self.function_type} type: ' + f'expected {expected_coeffs}, found {len(self.contraction_coefficients)}.' + ) + + # Split coefficients into separate lists for each type + self.coefficient_sets = { + t: self.contraction_coefficients[i::num_types] + for i, t in enumerate(self.function_type) + } + + # Debug: Log split coefficients + for t, coeffs in self.coefficient_sets.items(): + logger.info(f'{t}-type coefficients: {coeffs}') + else: + logger.warning( + f'No contraction coefficients provided for {self.function_type} type.' + ) + + # For single types, ensure coefficients match primitives + elif self.contraction_coefficients is not None: + if len(self.contraction_coefficients) != self.n_primitive: + raise ValueError( + f'Mismatch in contraction coefficients: expected {self.n_primitive}, ' + f'found {len(self.contraction_coefficients)}.' + ) class AtomCenteredBasisSet(BasisSetComponent): @@ -200,15 +310,50 @@ class AtomCenteredBasisSet(BasisSetComponent): Defines an atom-centered basis set. """ + basis_set = Quantity( + type=str, + description=""" + name of the basis set. + """, + ) + + type = Quantity( + type=MEnum( + 'STO', # Slater-type orbitals + 'GTO', # Gaussian-type orbitals + 'NAO', # Numerical atomic orbitals + 'cECP', # Capped effective core potentials + 'PC', # Point charges + ), + description=""" + Type of the basis set, e.g. STO or GTO. + """, + ) + + role = Quantity( + type=MEnum( + 'orbital', + 'auxiliary_scf', + 'auxiliary_post_hf', + 'cabs', # complementary auxiliary basis set + ), + description=""" + The role of the basis set. + """, + ) + + total_number_of_basis_functions = Quantity( + type=np.int32, + description='', + ) + functional_composition = SubSection( sub_section=AtomCenteredFunction.m_def, repeats=True - ) # TODO change name + ) def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: super().normalize(archive, logger) # self.name = self.m_def.name - # TODO: set name based on basis functions - # ? use basis set names from Basis Set Exchange class APWBaseOrbital(ArchiveSection): diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py index c7b143ff..c6700e1f 100644 --- a/src/nomad_simulations/schema_packages/model_method.py +++ b/src/nomad_simulations/schema_packages/model_method.py @@ -1221,3 +1221,259 @@ class DMFT(ModelMethodElectronic): def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: super().normalize(archive, logger) + + +class MolecularOrbital(ArchiveSection): + """ + A section representing a single molecular orbital. + """ + + energy = Quantity( + type=np.float64, + # unit='electron_volt', + description='Energy of the molecular orbital.', + ) + + occupation = Quantity( + type=np.float64, description='Occupation number of the molecular orbital.' + ) + + symmetry_label = Quantity( + type=str, description='Symmetry label of the molecular orbital.' + ) + + reference_orbital = SubSection( + sub_section=OrbitalsState.m_def, + description='Reference to the underlying atomic orbital state.', + ) + + +class LCAO(ArchiveSection): + """ + A base class for molecular orbital schemes used in quantum chemistry calculations. + Supports unrestricted (UHF, UKS), restricted (RHF, RKS), and restricted open-shell (ROHF, ROKS) schemes. + """ + + reference_type = Quantity( + type=MEnum('RHF', 'ROHF', 'UHF', 'RKS', 'ROKS', 'UKS'), + description=""" + Specifies the type of reference wavefunction: + - RHF: Restricted Hartree-Fock + - ROHF: Restricted Open-Shell Hartree-Fock + - UHF: Unrestricted Hartree-Fock + - RKS: Restricted Kohn-Sham + - ROKS: Restricted Open-Shell Kohn-Sham + - UKS: Unrestricted Kohn-Sham + """, + a_eln=ELNAnnotation(component='EnumEditQuantity'), + ) + + orbital_set = Quantity( + type=MEnum('canonical', 'natural', 'localized'), + description=""" + Specifies the type of orbitals used in the molecular orbital scheme: + - canonical: Default canonical molecular orbitals. + - natural: Natural orbitals obtained from the density matrix. + - localized: Localized orbitals such as Boys or Foster-Boys localization. + TODO: this will be later connected to many other things. + """, + a_eln=ELNAnnotation(component='EnumEditQuantity'), + ) + + n_alpha_electrons = Quantity( + type=np.int32, + description=""" + Number of alpha electrons (spin up) in the molecular system. + """, + ) + + n_beta_electrons = Quantity( + type=np.int32, + description=""" + Number of beta electrons (spin down) in the molecular system. + """, + ) + + n_molecular_orbitals = Quantity( + type=np.int32, + description=""" + Total number of molecular orbitals in the system. + """, + ) + + molecular_orbitals = SubSection( + sub_section=MolecularOrbital.m_def, + repeats=True, + description=""" + Detailed information about each molecular orbital, + including energy, occupation, and symmetry label. + """, + ) + + total_spin = Quantity( + type=np.float64, + description=""" + Total spin of the system defined as S = (n_alpha_electrons - n_beta_electrons) / 2. + Connect to the model system. + """, + ) + + spin_multiplicity = Quantity( + type=np.int32, + description=""" + Spin multiplicity of the system defined as 2S + 1. + """, + ) + + def resolve_spin_properties(self, logger: 'BoundLogger') -> None: + """ + Resolves the total spin and spin multiplicity of the system based on alpha and beta electrons. + """ + if self.n_alpha_electrons is not None and self.n_beta_electrons is not None: + self.total_spin = (self.n_alpha_electrons - self.n_beta_electrons) / 2 + self.spin_multiplicity = int(2 * self.total_spin + 1) + + def validate_scheme(self, logger: 'BoundLogger') -> bool: + """ + Validates the consistency of the molecular orbital scheme. + + Returns: + (bool): True if the scheme is consistent, False otherwise. + """ + if self.reference_type in ['RHF', 'RKS']: + if self.n_alpha_electrons != self.n_beta_electrons: + logger.error( + f'For {self.reference_type}, the number of alpha and beta electrons must be equal.' + ) + return False + if self.reference_type in ['ROHF', 'ROKS']: + if abs(self.n_alpha_electrons - self.n_beta_electrons) != 1: + logger.error( + f'For {self.reference_type}, there must be exactly one unpaired electron.' + ) + return False + return True + + def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: + super().normalize(archive, logger) + + # Resolve spin properties + self.resolve_spin_properties(logger) + + # Validate the molecular orbital scheme + if not self.validate_scheme(logger): + logger.error('Invalid molecular orbital scheme.') + + # Resolve the number of molecular orbitals + if self.n_molecular_orbitals is None and self.molecular_orbitals: + self.n_molecular_orbitals = len(self.molecular_orbitals) + + # Validate molecular orbital occupation + total_occupation = sum( + orbital.occupation + for orbital in self.molecular_orbitals + if orbital.occupation is not None + ) + expected_occupation = self.n_alpha_electrons + self.n_beta_electrons + if total_occupation != expected_occupation: + logger.warning( + f'The total occupation ({total_occupation}) does not match the expected value ({expected_occupation}).' + ) + + +class GTOIntegralDecomposition(BaseModelMethod): + """ + A general class for integral decomposition techniques for Coulomb and exchange integrals. + Examples: + Resolution of identity (RI-approximation): + RI + RIJK + .... + Chain-of-spheres (COSX) algorithm for exchange: doi:10.1016/j.chemphys.2008.10.036 + """ + + approximation_type = Quantity( + type=str, + description=""" + RIJ, RIK, RIJK, + """, + ) + + approximated_term = Quantity( + type=str, + description=""" + such as coulomb, exchange, explicit-correlation + to be converted to an MEnum later. + """, + ) + + def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: + super().normalize(archive, logger) + + +class HartreeFock(ModelMethodElectronic): + """ + A base section for Hartree Fock ansatz. + """ + + reference_determinant = Quantity( + type=MEnum('UHF', 'RHF', 'ROHF', 'UKS', 'RKS', 'ROKS'), + description=""" + the type of reference determinant. + """, + ) + + +class PerturbationMethod(ModelMethodElectronic): + type = Quantity( + type=MEnum('MP', 'RS', 'BW'), + description=""" + Perturbation approach. The abbreviations stand for: + | Abbreviation | Description | + | ------------ | ----------- | + | `'MP'` | Moller-Plesset | + | `'RS'` | Rayleigh-Schrödigner | + | `'BW'` | Brillouin-Wigner | + """, + a_eln=ELNAnnotation(component='EnumEditQuantity'), + ) # TODO: check if the special symbols are supported + + order = Quantity( + type=np.int32, + description=""" + Order up to which the perturbation is expanded. + """, + a_eln=ELNAnnotation(component='NumberEditQuantity'), + ) + + density = Quantity( + type=MEnum('relaxed', 'unrelaxed'), + description=""" + unrelaxed density: MP2 expectation value density + relaxed density : incorporates orbital relaxation + """, + ) + + +class LocalCorrelation(ArchiveSection): + """ + A base section used to define the parameters of a local correlation for the post-HF methods. + + It has a corresponding localization method. + """ + + method = Quantity( + type=MEnum('LMP2', 'LCCD', 'LCCSD', 'LCCSD(T)', 'DLPNO-CCSD(T)', 'LocalDFT'), + description=""" + The local correlation method applied. For example: + - LMP2: Local Møller-Plesset perturbation theory + - LCCD: Local Coupled-Cluster with Doubles + - LCCSD: Local Coupled-Cluster with Singles and Doubles + - LCCSD(T): Local Coupled-Cluster with Singles, Doubles, and Perturbative Triples + - DLPNO-CCSD(T): Domain-Based Local Pair Natural Orbital CCSD(T) + - LocalDFT: Local Density Functional Theory. + + # TODO: improve list! + """, + a_eln=ELNAnnotation(component='EnumEditQuantity'), + ) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 66364d5e..33ee56e7 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -1113,6 +1113,21 @@ class ModelSystem(System): """, ) + total_charge = Quantity( + type=np.int32, + description=""" + Total charge of the system. + """, + ) + + total_spin = Quantity( + type=np.int32, + description=""" + Total spin state of the system. + Not to be confused with the spin multiplicity 2S + 1. + """, + ) + model_system = SubSection(sub_section=SectionProxy('ModelSystem'), repeats=True) def resolve_system_type_and_dimensionality( diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py index 515deea7..48ec6dde 100644 --- a/src/nomad_simulations/schema_packages/numerical_settings.py +++ b/src/nomad_simulations/schema_packages/numerical_settings.py @@ -52,66 +52,52 @@ class Smearing(NumericalSettings): class Mesh(ArchiveSection): """ - A base section used to specify the settings of a sampling mesh. - It supports uniformly-spaced meshes and symmetry-reduced representations. + A base section used to define the mesh or space partitioning over which a discrete numerical integration is performed. """ - spacing = Quantity( - type=MEnum('Equidistant', 'Logarithmic', 'Tan'), - shape=['dimensionality'], + dimensionality = Quantity( + type=np.int32, + default=3, description=""" - Identifier for the spacing of the Mesh. Defaults to 'Equidistant' if not defined. It can take the values: - - | Name | Description | - | --------- | -------------------------------- | - | `'Equidistant'` | Equidistant grid (also known as 'Newton-Cotes') | - | `'Logarithmic'` | log distance grid | - | `'Tan'` | Non-uniform tan mesh for grids. More dense at low abs values of the points, while less dense for higher values | + Dimensionality of the mesh: 1, 2, or 3. Defaults to 3. """, ) - quadrature = Quantity( - type=MEnum( - 'Gauss-Legendre', - 'Gauss-Laguerre', - 'Clenshaw-Curtis', - 'Newton-Cotes', - 'Gauss-Hermite', - ), + mesh_type = Quantity( + type=MEnum('equidistant', 'logarithmic', 'tan'), + shape=['dimensionality'], description=""" - Quadrature rule used for integration of the Mesh. This quantity is relevant for 1D meshes: + Kind of mesh identifying the spacing in each of the dimensions specified by `dimensionality`. It can take the values: | Name | Description | | --------- | -------------------------------- | - | `'Gauss-Legendre'` | Quadrature rule for integration using Legendre polynomials | - | `'Gauss-Laguerre'` | Quadrature rule for integration using Laguerre polynomials | - | `'Clenshaw-Curtis'` | Quadrature rule for integration using Chebyshev polynomials using discrete cosine transformations | - | `'Gauss-Hermite'` | Quadrature rule for integration using Hermite polynomials | + | `'equidistant'` | Equidistant grid (also known as 'Newton-Cotes') | + | `'logarithmic'` | log distance grid | + | `'Tan'` | Non-uniform tan mesh for grids. More dense at low abs values of the points, while less dense for higher values | """, - ) # ! @JosePizarro3 I think that this is separate from the spacing + ) - n_points = Quantity( + grid = Quantity( type=np.int32, + shape=['dimensionality'], description=""" - Number of points in the mesh. + Number of points sampled along each axis of the mesh. """, ) - dimensionality = Quantity( + n_points = Quantity( type=np.int32, - default=3, description=""" - Dimensionality of the mesh: 1, 2, or 3. Defaults to 3. + Total number of points in the mesh. """, ) - grid = Quantity( - type=np.int32, + spacing = Quantity( + type=np.float64, shape=['dimensionality'], - description=""" - Amount of mesh point sampling along each axis. See `type` for the axes definition. + description="""Grid spacing for equidistant meshes. Ignored for other kinds. """, - ) # ? @JosePizzaro3: should the mesh also contain its boundary information + ) points = Quantity( type=np.complex128, @@ -126,25 +112,102 @@ class Mesh(ArchiveSection): shape=['n_points'], description=""" The amount of times the same point reappears. A value larger than 1, typically indicates - a symmetry operation that was applied to the `Mesh`. This quantity is equivalent to `weights`: + a symmetry operation that was applied to the `Mesh`. + """, + ) - multiplicities = n_points * weights + pruning = Quantity( + type=MEnum('fixed', 'adaptive'), + description=""" + Pruning method applied for reducing the amount of points in the Mesh. This is typically + used for numerical integration near the core levels in atoms. + In the fixed grid methods, the number of angular grid points is predetermined for + ranges of radial grid points, while in the adaptive methods, the angular grid is adjusted + on-the-fly for each radial point according to some accuracy criterion. """, ) - weights = Quantity( + def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: + super().normalize(archive, logger) + + if self.dimensionality not in [1, 2, 3]: + logger.error( + '`dimensionality` meshes different than 1, 2, or 3 are not supported.' + ) + + +class NumericalIntegration(NumericalSettings): + """ + Numerical integration settings used to resolve the following type of integrals by discrete + numerical integration: + + ```math + \int_{\vec{r}_a}^{\vec{r}_b} d^3 \vec{r} F(\vec{r}) \approx \sum_{n=a}^{b} w(\vec{r}_n) F(\vec{r}_n) + ``` + + Here, $F$ can be any type of function which would define the type of rules that can be applied + to solve such integral (e.g., 1D Gaussian quadrature rule or multi-dimensional `angular` rules like the + Lebedev quadrature rule). + + These multidimensional integral has a `Mesh` defined over which the integration is performed, i.e., the + $\vec{r}_n$ points. + """ + + mesh = SubSection(sub_section=Mesh.m_def) + + coordinate = Quantity( + type=MEnum('full', 'radial', 'angular'), + description=""" + Coordinate over which the integration is performed. `full` means the integration is performed in + entire space. `radial` and `angular` describe cases where the integration is performed for + functions which can be splitted into radial and angular distributions (e.g., orbital wavefunctions). + """, + ) + + integration_rule = Quantity( + type=str, # ? extend to MEnum? + description=""" + Integration rule used. This can be any 1D Gaussian quadrature rule or multi-dimensional `angular` rules, + e.g., Lebedev quadrature rule (see e.g., Becke, Chem. Phys. 88, 2547 (1988)). + """, + ) + + integration_thresh = Quantity( type=np.float64, - shape=['n_points'], description=""" - Weight of each point. A value smaller than 1, typically indicates a symmetry operation that was - applied to the mesh. This quantity is equivalent to `multiplicities`: + Accuracy threshold for integration grid. + GRIDTHR in Molpro. + BFCut in ORCA. + """, + ) - weights = multiplicities / n_points + weight_approximation = Quantity( + type=str, + description=""" + Approximation applied to the weight when doing the numerical integration. + See e.g., C. W. Murray, N. C. Handy + and G. J. Laming, Mol. Phys. 78, 997 (1993). + """, + ) + + weight_cutoff = Quantity( + type=np.float64, + description=""" + Threshold for discarding small weights during integration. + Grid points very close to the nucleus can have very small grid weights. + WEIGHT_CUT in Molpro. + Wcut in ORCA. """, ) def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: super().normalize(archive, logger) + valid_coordinates = ['full', 'radial', 'angular', None] + if self.coordinate not in valid_coordinates: + logger.warning( + f'Invalid coordinate value: {self.coordinate}. Resetting to None.' + ) + self.coordinate = None class KSpaceFunctionalities: @@ -887,3 +950,33 @@ def __init__(self, m_def: 'Section' = None, m_context: 'Context' = None, **kwarg def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: super().normalize(archive, logger) + + +class OrbitalLocalization(SelfConsistency): + """ + Numerical settings that control orbital localization. + """ + + localization_method = Quantity( + type=MEnum('FB', 'PM', 'IBO', 'IAOIBO', 'IAOBOYS', 'NEWBOYS', 'AHFB'), + description=""" + Name of the localization method. + """, + ) + + orbital_window = Quantity( + shape=['*'], + description=""" + the Molecular orbital range to be localized. + """, + ) + + core_threshold = Quantity( + type=np.float64, + description=""" + the energy window for the first occupied MO to be localized (in a.u.). + """, + ) + + def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: + super().normalize(archive, logger) diff --git a/tests/test_basis_set.py b/tests/test_basis_set.py index b1dcb03c..7e6dd427 100644 --- a/tests/test_basis_set.py +++ b/tests/test_basis_set.py @@ -15,6 +15,7 @@ APWOrbital, APWPlaneWaveBasisSet, AtomCenteredBasisSet, + AtomCenteredFunction, BasisSetContainer, MuffinTinRegion, PlaneWaveBasisSet, @@ -418,3 +419,134 @@ def test_quick_step() -> None: ], } # TODO: generate a QuickStep generator in the CP2K plugin + + +@pytest.mark.parametrize( + 'basis_set_name, basis_type, role, expected_name, expected_type, expected_role', + [ + ('cc-pVTZ', 'GTO', 'orbital', 'cc-pVTZ', 'GTO', 'orbital'), + ('def2-TZVP', 'GTO', 'auxiliary_scf', 'def2-TZVP', 'GTO', 'auxiliary_scf'), + ( + 'aug-cc-pVDZ', + 'STO', + 'auxiliary_post_hf', + 'aug-cc-pVDZ', + 'STO', + 'auxiliary_post_hf', + ), + ( + 'custom_basis', + None, + None, + 'custom_basis', + None, + None, + ), # Undefined type and role + ], +) +def test_atom_centered_basis_set_init( + basis_set_name, basis_type, role, expected_name, expected_type, expected_role +): + """Test initialization of AtomCenteredBasisSet.""" + bs = AtomCenteredBasisSet(basis_set=basis_set_name, type=basis_type, role=role) + assert bs.basis_set == expected_name + assert bs.type == expected_type + assert bs.role == expected_role + + +@pytest.mark.parametrize( + 'functions, expected_count', + [ + ( + [ + AtomCenteredFunction( + harmonic_type='spherical', + function_type='s', + n_primitive=3, + exponents=[1.0, 2.0, 3.0], + contraction_coefficients=[0.5, 0.3, 0.2], + ), + ], + 1, + ), + ( + [ + AtomCenteredFunction( + harmonic_type='cartesian', + function_type='p', + n_primitive=1, + exponents=[0.5], + contraction_coefficients=[1.0], + ), + AtomCenteredFunction( + harmonic_type='spherical', + function_type='d', + n_primitive=2, + exponents=[1.0, 2.0], + contraction_coefficients=[0.4, 0.6], + ), + ], + 2, + ), + ], +) +def test_atom_centered_basis_set_functional_composition(functions, expected_count): + """Test functional composition within AtomCenteredBasisSet.""" + bs = AtomCenteredBasisSet(functional_composition=functions) + assert len(bs.functional_composition) == expected_count + for f, ref_f in zip(bs.functional_composition, functions): + assert f.harmonic_type == ref_f.harmonic_type + assert f.function_type == ref_f.function_type + assert f.n_primitive == ref_f.n_primitive + assert np.allclose(f.exponents, ref_f.exponents) + assert np.allclose(f.contraction_coefficients, ref_f.contraction_coefficients) + + +def test_atom_centered_basis_set_normalize(): + """Test normalization of AtomCenteredBasisSet.""" + bs = AtomCenteredBasisSet( + basis_set='cc-pVTZ', + type='GTO', + role='orbital', + functional_composition=[ + AtomCenteredFunction( + harmonic_type='spherical', + function_type='s', + n_primitive=2, + exponents=[1.0, 2.0], + contraction_coefficients=[0.5, 0.5], + ) + ], + ) + bs.normalize(None, None) + # Add assertions for normalized attributes if needed + assert bs.basis_set == 'cc-pVTZ' + assert bs.type == 'GTO' + assert bs.role == 'orbital' + assert len(bs.functional_composition) == 1 + + +def test_atom_centered_basis_set_invalid_data(): + """Test behavior with missing or invalid data.""" + bs = AtomCenteredBasisSet( + basis_set='invalid_basis', + type=None, # Missing type + role=None, # Missing role + ) + assert bs.basis_set == 'invalid_basis' + assert bs.type is None + assert bs.role is None + + # Test functional composition with invalid data + invalid_function = AtomCenteredFunction( + harmonic_type='spherical', + function_type='s', + n_primitive=2, + exponents=[1.0], # Mismatched length + contraction_coefficients=[0.5, 0.5], + ) + bs.functional_composition = [invalid_function] + + # Call normalize to trigger validation + with pytest.raises(ValueError, match='Mismatch in number of exponents'): + invalid_function.normalize(None, None) diff --git a/tests/test_numerical_settings.py b/tests/test_numerical_settings.py index 6426b4cd..0ead803f 100644 --- a/tests/test_numerical_settings.py +++ b/tests/test_numerical_settings.py @@ -9,6 +9,8 @@ KLinePath, KMesh, KSpaceFunctionalities, + Mesh, + NumericalIntegration, ) from . import logger @@ -377,3 +379,82 @@ def test_resolve_points(self, k_line_path: KLinePath): ] ) assert np.allclose(k_line_path.points, points) + + +@pytest.mark.parametrize( + 'dimensionality, expected_warning', + [ + (3, None), # Valid case + (2, None), # Valid case + ( + 5, + '`dimensionality` meshes different than 1, 2, or 3 are not supported.', + ), # Invalid + ( + 0, + '`dimensionality` meshes different than 1, 2, or 3 are not supported.', + ), # Invalid + ], +) +def test_mesh_dimensionality_validation(dimensionality, expected_warning, caplog): + mesh = Mesh(dimensionality=dimensionality) + mesh.normalize(None, logger) + if expected_warning: + assert expected_warning in caplog.text + else: + assert caplog.text == '' + + +@pytest.mark.parametrize( + 'dimensionality, grid, points', + [ + (3, [10, 10, 10], None), # Valid grid, no points defined yet + (3, None, None), # Missing grid and points + ( + 3, + [10, 10, 10], + [[0, 0, 0], [1, 1, 1]], + ), # Valid points (though fewer than grid suggests) + ], +) +def test_mesh_grid_and_points(dimensionality, grid, points): + mesh = Mesh(dimensionality=dimensionality, grid=grid, points=points) + assert mesh.dimensionality == dimensionality + if grid is not None: + assert np.allclose(mesh.grid, grid) + else: + assert mesh.grid == grid + if points is not None: + assert np.allclose(mesh.points, points) + else: + assert mesh.points == points + + +def test_mesh_spacing_normalization(): + mesh = Mesh(dimensionality=3, grid=[10, 10, 10], spacing=[0.1, 0.1, 0.1]) + mesh.normalize(None, logger) + assert np.allclose(mesh.spacing, [0.1, 0.1, 0.1]) + + +def test_numerical_integration_mesh(): + mesh = Mesh(dimensionality=3, grid=[10, 10, 10]) + integration = NumericalIntegration(mesh=mesh) + assert integration.mesh.dimensionality == 3 + assert np.allclose(integration.mesh.grid, [10, 10, 10]) + + +@pytest.mark.parametrize( + 'integration_thresh, weight_cutoff', + [ + (1e-6, 1e-3), # Valid thresholds + (None, 1e-3), # Missing integration threshold + (1e-6, None), # Missing weight cutoff + (None, None), # Both thresholds missing + ], +) +def test_numerical_integration_thresholds(integration_thresh, weight_cutoff): + integration = NumericalIntegration( + integration_thresh=integration_thresh, weight_cutoff=weight_cutoff + ) + assert integration.integration_thresh == integration_thresh + assert integration.weight_cutoff == weight_cutoff