From d436e95deffe59272913fb561525f37ca1c750dc Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 24 Sep 2024 14:26:59 +0200
Subject: [PATCH 01/37] add main and auxiliary basis set quantities to
 AtomCenteredBasisSet

---
 .../schema_packages/basis_set.py              | 30 ++++++++++++++++++-
 1 file changed, 29 insertions(+), 1 deletion(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index bbd763e0..cea559ec 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -200,13 +200,41 @@ class AtomCenteredBasisSet(BasisSetComponent):
     Defines an atom-centered basis set.
     """
 
+    main_basis_set = Quantity(
+        type=str,
+        description="""
+        Name of the main basis set.
+        """,
+    )
+
+    aux_c_basis_set = Quantity(
+        type=str,
+        description="""
+        AuxC type of basis set.
+        """
+    )
+
+    aux_j_basis_set = Quantity(
+        type=str,
+        description="""
+        AuxJ type of basis set.
+        """
+    )
+
+    aux_jk_basis_set = Quantity(
+        type=str,
+        description="""
+        AuxJK type of basis set.
+        """
+    )
+
     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
+        self.name = self.m_def.name
         # TODO: set name based on basis functions
         # ? use basis set names from Basis Set Exchange
 

From 022c1fadba9188bdba0dccc0621a8fabb170dd5e Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 24 Sep 2024 14:32:22 +0200
Subject: [PATCH 02/37] formatting

---
 src/nomad_simulations/schema_packages/basis_set.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index cea559ec..90a97408 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -211,21 +211,21 @@ class AtomCenteredBasisSet(BasisSetComponent):
         type=str,
         description="""
         AuxC type of basis set.
-        """
+        """,
     )
 
     aux_j_basis_set = Quantity(
         type=str,
         description="""
         AuxJ type of basis set.
-        """
+        """,
     )
 
     aux_jk_basis_set = Quantity(
         type=str,
         description="""
         AuxJK type of basis set.
-        """
+        """,
     )
 
     functional_composition = SubSection(

From af88624e32f93fc46dba294a4b4284e73e5f9519 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 25 Sep 2024 14:04:34 +0200
Subject: [PATCH 03/37] add a draft for AtomCenteredFunction

---
 .../schema_packages/basis_set.py              | 35 +++++++++++++++++--
 1 file changed, 33 insertions(+), 2 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 90a97408..984326ef 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -190,9 +190,40 @@ class AtomCenteredFunction(ArchiveSection):
     Specifies a single function (term) in an atom-centered basis set.
     """
 
-    pass
+    function_type = Quantity(
+        type=str,
+        description="""
+        Type of the function (e.g. GTO for Gaussian, STO for Slater)
+        """,
+    )
+
+    exponents = Quantity(
+        type=np.float32,
+        shape=['*'],
+        description="""
+        List of exponents for the basis function.
+        """,
+    )
 
-    # TODO: design system for writing basis functions like gaussian or slater orbitals
+    contraction_coefficients = Quantity(
+        type=np.float32,
+        shape=['*'],
+        description="""
+        List of contraction coefficients corresponding to the exponents.
+        """,
+    )
+
+    atom_state = SubSection(sub_section=AtomsState.m_def, repeats=False)
+
+    def __init__(self, atom_state: AtomsState, function_type: str, exponents: list, contraction_coefficients: list):
+        self.atom_state = atom_state
+        self.function_type = function_type
+        self.exponents = exponents
+        self.contraction_coefficients = contraction_coefficients
+
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        super().normalize(archive, logger)
+        self.name = self.m_def.name
 
 
 class AtomCenteredBasisSet(BasisSetComponent):

From 0a9d98b19e2a154d7922ce4530bce10ef594a746 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 25 Sep 2024 15:09:10 +0200
Subject: [PATCH 04/37] add a quantity for number of primitives

---
 .../schema_packages/basis_set.py              | 25 +++++++++++++++----
 1 file changed, 20 insertions(+), 5 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 984326ef..ea165072 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -191,9 +191,16 @@ class AtomCenteredFunction(ArchiveSection):
     """
 
     function_type = Quantity(
-        type=str,
+        type=MEnum('S', 'P', 'D', 'F', 'G', 'H', 'I', 'J'),
+        description="""
+        the angular momentum of the shell to be added.
+        """,
+    )
+
+    n_primitive = Quantity(
+        type=int,
         description="""
-        Type of the function (e.g. GTO for Gaussian, STO for Slater)
+        Number of primitives.
         """,
     )
 
@@ -215,15 +222,23 @@ class AtomCenteredFunction(ArchiveSection):
 
     atom_state = SubSection(sub_section=AtomsState.m_def, repeats=False)
 
-    def __init__(self, atom_state: AtomsState, function_type: str, exponents: list, contraction_coefficients: list):
+    def __init__(
+        self,
+        atom_state: AtomsState,
+        function_type: str,
+        n_primitive: int,
+        exponents: list,
+        contraction_coefficients: list,
+    ):
         self.atom_state = atom_state
         self.function_type = function_type
+        self.n_primitive = n_primitive
         self.exponents = exponents
         self.contraction_coefficients = contraction_coefficients
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
-        self.name = self.m_def.name
+        # self.name = self.m_def.name
 
 
 class AtomCenteredBasisSet(BasisSetComponent):
@@ -265,7 +280,7 @@ class AtomCenteredBasisSet(BasisSetComponent):
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
-        self.name = self.m_def.name
+        # self.name = self.m_def.name
         # TODO: set name based on basis functions
         # ? use basis set names from Basis Set Exchange
 

From a204ec481e232b20c7ee334e7669bf0804d8d3f1 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 26 Sep 2024 13:10:18 +0200
Subject: [PATCH 05/37] add atoms_state reference to AtomCenteredBasisSet

---
 .../schema_packages/basis_set.py              | 28 +++++++------------
 1 file changed, 10 insertions(+), 18 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index ea165072..0941a628 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -206,7 +206,7 @@ class AtomCenteredFunction(ArchiveSection):
 
     exponents = Quantity(
         type=np.float32,
-        shape=['*'],
+        shape=['n_primitive'],
         description="""
         List of exponents for the basis function.
         """,
@@ -214,28 +214,12 @@ class AtomCenteredFunction(ArchiveSection):
 
     contraction_coefficients = Quantity(
         type=np.float32,
-        shape=['*'],
+        shape=['n_primitive'],
         description="""
         List of contraction coefficients corresponding to the exponents.
         """,
     )
 
-    atom_state = SubSection(sub_section=AtomsState.m_def, repeats=False)
-
-    def __init__(
-        self,
-        atom_state: AtomsState,
-        function_type: str,
-        n_primitive: int,
-        exponents: list,
-        contraction_coefficients: list,
-    ):
-        self.atom_state = atom_state
-        self.function_type = function_type
-        self.n_primitive = n_primitive
-        self.exponents = exponents
-        self.contraction_coefficients = contraction_coefficients
-
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
         # self.name = self.m_def.name
@@ -274,6 +258,14 @@ class AtomCenteredBasisSet(BasisSetComponent):
         """,
     )
 
+    atoms_ref = Quantity(
+        type=AtomsState,
+        shape=['*'],
+        description="""
+        References to the `AtomsState` sections that define the atoms this basis set applies to.
+        """,
+    )
+
     functional_composition = SubSection(
         sub_section=AtomCenteredFunction.m_def, repeats=True
     )  # TODO change name

From 8866ab7785a103cf11c50698f4fccd9637adf6d1 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 2 Oct 2024 11:03:52 +0200
Subject: [PATCH 06/37] assign JSON format to basis set

---
 .../schema_packages/basis_set.py              | 46 +++++++++++++++++--
 1 file changed, 43 insertions(+), 3 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 0941a628..97b01c4a 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 MEnum, Quantity, SubSection, JSON
 from nomad.units import ureg
 
 from nomad_simulations.schema_packages.atoms_state import AtomsState
@@ -191,7 +191,7 @@ class AtomCenteredFunction(ArchiveSection):
     """
 
     function_type = Quantity(
-        type=MEnum('S', 'P', 'D', 'F', 'G', 'H', 'I', 'J'),
+        type=MEnum('s', 'p', 'd', 'f', 'g', 'h', 'i', 'j'),
         description="""
         the angular momentum of the shell to be added.
         """,
@@ -224,7 +224,47 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
         # self.name = self.m_def.name
 
+class AtomCenteredBasisSet(BasisSetComponent):
+    """
+    Defines an atom-centered basis set, using a single JSON quantity for all basis set types,
+    while allowing different AtomState references for each basis set type.
+    """
+
+    basis_set_data = Quantity(
+        type=JSON,  # Use JSON to store basis set information, including atom references
+        description="""
+        JSON object containing all the basis set information along with atom references. Example:
+        {
+            "main_basis_set": {
+                "name": "cc-pVTZ",
+                "atoms_ref": [ref_to_atoms_1, ref_to_atoms_2]
+            },
+            "aux_c_basis_set": {
+                "name": "cc-pVTZ/C",
+                "atoms_ref": [ref_to_atoms_3]
+            },
+            "aux_j_basis_set": {
+                "name": "RIJ",
+                "atoms_ref": [ref_to_atoms_4]
+            },
+            "aux_jk_basis_set": {
+                "name": "aug-cc-pVTZ/JK",
+                "atoms_ref": [ref_to_atoms_1, ref_to_atoms_5]
+            }
+        }
+        """,
+    )
+
+    functional_composition = SubSection(
+        sub_section=AtomCenteredFunction.m_def, repeats=True
+    )
 
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        super().normalize(archive, logger)
+
+
+
+'''
 class AtomCenteredBasisSet(BasisSetComponent):
     """
     Defines an atom-centered basis set.
@@ -275,7 +315,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         # self.name = self.m_def.name
         # TODO: set name based on basis functions
         # ? use basis set names from Basis Set Exchange
-
+'''
 
 class APWBaseOrbital(ArchiveSection):
     """

From c7b078a7e589a4ea1db8c779114ac21dadf4f42f Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 8 Oct 2024 11:38:24 +0200
Subject: [PATCH 07/37] add type and auxiliary_type quantities to
 AtomCenteredBasisSet

---
 .../schema_packages/basis_set.py              | 92 +++++--------------
 1 file changed, 21 insertions(+), 71 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 97b01c4a..c3b8c517 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, JSON
+from nomad.metainfo import JSON, MEnum, Quantity, SubSection
 from nomad.units import ureg
 
 from nomad_simulations.schema_packages.atoms_state import AtomsState
@@ -187,13 +187,16 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 class AtomCenteredFunction(ArchiveSection):
     """
-    Specifies a single function (term) in an atom-centered basis set.
+    Specifies a single function (term) in a Gaussian-type basis set.
+    Cartesian Gaussian-type orbitals (GTOs)
     """
 
     function_type = Quantity(
         type=MEnum('s', 'p', 'd', 'f', 'g', 'h', 'i', 'j'),
         description="""
-        the angular momentum of the shell to be added.
+        the l value:
+        l = i + j + k
+        the angular momentum of GTO to be added.
         """,
     )
 
@@ -224,98 +227,45 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
         # self.name = self.m_def.name
 
-class AtomCenteredBasisSet(BasisSetComponent):
-    """
-    Defines an atom-centered basis set, using a single JSON quantity for all basis set types,
-    while allowing different AtomState references for each basis set type.
-    """
-
-    basis_set_data = Quantity(
-        type=JSON,  # Use JSON to store basis set information, including atom references
-        description="""
-        JSON object containing all the basis set information along with atom references. Example:
-        {
-            "main_basis_set": {
-                "name": "cc-pVTZ",
-                "atoms_ref": [ref_to_atoms_1, ref_to_atoms_2]
-            },
-            "aux_c_basis_set": {
-                "name": "cc-pVTZ/C",
-                "atoms_ref": [ref_to_atoms_3]
-            },
-            "aux_j_basis_set": {
-                "name": "RIJ",
-                "atoms_ref": [ref_to_atoms_4]
-            },
-            "aux_jk_basis_set": {
-                "name": "aug-cc-pVTZ/JK",
-                "atoms_ref": [ref_to_atoms_1, ref_to_atoms_5]
-            }
-        }
-        """,
-    )
-
-    functional_composition = SubSection(
-        sub_section=AtomCenteredFunction.m_def, repeats=True
-    )
-
-    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
-
-
-
-'''
 class AtomCenteredBasisSet(BasisSetComponent):
     """
     Defines an atom-centered basis set.
     """
 
-    main_basis_set = Quantity(
+    basis_set = Quantity(
         type=str,
         description="""
-        Name of the main basis set.
+        name of the basis set
         """,
     )
 
-    aux_c_basis_set = Quantity(
-        type=str,
+    type = Quantity(
+        type=MEnum('STO', 'GTO'),
         description="""
-        AuxC type of basis set.
+        Type of the basis set, e.g. STO or GTO.
         """,
     )
 
-    aux_j_basis_set = Quantity(
-        type=str,
-        description="""
-        AuxJ type of basis set.
-        """,
-    )
+    # TODO: connect RI approximation
 
-    aux_jk_basis_set = Quantity(
-        type=str,
+    auxiliary_type =  Quantity(
+        type=MEnum('AuxC', 'AuxJ', 'AuxJK'),
         description="""
-        AuxJK type of basis set.
-        """,
-    )
-
-    atoms_ref = Quantity(
-        type=AtomsState,
-        shape=['*'],
-        description="""
-        References to the `AtomsState` sections that define the atoms this basis set applies to.
+        the type of RI approximation.
+        AuxJ and AuxJK: Fock matrix construction.
+        AuxC: all other integral generation steps, e.g. post-HF methods.
+        Since a JK-type basis set can be assined to AuxJ, this quantity is needed.
         """,
     )
 
     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
-'''
+        #self.name = self.m_def.name
+
 
 class APWBaseOrbital(ArchiveSection):
     """

From 1208e88c4088d57fe02d8d7ac0fcbf33e2adadaa Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 8 Oct 2024 11:42:54 +0200
Subject: [PATCH 08/37] reformatted basis_set.py

---
 src/nomad_simulations/schema_packages/basis_set.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index c3b8c517..3deaafb4 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -227,6 +227,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
         # self.name = self.m_def.name
 
+
 class AtomCenteredBasisSet(BasisSetComponent):
     """
     Defines an atom-centered basis set.
@@ -248,7 +249,7 @@ class AtomCenteredBasisSet(BasisSetComponent):
 
     # TODO: connect RI approximation
 
-    auxiliary_type =  Quantity(
+    auxiliary_type = Quantity(
         type=MEnum('AuxC', 'AuxJ', 'AuxJK'),
         description="""
         the type of RI approximation.
@@ -264,7 +265,7 @@ class AtomCenteredBasisSet(BasisSetComponent):
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
-        #self.name = self.m_def.name
+        # self.name = self.m_def.name
 
 
 class APWBaseOrbital(ArchiveSection):

From fbc13f9ca069c4e19079102ac92d19d05271b717 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 8 Oct 2024 13:57:01 +0200
Subject: [PATCH 09/37] add NAO and point charges to basis set types

---
 .../schema_packages/basis_set.py              | 38 ++++++++++++++-----
 1 file changed, 28 insertions(+), 10 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 3deaafb4..58c1c40e 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -191,6 +191,19 @@ class AtomCenteredFunction(ArchiveSection):
     Cartesian Gaussian-type orbitals (GTOs)
     """
 
+    # TODO: add a quantity for spherical-harmonic or Cartesian angular functions
+    # Most of the codes use only spherical harmonic.
+
+    basis_type = Quantity(
+        type=MEnum(
+            'spherical',
+            'cartesian',
+        ),
+        description="""
+        spherical-harmonic or cartesian functions.
+        """,
+    )
+
     function_type = Quantity(
         type=MEnum('s', 'p', 'd', 'f', 'g', 'h', 'i', 'j'),
         description="""
@@ -236,26 +249,31 @@ class AtomCenteredBasisSet(BasisSetComponent):
     basis_set = Quantity(
         type=str,
         description="""
-        name of the basis set
+        name of the basis set.
         """,
     )
 
     type = Quantity(
-        type=MEnum('STO', 'GTO'),
+        type=MEnum(
+            'STO',  # Slater-type orbitals
+            'GTO',  # Gaussian-type orbitals
+            'NAO',  # Numerical atomic orbitals
+            'PC',  # Point charges
+        ),
         description="""
         Type of the basis set, e.g. STO or GTO.
         """,
     )
 
-    # TODO: connect RI approximation
-
-    auxiliary_type = Quantity(
-        type=MEnum('AuxC', 'AuxJ', 'AuxJK'),
+    role = Quantity(
+        type=MEnum(
+            'orbital',
+            'auxiliary_scf',
+            'auxiliary_post_hf',
+            'cabs',
+        ),
         description="""
-        the type of RI approximation.
-        AuxJ and AuxJK: Fock matrix construction.
-        AuxC: all other integral generation steps, e.g. post-HF methods.
-        Since a JK-type basis set can be assined to AuxJ, this quantity is needed.
+        role of the basis set
         """,
     )
 

From ddab789f0ce91cca39348e6fb841e85d02ed13f5 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 9 Oct 2024 15:51:53 +0200
Subject: [PATCH 10/37] add cECPs and pointcharges to AtomCenteredBasisSet

---
 .../schema_packages/basis_set.py                | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 58c1c40e..f2e68a8e 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -199,8 +199,9 @@ class AtomCenteredFunction(ArchiveSection):
             'spherical',
             'cartesian',
         ),
+        default='spherical',
         description="""
-        spherical-harmonic or cartesian functions.
+        Specifies whether the basis functions are spherical-harmonic or cartesian functions.
         """,
     )
 
@@ -214,7 +215,7 @@ class AtomCenteredFunction(ArchiveSection):
     )
 
     n_primitive = Quantity(
-        type=int,
+        type=np.int32,
         description="""
         Number of primitives.
         """,
@@ -236,6 +237,13 @@ class AtomCenteredFunction(ArchiveSection):
         """,
     )
 
+    point_charge = Quantity(
+        type=np.int32,
+        description="""
+        the value of the point charge.
+        """,
+    )
+
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
         # self.name = self.m_def.name
@@ -258,6 +266,7 @@ class AtomCenteredBasisSet(BasisSetComponent):
             'STO',  # Slater-type orbitals
             'GTO',  # Gaussian-type orbitals
             'NAO',  # Numerical atomic orbitals
+            'cECP',  # Capped effective core potentials
             'PC',  # Point charges
         ),
         description="""
@@ -270,10 +279,10 @@ class AtomCenteredBasisSet(BasisSetComponent):
             'orbital',
             'auxiliary_scf',
             'auxiliary_post_hf',
-            'cabs',
+            'cabs',  # complementary auxiliary basis set
         ),
         description="""
-        role of the basis set
+        The role of the basis set.
         """,
     )
 

From 8c6ae1a4b3bac5967c1f24e6f4ffa34e24b0179b Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Fri, 11 Oct 2024 14:55:26 +0200
Subject: [PATCH 11/37] fix point charge Quantity type

---
 src/nomad_simulations/schema_packages/basis_set.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index f2e68a8e..39079686 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -238,7 +238,7 @@ class AtomCenteredFunction(ArchiveSection):
     )
 
     point_charge = Quantity(
-        type=np.int32,
+        type=np.float32,
         description="""
         the value of the point charge.
         """,

From c94e995fcf2a00445702d9552b3f1118e752b998 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Mon, 14 Oct 2024 11:08:39 +0200
Subject: [PATCH 12/37] a bit of a cleanup

---
 src/nomad_simulations/schema_packages/basis_set.py | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 39079686..84c4c261 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -187,13 +187,9 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 class AtomCenteredFunction(ArchiveSection):
     """
-    Specifies a single function (term) in a Gaussian-type basis set.
-    Cartesian Gaussian-type orbitals (GTOs)
+    Specifies a single function (term) in an atom-centered basis set.
     """
 
-    # TODO: add a quantity for spherical-harmonic or Cartesian angular functions
-    # Most of the codes use only spherical harmonic.
-
     basis_type = Quantity(
         type=MEnum(
             'spherical',
@@ -206,11 +202,10 @@ class AtomCenteredFunction(ArchiveSection):
     )
 
     function_type = Quantity(
-        type=MEnum('s', 'p', 'd', 'f', 'g', 'h', 'i', 'j'),
+        type=MEnum('s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l'),
         description="""
-        the l value:
-        l = i + j + k
-        the angular momentum of GTO to be added.
+        L=a+b+c
+        The angular momentum of GTO to be added.
         """,
     )
 
@@ -218,6 +213,7 @@ class AtomCenteredFunction(ArchiveSection):
         type=np.int32,
         description="""
         Number of primitives.
+        Linear combinations of the primitive Gaussians are formed to approximate the radial extent of an STO.
         """,
     )
 

From 62de2433e45df887df9ac6308c11623b704ba918 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 24 Oct 2024 17:02:55 +0200
Subject: [PATCH 13/37] move GTOIntegralDecomposition to NumericalSettings

---
 .../schema_packages/numerical_settings.py     | 22 +++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 515deea7..d8d2f781 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -887,3 +887,25 @@ 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 GTOIntegralDecomposition(NumericalSettings):
+    """
+    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,
+        """,
+    )
+
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        super().normalize(archive, logger)
\ No newline at end of file

From 52d6b9ad6da931edc72d597b7f192bc5b540cedd Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 19 Nov 2024 09:24:53 +0100
Subject: [PATCH 14/37] merge develop

---
 src/nomad_simulations/schema_packages/numerical_settings.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index d8d2f781..fb005ff1 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -908,4 +908,4 @@ class GTOIntegralDecomposition(NumericalSettings):
     )
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
\ No newline at end of file
+        super().normalize(archive, logger)

From 646430cf67257fd980047b1715232da571720622 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 19 Nov 2024 17:11:30 +0100
Subject: [PATCH 15/37] add Mesh, NumericalIntegration and
 MolecularHamiltonianSubTerms

---
 .../schema_packages/model_method.py           |  13 ++
 .../schema_packages/numerical_settings.py     | 127 +++++++++++-------
 2 files changed, 91 insertions(+), 49 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index c7b143ff..cfdc487c 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1221,3 +1221,16 @@ class DMFT(ModelMethodElectronic):
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
+
+
+class MolecularHamiltonianSubTerms(BaseModelMethod):
+    type=Quantity(
+        type=MEnum('coulomb', 'exchange'),
+        description="""
+        Typical sub-terms of the molecular hamiltonian.
+        Relativistic effects, SOC, .....
+        """,
+    )
+
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        super().normalize(archive, logger)
\ No newline at end of file
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index fb005ff1..6a5b63a2 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -52,66 +52,45 @@ 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',
-        ),
+    kind = 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 |
-        """,
-    )  # ! @JosePizarro3 I think that this is separate from the spacing
-
-    n_points = Quantity(
-        type=np.int32,
-        description="""
-        Number of points in the mesh.
+        | `'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 = Quantity(
+    grid = Quantity(
         type=np.int32,
-        default=3,
+        shape=['dimensionality'],
         description="""
-        Dimensionality of the mesh: 1, 2, or 3. Defaults to 3.
+        Amount of mesh point sampling along each axis.
         """,
     )
 
-    grid = Quantity(
+    n_points = Quantity(
         type=np.int32,
-        shape=['dimensionality'],
         description="""
-        Amount of mesh point sampling along each axis. See `type` for the axes definition.
+        Number of points in the mesh.
         """,
-    )  # ? @JosePizzaro3: should the mesh also contain its boundary information
+    )
 
     points = Quantity(
         type=np.complex128,
@@ -126,23 +105,73 @@ 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`:
-
-            multiplicities = n_points * weights
+        a symmetry operation that was applied to the `Mesh`.
         """,
     )
 
-    weights = Quantity(
-        type=np.float64,
-        shape=['n_points'],
+    pruning = Quantity(
+        type=MEnum('fixed', 'adaptive'),
         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`:
+        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.
+        """
+    )
+
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        super().normalize(archive, logger)
 
-            weights = multiplicities / n_points
+        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.
+    """
+
+    coordinate = Quantity(
+        type=MEnum('all', 'radial', 'angular'),
+        description="""
+        Coordinate over which the integration is performed. `all` means the integration is performed in
+        all the 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)).
+        """
+    )
+
+    weight_partitioning = 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).
+        """
+    )
+
+    mesh = SubSection(sub_section=Mesh.m_def)
+
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
@@ -908,4 +937,4 @@ class GTOIntegralDecomposition(NumericalSettings):
     )
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
+        super().normalize(archive, logger)
\ No newline at end of file

From cca109f6ae371e0cc882b2709578ac086c68b836 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 19 Nov 2024 17:27:06 +0100
Subject: [PATCH 16/37] minor adjustments to Mesh and NumericalIntegration

---
 .../schema_packages/model_method.py           |  4 +--
 .../schema_packages/numerical_settings.py     | 26 ++++++++++++-------
 2 files changed, 19 insertions(+), 11 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index cfdc487c..9c1c727f 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1224,7 +1224,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 
 class MolecularHamiltonianSubTerms(BaseModelMethod):
-    type=Quantity(
+    type = Quantity(
         type=MEnum('coulomb', 'exchange'),
         description="""
         Typical sub-terms of the molecular hamiltonian.
@@ -1233,4 +1233,4 @@ class MolecularHamiltonianSubTerms(BaseModelMethod):
     )
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
\ No newline at end of file
+        super().normalize(archive, logger)
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 6a5b63a2..2986b674 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -92,6 +92,12 @@ class Mesh(ArchiveSection):
         """,
     )
 
+    spacing = Quantity(
+        type=np.float64,
+        shape=['dimensionality'],
+        description='Grid spacing for equidistant meshes. Ignored for other kinds.',
+    )
+
     points = Quantity(
         type=np.complex128,
         shape=['n_points', 'dimensionality'],
@@ -117,14 +123,16 @@ class Mesh(ArchiveSection):
         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.
-        """
+        """,
     )
 
     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.')
+            logger.error(
+                '`dimensionality` meshes different than 1, 2, or 3 are not supported.'
+            )
 
 
 class NumericalIntegration(NumericalSettings):
@@ -145,10 +153,10 @@ class NumericalIntegration(NumericalSettings):
     """
 
     coordinate = Quantity(
-        type=MEnum('all', 'radial', 'angular'),
+        type=MEnum('full', 'radial', 'angular'),
         description="""
-        Coordinate over which the integration is performed. `all` means the integration is performed in
-        all the space. `radial` and `angular` describe cases where the integration is performed for
+        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).
         """,
     )
@@ -158,16 +166,16 @@ class NumericalIntegration(NumericalSettings):
         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)).
-        """
+        """,
     )
 
-    weight_partitioning = Quantity(
+    weight_approximtion = 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).
-        """
+        """,
     )
 
     mesh = SubSection(sub_section=Mesh.m_def)
@@ -937,4 +945,4 @@ class GTOIntegralDecomposition(NumericalSettings):
     )
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
\ No newline at end of file
+        super().normalize(archive, logger)

From 42f4f38ac18c7079f6438ffc54afae94a248452b Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 20 Nov 2024 11:42:40 +0100
Subject: [PATCH 17/37] add integration_thresh and weight_cutoff to
 NumericalIntegration

---
 .../schema_packages/numerical_settings.py     | 22 ++++++++++++++++++-
 1 file changed, 21 insertions(+), 1 deletion(-)

diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 2986b674..63fb2ab4 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -169,7 +169,16 @@ class NumericalIntegration(NumericalSettings):
         """,
     )
 
-    weight_approximtion = Quantity(
+    integration_thresh = Quantity(
+        type=np.float64,
+        description="""
+        Accuracy threshold for integration grid.
+        GRIDTHR in Molpro.
+        BFCut in ORCA.
+        """,
+    )
+
+    weight_approximation = Quantity(
         type=str,
         description="""
         Approximation applied to the weight when doing the numerical integration.
@@ -178,6 +187,17 @@ class NumericalIntegration(NumericalSettings):
         """,
     )
 
+    weight_cutoff = Quantity(
+        type=np.float64,
+        description="""
+        Grid points very close to the nucleus can have very small grid weights.
+        These can be discarded with the option WEIGHT_CUT=thr, i.e., grid points with weights
+        smaller than thr will then not be used in the numerical integration anymore.
+        WEIGHT_CUT in Molpro.
+        Wcut in ORCA.
+        """,
+    )
+
     mesh = SubSection(sub_section=Mesh.m_def)
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:

From 0f92eb9ea0468f7417f30fd0edc0a2f76abfdbdc Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 20 Nov 2024 13:32:31 +0100
Subject: [PATCH 18/37] check whether n_primitive matches the lengths of
 exponents and contraction coefficients

---
 src/nomad_simulations/schema_packages/basis_set.py | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 84c4c261..5fd5ac03 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -242,7 +242,19 @@ class AtomCenteredFunction(ArchiveSection):
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
-        # self.name = self.m_def.name
+
+        # Validation: Check that n_primitive matches the lengths of exponents and contraction coefficients
+        if self.n_primitive is not None:
+            if len(self.exponents or []) != self.n_primitive:
+                logger.error(
+                    f"Mismatch in number of exponents: expected {self.n_primitive}, "
+                    f"found {len(self.exponents or [])}."
+                )
+            if len(self.contraction_coefficients or []) != self.n_primitive:
+                logger.error(
+                    f"Mismatch in number of contraction coefficients: expected {self.n_primitive}, "
+                    f"found {len(self.contraction_coefficients or [])}."
+                )
 
 
 class AtomCenteredBasisSet(BasisSetComponent):

From e8fb5ef559c83170c1e0facedc521104f0892fc5 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 20 Nov 2024 13:48:53 +0100
Subject: [PATCH 19/37] add tests for AtomCenteredBasisSet and
 AtomCenteredFunction

---
 .../schema_packages/basis_set.py              |  14 +--
 tests/test_basis_set.py                       | 109 ++++++++++++++++++
 2 files changed, 115 insertions(+), 8 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 5fd5ac03..4aaf55f6 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -245,15 +245,13 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
         # Validation: Check that n_primitive matches the lengths of exponents and contraction coefficients
         if self.n_primitive is not None:
-            if len(self.exponents or []) != self.n_primitive:
-                logger.error(
-                    f"Mismatch in number of exponents: expected {self.n_primitive}, "
-                    f"found {len(self.exponents or [])}."
+            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}, found {len(self.exponents)}."
                 )
-            if len(self.contraction_coefficients or []) != self.n_primitive:
-                logger.error(
-                    f"Mismatch in number of contraction coefficients: expected {self.n_primitive}, "
-                    f"found {len(self.contraction_coefficients or [])}."
+            if self.contraction_coefficients is not None and len(self.contraction_coefficients) != self.n_primitive:
+                raise ValueError(
+                f"Mismatch in number of contraction coefficients: expected {self.n_primitive}, found {len(self.contraction_coefficients)}."
                 )
 
 
diff --git a/tests/test_basis_set.py b/tests/test_basis_set.py
index b1dcb03c..e96c224b 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,111 @@ def test_quick_step() -> None:
         ],
     }
     # TODO: generate a QuickStep generator in the CP2K plugin
+
+
+@pytest.mark.parametrize(
+    'basis_set_name, basis_type, role',
+    [
+        ('cc-pVTZ', 'GTO', 'orbital'),
+        ('def2-TZVP', 'GTO', 'auxiliary_scf'),
+        ('aug-cc-pVDZ', 'STO', 'auxiliary_post_hf'),
+        ('custom_basis', None, None),  # Undefined type and role
+    ],
+)
+def test_atom_centered_basis_set_init(basis_set_name, basis_type, role) -> None:
+    """Test initialization of AtomCenteredBasisSet."""
+    bs = AtomCenteredBasisSet(basis_set=basis_set_name, type=basis_type, role=role)
+    assert bs.basis_set == basis_set_name
+    assert bs.type == basis_type
+    assert bs.role == role
+
+
+@pytest.mark.parametrize(
+    'functions',
+    [
+        [
+            AtomCenteredFunction(
+                basis_type='spherical',
+                function_type='s',
+                n_primitive=3,
+                exponents=[1.0, 2.0, 3.0],
+                contraction_coefficients=[0.5, 0.3, 0.2],
+            ),
+        ],
+        [
+            AtomCenteredFunction(
+                basis_type='cartesian',
+                function_type='p',
+                n_primitive=1,
+                exponents=[0.5],
+                contraction_coefficients=[1.0],
+            ),
+            AtomCenteredFunction(
+                basis_type='spherical',
+                function_type='d',
+                n_primitive=2,
+                exponents=[1.0, 2.0],
+                contraction_coefficients=[0.4, 0.6],
+            ),
+        ],
+    ],
+)
+def test_atom_centered_basis_set_functional_composition(functions) -> None:
+    """Test functional composition within AtomCenteredBasisSet."""
+    bs = AtomCenteredBasisSet(functional_composition=functions)
+    assert len(bs.functional_composition) == len(functions)
+    for f, ref_f in zip(bs.functional_composition, functions):
+        assert f.basis_type == ref_f.basis_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() -> None:
+    """Test normalization of AtomCenteredBasisSet."""
+    bs = AtomCenteredBasisSet(
+        basis_set='cc-pVTZ',
+        type='GTO',
+        role='orbital',
+        functional_composition=[
+            AtomCenteredFunction(
+                basis_type='spherical',
+                function_type='s',
+                n_primitive=2,
+                exponents=[1.0, 2.0],
+                contraction_coefficients=[0.5, 0.5],
+            )
+        ],
+    )
+    bs.normalize(None, logger)
+    # Add checks for normalized behavior, if any
+    assert bs.basis_set == 'cc-pVTZ'
+
+def test_atom_centered_basis_set_invalid_data() -> None:
+    """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(
+        basis_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, logger)
+
+
+

From 4826f0c5d3ab84a6722c3652c8eec3a4855a6829 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 20 Nov 2024 14:09:53 +0100
Subject: [PATCH 20/37] add tests for Mesh and NumericalIntegration

---
 .../schema_packages/basis_set.py              |  9 ++-
 .../schema_packages/numerical_settings.py     |  6 ++
 tests/test_basis_set.py                       |  6 +-
 tests/test_numerical_settings.py              | 81 +++++++++++++++++++
 4 files changed, 95 insertions(+), 7 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 4aaf55f6..0428c0c6 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -247,11 +247,14 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         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}, found {len(self.exponents)}."
+                    f'Mismatch in number of exponents: expected {self.n_primitive}, found {len(self.exponents)}.'
                 )
-            if self.contraction_coefficients is not None and len(self.contraction_coefficients) != self.n_primitive:
+            if (
+                self.contraction_coefficients is not None
+                and len(self.contraction_coefficients) != self.n_primitive
+            ):
                 raise ValueError(
-                f"Mismatch in number of contraction coefficients: expected {self.n_primitive}, found {len(self.contraction_coefficients)}."
+                    f'Mismatch in number of contraction coefficients: expected {self.n_primitive}, found {len(self.contraction_coefficients)}.'
                 )
 
 
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 63fb2ab4..37842f84 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -202,6 +202,12 @@ class NumericalIntegration(NumericalSettings):
 
     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:
diff --git a/tests/test_basis_set.py b/tests/test_basis_set.py
index e96c224b..eb2e20f0 100644
--- a/tests/test_basis_set.py
+++ b/tests/test_basis_set.py
@@ -500,6 +500,7 @@ def test_atom_centered_basis_set_normalize() -> None:
     # Add checks for normalized behavior, if any
     assert bs.basis_set == 'cc-pVTZ'
 
+
 def test_atom_centered_basis_set_invalid_data() -> None:
     """Test behavior with missing or invalid data."""
     bs = AtomCenteredBasisSet(
@@ -522,8 +523,5 @@ def test_atom_centered_basis_set_invalid_data() -> None:
     bs.functional_composition = [invalid_function]
 
     # Call normalize to trigger validation
-    with pytest.raises(ValueError, match="Mismatch in number of exponents"):
+    with pytest.raises(ValueError, match='Mismatch in number of exponents'):
         invalid_function.normalize(None, logger)
-
-
-
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

From c5141abd65be9ec2447ce7e70fee42fcb2aaee05 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 20 Nov 2024 16:56:10 +0100
Subject: [PATCH 21/37] modify Mesh

---
 .../schema_packages/numerical_settings.py        | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 37842f84..1a06b8ad 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -63,7 +63,7 @@ class Mesh(ArchiveSection):
         """,
     )
 
-    kind = Quantity(
+    mesh_type = Quantity(
         type=MEnum('equidistant', 'logarithmic', 'tan'),
         shape=['dimensionality'],
         description="""
@@ -81,21 +81,22 @@ class Mesh(ArchiveSection):
         type=np.int32,
         shape=['dimensionality'],
         description="""
-        Amount of mesh point sampling along each axis.
+        Number of points sampled along each axis of the mesh.
         """,
     )
 
     n_points = Quantity(
         type=np.int32,
         description="""
-        Number of points in the mesh.
+        Total number of points in the mesh.
         """,
     )
 
     spacing = Quantity(
         type=np.float64,
         shape=['dimensionality'],
-        description='Grid spacing for equidistant meshes. Ignored for other kinds.',
+        description="""Grid spacing for equidistant meshes. Ignored for other kinds.
+        """,
     )
 
     points = Quantity(
@@ -152,6 +153,8 @@ class NumericalIntegration(NumericalSettings):
     $\vec{r}_n$ points.
     """
 
+    mesh = SubSection(sub_section=Mesh.m_def)
+
     coordinate = Quantity(
         type=MEnum('full', 'radial', 'angular'),
         description="""
@@ -190,16 +193,13 @@ class NumericalIntegration(NumericalSettings):
     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.
-        These can be discarded with the option WEIGHT_CUT=thr, i.e., grid points with weights
-        smaller than thr will then not be used in the numerical integration anymore.
         WEIGHT_CUT in Molpro.
         Wcut in ORCA.
         """,
     )
 
-    mesh = SubSection(sub_section=Mesh.m_def)
-
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
         valid_coordinates = ['full', 'radial', 'angular', None]

From 23d7615dcc6984854d07192005f66990a979cdff Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 21 Nov 2024 11:15:08 +0100
Subject: [PATCH 22/37] MEnum for MolecularHamiltonianContributions

---
 .../schema_packages/model_method.py              | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index 9c1c727f..47d78939 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1223,12 +1223,22 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
 
-class MolecularHamiltonianSubTerms(BaseModelMethod):
+class MolecularHamiltonianContributions(BaseModelMethod):
     type = Quantity(
-        type=MEnum('coulomb', 'exchange'),
+            type=MEnum(
+            'coulomb',
+            'exchange',
+            'spin_orbit_coupling',
+            'scalar_relativistic',
+            'ri_approximation',
+            'cosx_approximation',
+            'density_fitting',
+            'local_correlation',
+            'explicit_correlation',
+        ), # TODO: expand this
         description="""
+        TODO: Name is only a placeholder.
         Typical sub-terms of the molecular hamiltonian.
-        Relativistic effects, SOC, .....
         """,
     )
 

From 9ef13ea5aedcc8f70db3bb0e7e20bfc55c887a70 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 21 Nov 2024 14:57:41 +0100
Subject: [PATCH 23/37] remove contributions

---
 .../schema_packages/model_method.py           | 23 -------------------
 1 file changed, 23 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index 47d78939..c7b143ff 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1221,26 +1221,3 @@ class DMFT(ModelMethodElectronic):
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
-
-
-class MolecularHamiltonianContributions(BaseModelMethod):
-    type = Quantity(
-            type=MEnum(
-            'coulomb',
-            'exchange',
-            'spin_orbit_coupling',
-            'scalar_relativistic',
-            'ri_approximation',
-            'cosx_approximation',
-            'density_fitting',
-            'local_correlation',
-            'explicit_correlation',
-        ), # TODO: expand this
-        description="""
-        TODO: Name is only a placeholder.
-        Typical sub-terms of the molecular hamiltonian.
-        """,
-    )
-
-    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)

From 80c8b642dc3031e0a314ce4b46df8ae50bd99fa0 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 28 Nov 2024 16:44:30 +0100
Subject: [PATCH 24/37] add a normalizer function for the AtomCenteredFunction
 to handle combined orbital types

---
 .../schema_packages/basis_set.py              | 58 +++++++++++++++----
 1 file changed, 48 insertions(+), 10 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 0428c0c6..7ee5ba87 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -190,7 +190,7 @@ class AtomCenteredFunction(ArchiveSection):
     Specifies a single function (term) in an atom-centered basis set.
     """
 
-    basis_type = Quantity(
+    harmonic_type = Quantity(
         type=MEnum(
             'spherical',
             'cartesian',
@@ -202,7 +202,9 @@ class AtomCenteredFunction(ArchiveSection):
     )
 
     function_type = Quantity(
-        type=MEnum('s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l'),
+        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.
@@ -227,7 +229,7 @@ class AtomCenteredFunction(ArchiveSection):
 
     contraction_coefficients = Quantity(
         type=np.float32,
-        shape=['n_primitive'],
+        shape=['*'], # Flexible shape to handle combined types (e.g. SP, SPD..)
         description="""
         List of contraction coefficients corresponding to the exponents.
         """,
@@ -241,20 +243,51 @@ class AtomCenteredFunction(ArchiveSection):
     )
 
     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)
 
-        # Validation: Check that n_primitive matches the lengths of exponents and contraction coefficients
+        # 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}, found {len(self.exponents)}.'
+                    f"Mismatch in number of exponents: expected {self.n_primitive}, "
+                    f"found {len(self.exponents)}."
                 )
-            if (
-                self.contraction_coefficients is not None
-                and len(self.contraction_coefficients) != self.n_primitive
-            ):
+
+        # 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 number of contraction coefficients: expected {self.n_primitive}, found {len(self.contraction_coefficients)}.'
+                    f"Mismatch in contraction coefficients: expected {self.n_primitive}, "
+                    f"found {len(self.contraction_coefficients)}."
                 )
 
 
@@ -295,6 +328,11 @@ class AtomCenteredBasisSet(BasisSetComponent):
         """,
     )
 
+    total_number_of_basis_functions = Quantity(
+        type=np.int32,
+        description="",
+    )
+
     functional_composition = SubSection(
         sub_section=AtomCenteredFunction.m_def, repeats=True
     )

From 4ac10765abe0581df19aa3f97358d0ccecda353a Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 4 Dec 2024 08:59:59 +0100
Subject: [PATCH 25/37] fix test_basis_set.py

---
 tests/test_basis_set.py | 102 ++++++++++++++++++++++------------------
 1 file changed, 57 insertions(+), 45 deletions(-)

diff --git a/tests/test_basis_set.py b/tests/test_basis_set.py
index eb2e20f0..b209e92a 100644
--- a/tests/test_basis_set.py
+++ b/tests/test_basis_set.py
@@ -422,65 +422,73 @@ def test_quick_step() -> None:
 
 
 @pytest.mark.parametrize(
-    'basis_set_name, basis_type, role',
+    'basis_set_name, basis_type, role, expected_name, expected_type, expected_role',
     [
-        ('cc-pVTZ', 'GTO', 'orbital'),
-        ('def2-TZVP', 'GTO', 'auxiliary_scf'),
-        ('aug-cc-pVDZ', 'STO', 'auxiliary_post_hf'),
-        ('custom_basis', None, None),  # Undefined type and 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) -> None:
+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 == basis_set_name
-    assert bs.type == basis_type
-    assert bs.role == role
+    assert bs.basis_set == expected_name
+    assert bs.type == expected_type
+    assert bs.role == expected_role
 
 
 @pytest.mark.parametrize(
-    'functions',
+    'functions, expected_count',
     [
-        [
-            AtomCenteredFunction(
-                basis_type='spherical',
-                function_type='s',
-                n_primitive=3,
-                exponents=[1.0, 2.0, 3.0],
-                contraction_coefficients=[0.5, 0.3, 0.2],
-            ),
-        ],
-        [
-            AtomCenteredFunction(
-                basis_type='cartesian',
-                function_type='p',
-                n_primitive=1,
-                exponents=[0.5],
-                contraction_coefficients=[1.0],
-            ),
-            AtomCenteredFunction(
-                basis_type='spherical',
-                function_type='d',
-                n_primitive=2,
-                exponents=[1.0, 2.0],
-                contraction_coefficients=[0.4, 0.6],
-            ),
-        ],
+        (
+            [
+                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) -> None:
+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) == len(functions)
+    assert len(bs.functional_composition) == expected_count
     for f, ref_f in zip(bs.functional_composition, functions):
-        assert f.basis_type == ref_f.basis_type
+        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() -> None:
+def test_atom_centered_basis_set_normalize():
     """Test normalization of AtomCenteredBasisSet."""
     bs = AtomCenteredBasisSet(
         basis_set='cc-pVTZ',
@@ -488,7 +496,7 @@ def test_atom_centered_basis_set_normalize() -> None:
         role='orbital',
         functional_composition=[
             AtomCenteredFunction(
-                basis_type='spherical',
+                harmonic_type='spherical',
                 function_type='s',
                 n_primitive=2,
                 exponents=[1.0, 2.0],
@@ -496,12 +504,15 @@ def test_atom_centered_basis_set_normalize() -> None:
             )
         ],
     )
-    bs.normalize(None, logger)
-    # Add checks for normalized behavior, if any
+    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() -> None:
+def test_atom_centered_basis_set_invalid_data():
     """Test behavior with missing or invalid data."""
     bs = AtomCenteredBasisSet(
         basis_set='invalid_basis',
@@ -514,7 +525,7 @@ def test_atom_centered_basis_set_invalid_data() -> None:
 
     # Test functional composition with invalid data
     invalid_function = AtomCenteredFunction(
-        basis_type='spherical',
+        harmonic_type='spherical',
         function_type='s',
         n_primitive=2,
         exponents=[1.0],  # Mismatched length
@@ -524,4 +535,5 @@ def test_atom_centered_basis_set_invalid_data() -> None:
 
     # Call normalize to trigger validation
     with pytest.raises(ValueError, match='Mismatch in number of exponents'):
-        invalid_function.normalize(None, logger)
+        invalid_function.normalize(None, None)
+

From 10e782cec730a8d987c8d7d417fac2b8886c752f Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 4 Dec 2024 09:51:05 +0100
Subject: [PATCH 26/37] add OrbitalLocalization to numerical_settings.py

---
 .../schema_packages/basis_set.py              | 40 +++++++++++++------
 .../schema_packages/numerical_settings.py     | 36 +++++++++++++++++
 tests/test_basis_set.py                       | 19 +++++++--
 3 files changed, 79 insertions(+), 16 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 7ee5ba87..b19f6618 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -202,9 +202,21 @@ class AtomCenteredFunction(ArchiveSection):
     )
 
     function_type = Quantity(
-        type=MEnum('s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l',
-                   'sp', 'spd', 'spdf',
-                   ),
+        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.
@@ -229,7 +241,7 @@ class AtomCenteredFunction(ArchiveSection):
 
     contraction_coefficients = Quantity(
         type=np.float32,
-        shape=['*'], # Flexible shape to handle combined types (e.g. SP, SPD..)
+        shape=['*'],  # Flexible shape to handle combined types (e.g. SP, SPD..)
         description="""
         List of contraction coefficients corresponding to the exponents.
         """,
@@ -255,8 +267,8 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         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)}."
+                    f'Mismatch in number of exponents: expected {self.n_primitive}, '
+                    f'found {len(self.exponents)}.'
                 )
 
         # Resolve combined types
@@ -266,8 +278,8 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> 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)}."
+                        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
@@ -278,16 +290,18 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
                 # Debug: Log split coefficients
                 for t, coeffs in self.coefficient_sets.items():
-                    logger.info(f"{t}-type coefficients: {coeffs}")
+                    logger.info(f'{t}-type coefficients: {coeffs}')
             else:
-                logger.warning(f"No contraction coefficients provided for {self.function_type} type.")
+                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)}."
+                    f'Mismatch in contraction coefficients: expected {self.n_primitive}, '
+                    f'found {len(self.contraction_coefficients)}.'
                 )
 
 
@@ -330,7 +344,7 @@ class AtomCenteredBasisSet(BasisSetComponent):
 
     total_number_of_basis_functions = Quantity(
         type=np.int32,
-        description="",
+        description='',
     )
 
     functional_composition = SubSection(
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 1a06b8ad..4abe1471 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -952,6 +952,42 @@ 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)
+
+
 class GTOIntegralDecomposition(NumericalSettings):
     """
     A general class for integral decomposition techniques for Coulomb and exchange integrals.
diff --git a/tests/test_basis_set.py b/tests/test_basis_set.py
index b209e92a..7e6dd427 100644
--- a/tests/test_basis_set.py
+++ b/tests/test_basis_set.py
@@ -426,8 +426,22 @@ def test_quick_step() -> None:
     [
         ('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
+        (
+            '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(
@@ -536,4 +550,3 @@ def test_atom_centered_basis_set_invalid_data():
     # Call normalize to trigger validation
     with pytest.raises(ValueError, match='Mismatch in number of exponents'):
         invalid_function.normalize(None, None)
-

From 5c15e97c84ed5f16c99c20fc999df0652a805258 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 4 Dec 2024 13:36:25 +0100
Subject: [PATCH 27/37] add method to LocalCorrelation

---
 .../schema_packages/model_method.py           | 55 +++++++++++++++++++
 .../schema_packages/numerical_settings.py     | 24 +++-----
 2 files changed, 64 insertions(+), 15 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index c7b143ff..f66955dc 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1221,3 +1221,58 @@ class DMFT(ModelMethodElectronic):
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
+
+
+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/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 4abe1471..34d59402 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -957,31 +957,25 @@ class OrbitalLocalization(SelfConsistency):
     Numerical settings that control orbital localization.
     """
 
-    localization_method = (
-        Quantity(
-            type=MEnum('FB', 'PM', 'IBO', 'IAOIBO', 'IAOBOYS' 'NEWBOYS' 'AHFB'),
-            description="""
+    localization_method = Quantity(
+        type=MEnum('FB', 'PM', 'IBO', 'IAOIBO', 'IAOBOYS' 'NEWBOYS' 'AHFB'),
+        description="""
         Name of the localization method.
         """,
-        ),
     )
 
-    orbital_window = (
-        Quantity(
-            shape=['*'],
-            description="""
+    orbital_window = Quantity(
+        shape=['*'],
+        description="""
         the Molecular orbital range to be localized.
         """,
-        ),
     )
 
-    core_threshold = (
-        Quantity(
-            type=np.float64,
-            description="""
+    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:

From 2fd6398a939aa1ce9ae03a1f0fa33f4206b0e2b6 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 5 Dec 2024 15:58:58 +0100
Subject: [PATCH 28/37] add total_charge and total_spin to ModelSystem

---
 .../schema_packages/model_system.py               | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

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(

From c9d6118fca0cbab4e325293542a009d4bee706f2 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 10 Dec 2024 14:42:49 +0100
Subject: [PATCH 29/37] add a simple HF class

---
 .../schema_packages/model_method.py           | 30 +++++++++++++++++++
 .../schema_packages/numerical_settings.py     |  2 +-
 2 files changed, 31 insertions(+), 1 deletion(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index f66955dc..1a2a29c5 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1223,6 +1223,36 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
 
+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 PerturbationMethod(ModelMethodElectronic):
     type = Quantity(
         type=MEnum('MP', 'RS', 'BW'),
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 34d59402..477eadd3 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -958,7 +958,7 @@ class OrbitalLocalization(SelfConsistency):
     """
 
     localization_method = Quantity(
-        type=MEnum('FB', 'PM', 'IBO', 'IAOIBO', 'IAOBOYS' 'NEWBOYS' 'AHFB'),
+        type=MEnum('FB', 'PM', 'IBO', 'IAOIBO', 'IAOBOYS', 'NEWBOYS', 'AHFB'),
         description="""
         Name of the localization method.
         """,

From 816c7baf1b8512d473a033d5a6e22545fa8d12da Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 11 Dec 2024 15:20:13 +0100
Subject: [PATCH 30/37] a placeholder for MO and LCAO

---
 .../schema_packages/model_method.py           | 171 ++++++++++++++++++
 .../schema_packages/numerical_settings.py     |  22 ---
 2 files changed, 171 insertions(+), 22 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index 1a2a29c5..c6700e1f 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1223,6 +1223,164 @@ 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.
@@ -1253,6 +1411,19 @@ 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'),
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 477eadd3..48ec6dde 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -980,25 +980,3 @@ class OrbitalLocalization(SelfConsistency):
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
-
-
-class GTOIntegralDecomposition(NumericalSettings):
-    """
-    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,
-        """,
-    )
-
-    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)

From 09b0dfc740357b1d44e95a03d259c5aac816155f Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 2 Jan 2025 09:28:45 +0100
Subject: [PATCH 31/37] add MEnum for GTOIntegralDecomposition

---
 .../schema_packages/model_method.py           | 28 +++++++++----------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index c6700e1f..6fb99f95 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1242,7 +1242,7 @@ class MolecularOrbital(ArchiveSection):
         type=str, description='Symmetry label of the molecular orbital.'
     )
 
-    reference_orbital = SubSection(
+    orbital_ref = SubSection(
         sub_section=OrbitalsState.m_def,
         description='Reference to the underlying atomic orbital state.',
     )
@@ -1275,7 +1275,7 @@ class LCAO(ArchiveSection):
         - 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.
+        TODO: this will be later connected to MCSCF.
         """,
         a_eln=ELNAnnotation(component='EnumEditQuantity'),
     )
@@ -1346,7 +1346,7 @@ def validate_scheme(self, logger: 'BoundLogger') -> bool:
                     f'For {self.reference_type}, the number of alpha and beta electrons must be equal.'
                 )
                 return False
-        if self.reference_type in ['ROHF', 'ROKS']:
+        elif 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.'
@@ -1383,7 +1383,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 class GTOIntegralDecomposition(BaseModelMethod):
     """
-    A general class for integral decomposition techniques for Coulomb and exchange integrals.
+    A general class for integral decomposition techniques for Coulomb and exchange integrals to speed up the calculation.
     Examples:
     Resolution of identity (RI-approximation):
     RI
@@ -1393,23 +1393,22 @@ class GTOIntegralDecomposition(BaseModelMethod):
     """
 
     approximation_type = Quantity(
-        type=str,
+        type=MEnum('RIJ', 'RIJK', 'RIK', 'SENEX', 'RIJCOSX'),
         description="""
-        RIJ, RIK, RIJK,
+        RIJ     : also known as RIJONX, where only Coulomb integrals are approximated.
+        RIJK    : both Coulomb and exchange integrals.
+        RIJCOSX : RIJ for Coulomb and COSX for HF exchange.
         """,
     )
 
     approximated_term = Quantity(
         type=str,
         description="""
-        such as coulomb, exchange, explicit-correlation
-        to be converted to an MEnum later.
+        Such as coulomb, exchange, explicit-correlation, MP2 integrals and so on.
+        TODO: MEnum.
         """,
     )
 
-    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
-
 
 class HartreeFock(ModelMethodElectronic):
     """
@@ -1427,6 +1426,7 @@ class HartreeFock(ModelMethodElectronic):
 class PerturbationMethod(ModelMethodElectronic):
     type = Quantity(
         type=MEnum('MP', 'RS', 'BW'),
+        default='MP',
         description="""
         Perturbation approach. The abbreviations stand for:
         | Abbreviation | Description |
@@ -1435,8 +1435,8 @@ class PerturbationMethod(ModelMethodElectronic):
         | `'RS'`       | Rayleigh-Schrödigner |
         | `'BW'`       | Brillouin-Wigner |
         """,
-        a_eln=ELNAnnotation(component='EnumEditQuantity'),
-    )  # TODO: check if the special symbols are supported
+        # a_eln=ELNAnnotation(component='EnumEditQuantity'),
+    )
 
     order = Quantity(
         type=np.int32,
@@ -1475,5 +1475,5 @@ class LocalCorrelation(ArchiveSection):
 
             # TODO: improve list!
         """,
-        a_eln=ELNAnnotation(component='EnumEditQuantity'),
+        # a_eln=ELNAnnotation(component='EnumEditQuantity'),
     )

From 271f69c8ee7e409e8dc74fb372218f5badbc4e57 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Thu, 2 Jan 2025 10:53:39 +0100
Subject: [PATCH 32/37] modify reference determinants in HF base class

---
 src/nomad_simulations/schema_packages/model_method.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index 6fb99f95..e565d264 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1416,9 +1416,9 @@ class HartreeFock(ModelMethodElectronic):
     """
 
     reference_determinant = Quantity(
-        type=MEnum('UHF', 'RHF', 'ROHF', 'UKS', 'RKS', 'ROKS'),
+        type=MEnum('UHF', 'RHF', 'ROHF'),
         description="""
-        the type of reference determinant.
+        Type of reference determinant.
         """,
     )
 

From d90eac1729d7c6cee2ee4c62366146b80c41e3ff Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Mon, 13 Jan 2025 09:26:51 +0100
Subject: [PATCH 33/37] fix ROKS issue, minor modifications in model_method

---
 .../schema_packages/basis_set.py              |  2 +-
 .../schema_packages/model_method.py           | 77 ++++++++++++++++++-
 2 files changed, 75 insertions(+), 4 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index b19f6618..219eb5ad 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -344,7 +344,7 @@ class AtomCenteredBasisSet(BasisSetComponent):
 
     total_number_of_basis_functions = Quantity(
         type=np.int32,
-        description='',
+        description=""""the total number of basis functions""",
     )
 
     functional_composition = SubSection(
diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index e565d264..d867c52f 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1349,7 +1349,7 @@ def validate_scheme(self, logger: 'BoundLogger') -> bool:
         elif 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.'
+                    f'For {self.reference_type}, there must be unpaired electron(s).'
                 )
                 return False
         return True
@@ -1435,7 +1435,6 @@ class PerturbationMethod(ModelMethodElectronic):
         | `'RS'`       | Rayleigh-Schrödigner |
         | `'BW'`       | Brillouin-Wigner |
         """,
-        # a_eln=ELNAnnotation(component='EnumEditQuantity'),
     )
 
     order = Quantity(
@@ -1475,5 +1474,77 @@ class LocalCorrelation(ArchiveSection):
 
             # TODO: improve list!
         """,
-        # a_eln=ELNAnnotation(component='EnumEditQuantity'),
     )
+
+
+class CoupledCluster(ModelMethodElectronic):
+    """
+    A base section used to define the parameters of a Coupled Cluster calculation.
+    A standard schema is defined, though the most common cases can be summarized in the `type` quantity.
+    """
+
+    type = Quantity(
+        type=str,
+        description="""
+        Coupled Cluster flavor.
+        Examples: CC2, CC3, CCD, CCSD, BCCD, QCCD and so on.
+        The perturbative corrections are not included.
+        """,
+        a_eln=ELNAnnotation(component='StringEditQuantity'),
+    )
+
+    excitation_order = Quantity(
+        type=np.int32,
+        shape=['*'],
+        description="""
+        Orders at which the excitation are used.
+        1 = single, 2 = double, 3 = triple, 4 = quadruple, etc.
+        """
+    )
+
+    reference_determinant = Quantity(
+        type=MEnum('UHF','RHF','ROHF',
+                   'UKS', 'RKS', 'ROKS'),
+        description="""
+        The type of reference determinant.
+        """,
+    )
+
+    perturbation_method = SubSection(sub_section=PerturbationMethod.m_def)
+
+    local_correlation = SubSection(sub_section=LocalCorrelation.m_def)
+
+    perturbative_correction = Quantity(
+        type=MEnum('(T)', '[T]',
+                   '(T0)', '[T0]',
+                   '(Q)'),
+        description="""
+        The type of perturbative corrections.
+        A perturbative correction is different than a perturbation method.
+        """
+    )
+
+    explicit_correlation = Quantity(
+        type=MEnum('F12', 'F12a', 'F12b', 'F12c',
+                   'R12', ''),
+        default='',
+        description="""
+        Explicit correlation treatment.
+        These methods introduce the interelectronic distance coordinate
+        directly into the wavefunction to treat dynamical electron correlation.
+        It can be added linearly (R12) or exponentially (F12).
+        """,
+    )
+
+    is_frozencore = Quantity(
+        type=bool,
+        description="""
+        Flag for frozencore approximation.
+        In post-HF calculation only the valence electrons are typically correlated.
+        The others are kept frozen.
+        FC approximations differ between quantum chemistry codes.
+        """,
+    )
+
+
+

From 95e319c9f008c5241d49017c0006618d3caa5fc6 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 14 Jan 2025 11:24:13 +0100
Subject: [PATCH 34/37] enhance LocalCorrelation and CoupledCluster

---
 .../schema_packages/basis_set.py              |  6 +-
 .../schema_packages/model_method.py           | 74 +++++++-------
 .../schema_packages/numerical_settings.py     | 98 +++++++++++++++++++
 3 files changed, 139 insertions(+), 39 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index 219eb5ad..e9685d8b 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -344,17 +344,13 @@ class AtomCenteredBasisSet(BasisSetComponent):
 
     total_number_of_basis_functions = Quantity(
         type=np.int32,
-        description=""""the total number of basis functions""",
+        description=""""The total number of basis functions.""",
     )
 
     functional_composition = SubSection(
         sub_section=AtomCenteredFunction.m_def, repeats=True
     )
 
-    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
-        # self.name = self.m_def.name
-
 
 class APWBaseOrbital(ArchiveSection):
     """
diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index d867c52f..1c89b94c 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1223,7 +1223,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
 
-class MolecularOrbital(ArchiveSection):
+class BaseMolecularOrbital(ArchiveSection):
     """
     A section representing a single molecular orbital.
     """
@@ -1248,9 +1248,9 @@ class MolecularOrbital(ArchiveSection):
     )
 
 
-class LCAO(ArchiveSection):
+class MolecularOrbitals(ArchiveSection):
     """
-    A base class for molecular orbital schemes used in quantum chemistry calculations.
+    A section for molecular orbital schemes used in quantum chemistry calculations.
     Supports unrestricted (UHF, UKS), restricted (RHF, RKS), and restricted open-shell (ROHF, ROKS) schemes.
     """
 
@@ -1302,7 +1302,7 @@ class LCAO(ArchiveSection):
     )
 
     molecular_orbitals = SubSection(
-        sub_section=MolecularOrbital.m_def,
+        sub_section=BaseMolecularOrbital.m_def,
         repeats=True,
         description="""
         Detailed information about each molecular orbital,
@@ -1398,6 +1398,7 @@ class GTOIntegralDecomposition(BaseModelMethod):
         RIJ     : also known as RIJONX, where only Coulomb integrals are approximated.
         RIJK    : both Coulomb and exchange integrals.
         RIJCOSX : RIJ for Coulomb and COSX for HF exchange.
+        SENEX   : Similar to COSX, relevant for Turbomole.
         """,
     )
 
@@ -1448,8 +1449,8 @@ class PerturbationMethod(ModelMethodElectronic):
     density = Quantity(
         type=MEnum('relaxed', 'unrelaxed'),
         description="""
-        unrelaxed density: MP2 expectation value density
-        relaxed density  : incorporates orbital relaxation
+        unrelaxed density: MP2 expectation value density.
+        relaxed density  : incorporates orbital relaxation.
         """,
     )
 
@@ -1461,18 +1462,17 @@ class LocalCorrelation(ArchiveSection):
     It has a corresponding localization method.
     """
 
-    method = Quantity(
-        type=MEnum('LMP2', 'LCCD', 'LCCSD', 'LCCSD(T)', 'DLPNO-CCSD(T)', 'LocalDFT'),
+    type = Quantity(
+        type=MEnum('PNO', 'domain'),
         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.
+        the type of local correlation.
+        """,
+    )
 
-            # TODO: improve list!
+    ansatz = Quantity(
+        type=MEnum('LMP2', 'LCC', 'LocalDFT'),
+        description="""
+        The underlying ansatz for the local correlation treatment.
         """,
     )
 
@@ -1493,40 +1493,49 @@ class CoupledCluster(ModelMethodElectronic):
         a_eln=ELNAnnotation(component='StringEditQuantity'),
     )
 
+    reference_determinant = Quantity(
+        type=MEnum('UHF', 'RHF', 'ROHF', 'UKS', 'RKS', 'ROKS'),
+        description="""
+        The type of reference determinant.
+        """,
+    )
+
     excitation_order = Quantity(
         type=np.int32,
         shape=['*'],
         description="""
         Orders at which the excitation are used.
         1 = single, 2 = double, 3 = triple, 4 = quadruple, etc.
-        """
+        """,
     )
 
-    reference_determinant = Quantity(
-        type=MEnum('UHF','RHF','ROHF',
-                   'UKS', 'RKS', 'ROKS'),
+    perturbative_correction_order = Quantity(
+        type=np.int32,
+        shape=['*'],
         description="""
-        The type of reference determinant.
+        Orders at which the excitation are used.
+        1 = single, 2 = double, 3 = triple, 4 = quadruple, etc.
         """,
     )
 
-    perturbation_method = SubSection(sub_section=PerturbationMethod.m_def)
-
-    local_correlation = SubSection(sub_section=LocalCorrelation.m_def)
-
     perturbative_correction = Quantity(
-        type=MEnum('(T)', '[T]',
-                   '(T0)', '[T0]',
-                   '(Q)'),
+        type=MEnum('(T)', '[T]', '(T0)', '[T0]', '(Q)'),
         description="""
         The type of perturbative corrections.
         A perturbative correction is different than a perturbation method.
-        """
+        """,
+    )
+
+    perturbation_method = SubSection(sub_section=PerturbationMethod.m_def)
+
+    local_correlation = SubSection(sub_section=LocalCorrelation.m_def, repeats=True)
+
+    gto_integral_decomposition = SubSection(
+        sub_section=GTOIntegralDecomposition.m_def, repeats=True
     )
 
     explicit_correlation = Quantity(
-        type=MEnum('F12', 'F12a', 'F12b', 'F12c',
-                   'R12', ''),
+        type=MEnum('F12', 'F12a', 'F12b', 'F12c', 'R12', ''),
         default='',
         description="""
         Explicit correlation treatment.
@@ -1545,6 +1554,3 @@ class CoupledCluster(ModelMethodElectronic):
         FC approximations differ between quantum chemistry codes.
         """,
     )
-
-
-
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 48ec6dde..f646adb4 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -966,6 +966,7 @@ class OrbitalLocalization(SelfConsistency):
 
     orbital_window = Quantity(
         shape=['*'],
+        type=str,
         description="""
         the Molecular orbital range to be localized.
         """,
@@ -978,5 +979,102 @@ class OrbitalLocalization(SelfConsistency):
         """,
     )
 
+
+class PNOSettings(NumericalSettings):
+    """Numerical settings that control pair natural orbitals (PNOs).
+    The nomenclature has been adapted from Molpro.
+    """
+
+    domain_connectivity = Quantity(
+        type=int,
+        description="""
+        the connectivity parameter for domain extension.
+        """,
+    )
+
+    domain_radius = Quantity(
+        type=int,
+        description="""
+        the radius parameter for domain extension.
+        """,
+    )
+
+    t_domain_osv_occ = Quantity(
+        type=np.float32,
+        description="""
+        OSV domain occupation number threshold.
+        """,
+    )
+
+    t_occ_lmp2 = Quantity(
+        type=np.float32,
+        description="""
+        LMP2 PNO domains (occ. number threshold).
+        """,
+    )
+
+    t_energy_lmp2 = Quantity(
+        type=np.float32,
+        description="""
+        LMP2 PNO domains (energy threshold).
+        """,
+    )
+
+    t_occ_lccsd = Quantity(
+        type=np.float32,
+        description="""
+        LCCSD PNO domains (occ. number threshold).
+        """,
+    )
+
+    t_energy_lccsd = Quantity(
+        type=np.float32,
+        description="""
+        LCCSD PNO domains (energy threshold).
+        """,
+    )
+
+    t_close_pair = Quantity(
+        type=str,
+        description="""
+        close pair energy threshold.
+        """,
+    )
+
+    t_weak_pair = Quantity(
+        type=np.float32,
+        description="""
+        weak pair energy threshold.
+        """,
+    )
+
+    t_distant_pair = Quantity(
+        type=np.float32,
+        description="""
+        distant pair energy threshold
+        """,
+    )
+
+    t_verydistant_pair = Quantity(
+        type=np.float32,
+        description="""
+        very distant pair energy threshold
+        """,
+    )
+
+    t_triples_preselection = Quantity(
+        type=np.float32,
+        description="""
+        preselection of triples list.
+        """,
+    )
+
+    t_triples_iteration = Quantity(
+        type=np.float32,
+        description="""
+        selection of triples for iterations
+        """,
+    )
+
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)

From 6f3e30b13f9ba8cbd25273ed427a6541dbf6add4 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Wed, 15 Jan 2025 13:54:51 +0100
Subject: [PATCH 35/37] add a placeholder MolecularOrbitalsState

---
 .../schema_packages/atoms_state.py            | 69 ++++++++++++++++
 .../schema_packages/model_method.py           | 82 ++++---------------
 2 files changed, 86 insertions(+), 65 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/atoms_state.py b/src/nomad_simulations/schema_packages/atoms_state.py
index 32fbdd11..8e2c7327 100644
--- a/src/nomad_simulations/schema_packages/atoms_state.py
+++ b/src/nomad_simulations/schema_packages/atoms_state.py
@@ -641,3 +641,72 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
             self.chemical_symbol = self.resolve_chemical_symbol(logger=logger)
         if self.atomic_number is None:
             self.atomic_number = self.resolve_atomic_number(logger=logger)
+
+
+class MolecularOrbitalsState(Entity):
+    """
+    A base section to define molecular orbitals.
+    """
+
+    symmetry_label = Quantity(
+        type=str,
+        description="""
+        Symmetry label of the molecular orbital (e.g., 'sigma', 'pi', 'delta').
+        """,
+    )
+
+    energy = Quantity(
+        type=np.float64,
+        #unit='eV',
+        description="""
+        Energy of the molecular orbital.
+        """,
+    )
+
+    occupation = Quantity(
+        type=np.float64,
+        description="""
+        Occupation of the molecular orbital. This value is typically an integer (0 or 2)
+        in closed-shell systems, but can be fractional in open-shell or spin-polarized
+        calculations.
+        """,
+    )
+
+    spin = Quantity(
+        type=MEnum('alpha', 'beta'),
+        description="""
+        Spin of the molecular orbital. 'alpha' corresponds to spin-up, 'beta' corresponds
+        to spin-down.
+        """,
+    )
+
+    coefficients = Quantity(
+        type=np.float64,
+        shape=['number_of_atoms', 'number_of_basis_functions'],
+        description="""
+        Coefficients of the molecular orbital expressed as a linear combination of atomic orbitals.
+        The shape corresponds to the number of atoms and their associated basis functions.
+        """,
+    )
+
+    atom_contributions = SubSection(
+        sub_section=AtomsState.m_def,
+        repeats=True,
+        description="""
+        Contribution of each atom to the molecular orbital, as defined by its basis functions.
+        """,
+    )
+
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        super().normalize(archive, logger)
+
+        # Validation: Ensure occupation values are consistent
+        if self.occupation is not None and (self.occupation < 0 or self.occupation > 2):
+            logger.error("The molecular orbital occupation must be between 0 and 2.")
+
+        # Validation: Ensure coefficients are provided if atom contributions are defined
+        if self.atom_contributions and self.coefficients is None:
+            logger.error(
+                "Coefficients must be defined when atom contributions are provided."
+            )
+
diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index 1c89b94c..29266538 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1223,31 +1223,6 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
 
-class BaseMolecularOrbital(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.'
-    )
-
-    orbital_ref = SubSection(
-        sub_section=OrbitalsState.m_def,
-        description='Reference to the underlying atomic orbital state.',
-    )
-
-
 class MolecularOrbitals(ArchiveSection):
     """
     A section for molecular orbital schemes used in quantum chemistry calculations.
@@ -1268,18 +1243,6 @@ class MolecularOrbitals(ArchiveSection):
         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 MCSCF.
-        """,
-        a_eln=ELNAnnotation(component='EnumEditQuantity'),
-    )
-
     n_alpha_electrons = Quantity(
         type=np.int32,
         description="""
@@ -1301,16 +1264,16 @@ class MolecularOrbitals(ArchiveSection):
         """,
     )
 
-    molecular_orbitals = SubSection(
-        sub_section=BaseMolecularOrbital.m_def,
-        repeats=True,
-        description="""
-        Detailed information about each molecular orbital,
-        including energy, occupation, and symmetry label.
-        """,
-    )
+    # molecular_orbitals = SubSection(
+    #     sub_section=BaseMolecularOrbital.m_def,
+    #     repeats=True,
+    #     description="""
+    #     Detailed information about each molecular orbital,
+    #     including energy, occupation, and symmetry label.
+    #     """,
+    # )
 
-    total_spin = Quantity(
+    sz_projection = Quantity(
         type=np.float64,
         description="""
         Total spin of the system defined as S = (n_alpha_electrons - n_beta_electrons) / 2.
@@ -1364,24 +1327,8 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         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):
+class IntegralDecomposition(BaseModelMethod):
     """
     A general class for integral decomposition techniques for Coulomb and exchange integrals to speed up the calculation.
     Examples:
@@ -1463,16 +1410,20 @@ class LocalCorrelation(ArchiveSection):
     """
 
     type = Quantity(
-        type=MEnum('PNO', 'domain'),
+        type=MEnum('PNO', 'LPNO', 'DLPNO'),
         description="""
         the type of local correlation.
         """,
     )
 
     ansatz = Quantity(
-        type=MEnum('LMP2', 'LCC', 'LocalDFT'),
+        type=MEnum('MP2', 'CC', 'DH-DFT'),
         description="""
         The underlying ansatz for the local correlation treatment.
+        LMP2
+        LCC
+        DH-DFT: double-hybrid DFT
+        ...
         """,
     )
 
@@ -1493,6 +1444,7 @@ class CoupledCluster(ModelMethodElectronic):
         a_eln=ELNAnnotation(component='StringEditQuantity'),
     )
 
+    # add a real reference determinant reference
     reference_determinant = Quantity(
         type=MEnum('UHF', 'RHF', 'ROHF', 'UKS', 'RKS', 'ROKS'),
         description="""

From 92fd95e91a612f05d5742a85b69fe98ac841d04b Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 28 Jan 2025 10:28:48 +0100
Subject: [PATCH 36/37] improved class descriptions, new MolecularOrbitals
 class, improved HartreeFock section

---
 .../schema_packages/atoms_state.py            | 103 ++++---
 .../schema_packages/basis_set.py              | 116 +++++++-
 .../schema_packages/model_method.py           | 275 +++++++++---------
 .../schema_packages/numerical_settings.py     |  60 +++-
 4 files changed, 359 insertions(+), 195 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/atoms_state.py b/src/nomad_simulations/schema_packages/atoms_state.py
index 8e2c7327..b2b97b06 100644
--- a/src/nomad_simulations/schema_packages/atoms_state.py
+++ b/src/nomad_simulations/schema_packages/atoms_state.py
@@ -643,70 +643,105 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
             self.atomic_number = self.resolve_atomic_number(logger=logger)
 
 
-class MolecularOrbitalsState(Entity):
+class MolecularOrbitals(Entity):
     """
-    A base section to define molecular orbitals.
+    This class stores all molecular orbitals (MO) in a single container, with each Quantity using
+    arrays indexed by mo_num and ao_num.
+
+        Comparison to TREXIO:
+      - mo/type           -> mo_type
+      - mo/num            -> mo_num
+      - mo/coefficient    -> coefficient
+      - mo/coefficient_im -> coefficient_im
+      - mo/symmetry       -> symmetry
+      - mo/occupation     -> occupation
+      - mo/energy         -> energy
+      - mo/spin           -> spin
+
     """
 
-    symmetry_label = Quantity(
+    mo_type = Quantity(
         type=str,
+        shape=['mo_num'],
         description="""
-        Symmetry label of the molecular orbital (e.g., 'sigma', 'pi', 'delta').
+        Type of the molecular orbitals
+         e.g. 'canonical', 'localized'.
+        In case of CASSCF calculations, there will be orbital subspaces of different nature.
+        E.g. : 
+        Internal orbitals : canonical
+        Active orbitals   : natural
+        Virtual orbitals  : canonical
         """,
     )
 
-    energy = Quantity(
+    mo_num = Quantity(
+        type=np.int32,
+        description="""
+        Number of molecular orbitals.
+        """,
+    )
+
+    ao_num = Quantity(
+        type=np.int32,
+        description="""
+        Number of atomic orbitals or basis functions (often needed for coefficient shape).
+        Corresponds to the 'ao.num' dimension in TREXIO.
+        """,
+    )
+
+    coefficient = Quantity(
         type=np.float64,
-        #unit='eV',
+        shape=['mo_num', 'ao_num'],
         description="""
-        Energy of the molecular orbital.
+        Real part of the MO coefficients. The shape is
+        [mo.num, ao.num], meaning each row corresponds to one MO, and each column
+        to one atomic orbital (or basis function).
         """,
     )
 
-    occupation = Quantity(
+    coefficient_im = Quantity(
         type=np.float64,
+        shape=['mo_num', 'ao_num'],
         description="""
-        Occupation of the molecular orbital. This value is typically an integer (0 or 2)
-        in closed-shell systems, but can be fractional in open-shell or spin-polarized
-        calculations.
+        Imaginary part of the MO coefficients. The shape is
+        [mo.num, ao.num]. This array may be omitted or set to zero if the orbitals
+        are purely real.
         """,
     )
 
-    spin = Quantity(
-        type=MEnum('alpha', 'beta'),
+    symmetry = Quantity(
+        type=str,
+        shape=['mo_num'],
         description="""
-        Spin of the molecular orbital. 'alpha' corresponds to spin-up, 'beta' corresponds
-        to spin-down.
+        Symmetry label for each MO, e.g. group-theory labels or
+        simpler 'sigma', 'pi', 'delta'.
         """,
     )
 
-    coefficients = Quantity(
+    occupation = Quantity(
         type=np.float64,
-        shape=['number_of_atoms', 'number_of_basis_functions'],
+        shape=['mo_num'],
         description="""
-        Coefficients of the molecular orbital expressed as a linear combination of atomic orbitals.
-        The shape corresponds to the number of atoms and their associated basis functions.
+        Occupation numbers for each MO. Typically in [0, 2]
+        for closed-shell systems, but might be fractional in open-shell systems or multi-reference calculations.
         """,
     )
 
-    atom_contributions = SubSection(
-        sub_section=AtomsState.m_def,
-        repeats=True,
+    energy = Quantity(
+        type=np.float64,
+        shape=['mo_num'],
         description="""
-        Contribution of each atom to the molecular orbital, as defined by its basis functions.
+        Orbital energies for each MO.
         """,
     )
 
-    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)
-
-        # Validation: Ensure occupation values are consistent
-        if self.occupation is not None and (self.occupation < 0 or self.occupation > 2):
-            logger.error("The molecular orbital occupation must be between 0 and 2.")
+    spin = Quantity(
+        type=np.int32,
+        shape=['mo_num'],
+        description="""
+        Spin channel for each MO if this is an unrestricted open-shell set. 
+        Typically 0 for alpha, 1 for beta.
+        """,
+    )
 
-        # Validation: Ensure coefficients are provided if atom contributions are defined
-        if self.atom_contributions and self.coefficients is None:
-            logger.error(
-                "Coefficients must be defined when atom contributions are provided."
-            )
 
diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index e9685d8b..b66d2f50 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -187,7 +187,24 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 class AtomCenteredFunction(ArchiveSection):
     """
-    Specifies a single function (term) in an atom-centered basis set.
+    Specifies a single contracted basis function in an atom-centered basis set. 
+    
+    In many quantum-chemistry codes, an atom-centered basis set is composed of 
+    several "shells," each shell containing one or more basis functions of a certain 
+    angular momentum. For instance, a shell of p-type orbitals (L=1) typically 
+    consists of 3 degenerate functions (p_x, p_y, p_z) if `harmonic_type='cartesian'`
+    or 3 spherical harmonics if `harmonic_type='spherical'`.
+    
+    A single "atom-centered function" can be a linear combination of multiple 
+    primitive Gaussians (or Slater-type orbitals, STOs). 
+    In practice, these contract together to form the final basis function used by 
+    the SCF or post-SCF method. Often, each contraction is labeled by its 
+    angular momentum (e.g., s, p, d, f) and a set of exponents and coefficients.
+
+    **References**:
+      - T. Helgaker, P. Jørgensen, J. Olsen, *Molecular Electronic-Structure Theory*, Wiley (2000).
+      - F. Jensen, *Introduction to Computational Chemistry*, 2nd ed., Wiley (2007).
+      - J. B. Foresman, Æ. Frisch, *Exploring Chemistry with Electronic Structure Methods*, Gaussian Inc.
     """
 
     harmonic_type = Quantity(
@@ -197,7 +214,13 @@ class AtomCenteredFunction(ArchiveSection):
         ),
         default='spherical',
         description="""
-        Specifies whether the basis functions are spherical-harmonic or cartesian functions.
+        Specifies whether the basis functions are expanded in **spherical** (pure) 
+        harmonics or **cartesian** harmonics. Many modern quantum-chemistry codes 
+        default to *spherical harmonics* for d, f, g..., which eliminates the 
+        redundant functions found in the cartesian sets.
+        
+        - `'spherical'` : (2l+1) functions for a shell of angular momentum l
+        - `'cartesian'` : (l+1)(l+2)/2 functions for that shell (extra functions appear)
         """,
     )
 
@@ -218,16 +241,28 @@ class AtomCenteredFunction(ArchiveSection):
             'spdf',
         ),
         description="""
-        L=a+b+c
-        The angular momentum of GTO to be added.
+        Symbolic label for the **angular momentum** of this contracted function. 
+        Typical single-letter labels:
+          - 's' => L=0
+          - 'p' => L=1
+          - 'd' => L=2
+          - 'f' => L=3
+          - 'g' => L=4
+          - 'h', 'i', etc. => still higher angular momenta
+        Combined labels like 'sp' or 'spdf' indicate a **combined shell** in which 
+        multiple angular momenta share exponents. For example, in some older Pople 
+        basis sets, an 'sp' shell has an s- and p-type function sharing the same 
+        exponents but different contraction coefficients.
         """,
     )
 
     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.
+        Number of **primitive** functions in this contracted basis function. 
+        For example, in a contracted Gaussian-type orbital (GTO) approach, each basis 
+        function might be built from a sum of `n_primitive` Gaussians with different 
+        exponents, each scaled by a contraction coefficient.
         """,
     )
 
@@ -235,7 +270,10 @@ class AtomCenteredFunction(ArchiveSection):
         type=np.float32,
         shape=['n_primitive'],
         description="""
-        List of exponents for the basis function.
+        The **exponents** of each primitive basis function. 
+        In a Gaussian basis set (GTO), these are the alpha_i in 
+        exp(-alpha_i * r^2). In a Slater-type basis (STO), they'd be 
+        exp(-zeta_i * r). Typically sorted from largest to smallest.
         """,
     )
 
@@ -243,14 +281,21 @@ class AtomCenteredFunction(ArchiveSection):
         type=np.float32,
         shape=['*'],  # Flexible shape to handle combined types (e.g. SP, SPD..)
         description="""
-        List of contraction coefficients corresponding to the exponents.
+        The **contraction coefficients** associated with each primitive exponent. 
+        In the simplest case (pure s- or p-function), this array has length 
+        equal to `n_primitive`. For combined shells (like 'sp'), the length 
+        might be `2 * n_primitive`, because you have separate coefficients 
+        for the s-part and the p-part.
         """,
     )
 
     point_charge = Quantity(
         type=np.float32,
         description="""
-        the value of the point charge.
+        If using a basis function that explicitly includes a point-charge or an 
+        ECP-like pseudo-component, this field can store that net charge. 
+        Typically 0 for standard GTO or STO expansions. 
+        Some extended basis concepts (or embedded charges) might set it differently.
         """,
     )
 
@@ -271,7 +316,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
                     f'found {len(self.exponents)}.'
                 )
 
-        # Resolve combined types
+        # For combined shells (like 'sp', 'spd', etc.), ensure the coefficient array is large enough
         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:
@@ -307,13 +352,36 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 class AtomCenteredBasisSet(BasisSetComponent):
     """
-    Defines an atom-centered basis set.
+    Defines an **atom-centered basis set** for quantum chemistry calculations.
+    Unlike plane-wave methods, these expansions are typically built around each atom's 
+    position, using either:
+      - Slater-type orbitals (STO)
+      - Gaussian-type orbitals (GTO)
+      - Numerical atomic orbitals (NAO)
+      - Effective-core potentials or point-charges (PC, cECP, etc.)
+    
+    This section references multiple `AtomCenteredFunction` objects, each describing
+    a single contracted function or shell. Additionally, one can label the overall 
+    basis set name (e.g., "cc-pVTZ", "def2-SVP", "6-31G**") and specify the high-level 
+    role of the basis set in the calculation.
+
+    **Common examples**:
+      - **Pople basis** (3-21G, 6-31G(d), 6-311++G(2df,2pd), etc.)
+      - **Dunning correlation-consistent** (cc-pVDZ, cc-pVTZ, aug-cc-pVTZ, etc.)
+      - **Slater basis** used in ADF, for instance
+      - **ECP** (Effective Core Potential) expansions like LANL2DZ or SDD for transition metals
+
+    **References**:
+      - F. Jensen, *Introduction to Computational Chemistry*, 2nd ed., Wiley (2007).
+      - A. Szabo, N. S. Ostlund, *Modern Quantum Chemistry*, McGraw-Hill (1989).
+      - T. H. Dunning Jr., J. Chem. Phys. 90, 1007 (1989) for correlation-consistent basis sets.
     """
 
     basis_set = Quantity(
         type=str,
         description="""
-        name of the basis set.
+        **Name** or label of the basis set as recognized by the code or standard 
+        library. Examples: "6-31G*", "cc-pVTZ", "def2-SVP", "STO-3G", "LANL2DZ" (ECP).
         """,
     )
 
@@ -326,7 +394,15 @@ class AtomCenteredBasisSet(BasisSetComponent):
             'PC',  # Point charges
         ),
         description="""
-        Type of the basis set, e.g. STO or GTO.
+        The **functional form** of the basis set:
+          - 'STO': Slater-type orbitals
+          - 'GTO': Gaussian-type orbitals
+          - 'NAO': Numerical atomic orbitals
+          - 'cECP': Some variant of a "capped" or shape-consistent ECP 
+          - 'PC': Point charges (or ghost basis centers)
+
+        If a code uses a mixture (e.g., GTO + ECP), it might either store them 
+        as separate `AtomCenteredBasisSet` sections or unify them if the code does so internally.
         """,
     )
 
@@ -335,16 +411,24 @@ class AtomCenteredBasisSet(BasisSetComponent):
             'orbital',
             'auxiliary_scf',
             'auxiliary_post_hf',
-            'cabs',  # complementary auxiliary basis set
+            'cabs',  
         ),
         description="""
-        The role of the basis set.
+        The role of this basis set in the calculation:
+          - 'orbital': main orbital basis for the SCF
+          - 'auxiliary_scf': used for RI-J or density fitting in SCF
+          - 'auxiliary_post_hf': used in MP2, CC, etc. 
+          - 'cabs': complementary auxiliary basis for explicitly correlated (F12) methods.
         """,
     )
 
     total_number_of_basis_functions = Quantity(
         type=np.int32,
-        description=""""The total number of basis functions.""",
+        description="""
+        The **total** number of contracted basis functions in this entire set. 
+        This is typically the sum of all `(2l+1)` or cartesian expansions across
+        all shells on all relevant atoms (within the scope of this section).
+        """,
     )
 
     functional_composition = SubSection(
diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index 29266538..d3055ffc 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -11,7 +11,11 @@
     from nomad.metainfo import Context
     from structlog.stdlib import BoundLogger
 
-from nomad_simulations.schema_packages.atoms_state import CoreHole, OrbitalsState
+from nomad_simulations.schema_packages.atoms_state import (
+    CoreHole,
+    MolecularOrbitals,
+    OrbitalsState,
+)
 from nomad_simulations.schema_packages.model_system import ModelSystem
 from nomad_simulations.schema_packages.numerical_settings import NumericalSettings
 from nomad_simulations.schema_packages.utils import is_not_representative
@@ -1223,120 +1227,24 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
 
-class MolecularOrbitals(ArchiveSection):
-    """
-    A section 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'),
-    )
-
-    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=BaseMolecularOrbital.m_def,
-    #     repeats=True,
-    #     description="""
-    #     Detailed information about each molecular orbital,
-    #     including energy, occupation, and symmetry label.
-    #     """,
-    # )
-
-    sz_projection = 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
-        elif 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 unpaired electron(s).'
-                )
-                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.')
-
-
 class IntegralDecomposition(BaseModelMethod):
     """
-    A general class for integral decomposition techniques for Coulomb and exchange integrals to speed up the calculation.
-    Examples:
-    Resolution of identity (RI-approximation):
-    RI
-    RIJK
-    ....
-    Chain-of-spheres (COSX) algorithm for exchange: doi:10.1016/j.chemphys.2008.10.036
+    A general class for integral decomposition techniques that approximate
+    Coulomb and/or exchange integrals to reduce computational cost in quantum
+    chemistry. Examples include:
+
+      - Resolution of the Identity (RI, a.k.a. density fitting)
+      - Chain-of-Spheres exchange (COSX)
+      - Cholesky Decomposition (CD)
+      - Other domain-based or rank-reduced approximations
+
+    Typical references:
+      - F. Weigend, M. Häser, The RI-MP2 method: Algorithmic
+        implementation of efficient, approximate MP2 theories,
+        Theor. Chem. Acc. 97, 331-340 (1997).
+      - S. Hättig, F. Weigend, J. Chem. Phys. 113, 5154 (2000). (RI-J)
+      - Neese et al., “Chain-of-spheres algorithms for HF exchange,”
+        Chem. Phys. 356 (2008), 98-109.
     """
 
     approximation_type = Quantity(
@@ -1350,17 +1258,34 @@ class IntegralDecomposition(BaseModelMethod):
     )
 
     approximated_term = Quantity(
-        type=str,
+        type=MEnum('coulomb', 'exchange', 'mp2', 'cc', 'explicit_correlation', 'other'),
         description="""
-        Such as coulomb, exchange, explicit-correlation, MP2 integrals and so on.
-        TODO: MEnum.
+        Which terms are approximated by this method:
+          - 'coulomb': only the J integrals
+          - 'exchange': only K integrals
+          - 'mp2': MP2 integrals
+          - 'cc': Coupled Cluster integrals
+          - 'explicit_correlation': e.g. R12, F12 integrals
         """,
     )
 
 
 class HartreeFock(ModelMethodElectronic):
     """
-    A base section for Hartree Fock ansatz.
+    Defines a Hartree–Fock (HF) calculation using a specified reference determinant 
+    (RHF, UHF, or ROHF).
+
+    In HF theory:
+      - RHF  = Restricted Hartree–Fock, for closed-shell systems.
+      - UHF  = Unrestricted Hartree–Fock, allows different orbitals for alpha/beta spin.
+      - ROHF = Restricted Open-Shell Hartree–Fock, a partially restricted approach for 
+               open-shell systems.
+
+    **References**:
+      - Roothaan, C. C. J. (1951). "New Developments in Molecular Orbital Theory." 
+        Rev. Mod. Phys. 23, 69.
+      - Szabo, A., & Ostlund, N. S. (1989). *Modern Quantum Chemistry*. McGraw-Hill.
+      - Jensen, F. (2007). *Introduction to Computational Chemistry*. 2nd ed., Wiley.
     """
 
     reference_determinant = Quantity(
@@ -1370,6 +1295,42 @@ class HartreeFock(ModelMethodElectronic):
         """,
     )
 
+    molecular_orbitals_ref = Quantity(
+        type=MolecularOrbitals,
+        description="""
+        Final, converged molecular orbitals from the Hartree–Fock SCF procedure. This includes
+        orbital energies, occupations, symmetry labels, and spin channels if applicable.
+        """,
+    )
+
+    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
+        """
+        Perform minimal consistency checks between the HF reference determinant
+        and the final molecular orbitals spin array (if available).
+        """
+        super().normalize(archive, logger)
+
+        if self.molecular_orbitals is not None:
+            # If the user is storing a spin array for MOs, check for consistency with the determinant
+            mo_spin = self.molecular_orbitals.spin
+            if mo_spin is not None and len(mo_spin) > 0:
+                unique_spins = np.unique(mo_spin)
+                if self.reference_determinant == 'RHF':
+                    # For RHF, we typically expect spin=0 (alpha) only
+                    if not np.all(mo_spin == 0):
+                        logger.warning(
+                            "RHF reference used, but molecular_orbitals.spin contains non-zero spin indices."
+                        )
+                elif self.reference_determinant == 'UHF':
+                    # UHF often has spin=0 for alpha and spin=1 for beta
+                    # If we only see spin=0, that's effectively no spin polarization
+                    if len(unique_spins) == 1 and unique_spins[0] == 0:
+                        logger.info(
+                            "UHF reference chosen, but only alpha spin found in MOs (spin=0). "
+                            "This might still be valid if spin polarization is zero."
+                        )
+                # For ROHF, spin indexing can vary across codes, so no strict check here.
+
 
 class PerturbationMethod(ModelMethodElectronic):
     type = Quantity(
@@ -1379,7 +1340,7 @@ class PerturbationMethod(ModelMethodElectronic):
         Perturbation approach. The abbreviations stand for:
         | Abbreviation | Description |
         | ------------ | ----------- |
-        | `'MP'`       | Moller-Plesset |
+        | `'MP'`       | Møller-Plesset |
         | `'RS'`       | Rayleigh-Schrödigner |
         | `'BW'`       | Brillouin-Wigner |
         """,
@@ -1401,18 +1362,43 @@ class PerturbationMethod(ModelMethodElectronic):
         """,
     )
 
+    spin_component_scaling = Quantity(
+        type=MEnum('none', 'SCS', 'SOS', 'custom'),
+        default='none',
+        description="""
+        Spin-component scaling approach for perturbation methods:
+          - none  : no spin-component scaling
+          - SCS   : spin-component scaled (Grimme's approach, https://doi.org/10.1002/wcms.1110)
+          - SOS   : spin-opposite scaled
+          - custom: user-defined scaling factors
+        This is typically relevant only for MP2 calculations.
+        """,
+    )
+
+
 
 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.
+    A base section used to define the parameters of a local correlation method for
+    post-HF calculations, e.g. LMP2, LCC, or domain-based local pair natural orbitals
+    (PNO, LPNO, DLPNO) in coupled cluster or double-hybrid DFT.
+
+    Typical references:
+      - Pulay, Chem. Phys. Lett. 100, 151 (1983) (LMP2 concept).
+      - G. Knizia, G. K.-L. Chan, “Density Matrix Embedding,” J. Chem. Theory Comput. 9, 1428 (2013).
+      - F. Neese, “The ORCA program system,” WIREs Comput. Mol. Sci. 2, 73-78 (2012)
+        for DLPNO approaches.
     """
 
     type = Quantity(
-        type=MEnum('PNO', 'LPNO', 'DLPNO'),
+        type=MEnum('PNO', 'LPNO', 'DLPNO', 'LMP2', 'other'),
         description="""
-        the type of local correlation.
+        The local correlation flavor:
+          - 'PNO'   : Pair Natural Orbitals in generic form
+          - 'LPNO'  : Local PNO approach
+          - 'DLPNO' : Domain-based Local PNO approach
+          - 'LMP2'  : Local MP2 (Pulay approach)
+          - 'other' : Another local correlation scheme
         """,
     )
 
@@ -1437,18 +1423,24 @@ class CoupledCluster(ModelMethodElectronic):
     type = Quantity(
         type=str,
         description="""
-        Coupled Cluster flavor.
-        Examples: CC2, CC3, CCD, CCSD, BCCD, QCCD and so on.
-        The perturbative corrections are not included.
+        String labeling the Coupled Cluster flavor (e.g., CC2, CC3, CCD, CCSD, CCSDT, etc.).
+        If a known standard approach, it might match these examples:
+          - CC2, CC3  : approximate methods for excited states
+          - CCD       : Coupled Cluster Doubles
+          - CCSD      : Singles and Doubles
+          - CCSDT     : Singles, Doubles, and Triples
+          - CCSDTQ    : Singles, Doubles, Triples, and Quadruples
+        By default, the "perturbative corrections" like (T) are not included in this string.
         """,
         a_eln=ELNAnnotation(component='StringEditQuantity'),
     )
 
-    # add a real reference determinant reference
     reference_determinant = Quantity(
         type=MEnum('UHF', 'RHF', 'ROHF', 'UKS', 'RKS', 'ROKS'),
         description="""
-        The type of reference determinant.
+        The type of reference determinant on which the Coupled Cluster expansion is built.
+          - UHF / RHF / ROHF : common for wavefunction-based reference
+          - UKS / RKS / ROKS : if a KS-DFT reference is used (rare but possible).        
         """,
     )
 
@@ -1456,8 +1448,13 @@ class CoupledCluster(ModelMethodElectronic):
         type=np.int32,
         shape=['*'],
         description="""
-        Orders at which the excitation are used.
-        1 = single, 2 = double, 3 = triple, 4 = quadruple, etc.
+        The excitation orders explicitly included in the cluster operator, e.g. [1,2]
+        for CCSD. 
+        - 1 = singles 
+        - 2 = doubles 
+        - 3 = triples 
+        - 4 = quadruples, etc.
+        Example: CCSDT => [1, 2, 3].
         """,
     )
 
@@ -1465,16 +1462,22 @@ class CoupledCluster(ModelMethodElectronic):
         type=np.int32,
         shape=['*'],
         description="""
-        Orders at which the excitation are used.
-        1 = single, 2 = double, 3 = triple, 4 = quadruple, etc.
+        The excitation orders included only in a perturbative manner.
+        For instance, in CCSD(T), singles and doubles are solved iteratively,
+        while triples appear as a perturbative correction => [3].
         """,
     )
 
     perturbative_correction = Quantity(
         type=MEnum('(T)', '[T]', '(T0)', '[T0]', '(Q)'),
         description="""
-        The type of perturbative corrections.
-        A perturbative correction is different than a perturbation method.
+        Label for the perturbative corrections:
+          - '(T)'   : standard perturbative triples
+          - '[T]'   : Brueckner-based or other variant
+          - '(T0)'  : approximate version of (T)
+          - '[T0]'  : approximate, typically for Brueckner references
+          - '(Q)'   : perturbative quadruples, e.g., CCSDT(Q)
+          - 'none'  : no perturbative correction
         """,
     )
 
@@ -1482,8 +1485,8 @@ class CoupledCluster(ModelMethodElectronic):
 
     local_correlation = SubSection(sub_section=LocalCorrelation.m_def, repeats=True)
 
-    gto_integral_decomposition = SubSection(
-        sub_section=GTOIntegralDecomposition.m_def, repeats=True
+    integral_decomposition = SubSection(
+        sub_section=IntegralDecomposition.m_def, repeats=True
     )
 
     explicit_correlation = Quantity(
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index f646adb4..eca33e56 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -175,9 +175,9 @@ class NumericalIntegration(NumericalSettings):
     integration_thresh = Quantity(
         type=np.float64,
         description="""
-        Accuracy threshold for integration grid.
-        GRIDTHR in Molpro.
-        BFCut in ORCA.
+        An accuracy threshold for the integration grid, controlling how fine the 
+        discretization is. Some programs label it "integral accuracy" or "grid accuracy".
+        For instance, GRIDTHR in Molpro or BFCut in ORCA.
         """,
     )
 
@@ -954,13 +954,37 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 class OrbitalLocalization(SelfConsistency):
     """
-    Numerical settings that control orbital localization.
+    Numerical settings that control orbital localization, typically applied
+    to transform canonical molecular orbitals into localized orbitals. 
+    These localized orbitals can then be used for:
+      - Local correlation methods (e.g., LMP2, local CC)
+      - Interpretable chemical analysis (e.g., identifying bonds, lone pairs)
+      - Faster post-HF or excited-state calculations
+
+    Inherit from `SelfConsistency` because some localization algorithms are
+    iterative, requiring thresholds and iteration limits akin to an SCF loop.
+
+    References:
+      - R. F. Boys, "Construction of Some Molecular Orbitals to Be Approximately Invariant 
+        for Changes in Molecular Conformation," Rev. Mod. Phys. 32, 296 (1960). [Boys]
+      - J. Pipek, P. G. Mezey, J. Chem. Phys. 90, 4916 (1989). [PM method]
+      - J. L. Knizia, "Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory 
+        and Chemical Concepts," J. Chem. Theory Comput. 9, 4834 (2013). [IBO method]
+      - P. Pulay, Chem. Phys. Lett. 100, 151 (1983). [Local correlation context]
     """
 
     localization_method = Quantity(
         type=MEnum('FB', 'PM', 'IBO', 'IAOIBO', 'IAOBOYS', 'NEWBOYS', 'AHFB'),
         description="""
-        Name of the localization method.
+        The chosen localization algorithm:
+          - FB     : Foster-Boys method (a variant of Boys localization)
+          - PM     : Pipek–Mezey localization
+          - IBO    : Intrinsic Bond Orbitals (Knizia)
+          - IAOIBO : Combination of Intrinsic Atomic Orbitals + Intrinsic Bond Orbitals
+          - IAOBOYS: IAO-based Boys approach (a custom variant)
+          - NEWBOYS: Another specialized Boys-based approach
+          - AHFB   : Augmented Hessian Foster-Boys, an alternative with different
+                     optimization steps
         """,
     )
 
@@ -968,21 +992,39 @@ class OrbitalLocalization(SelfConsistency):
         shape=['*'],
         type=str,
         description="""
-        the Molecular orbital range to be localized.
+        The set of molecular orbitals to be localized (e.g., 'occupied only', 'valence only',
+        or a range like `[5..20]`). This can be code- or user-defined. For instance,
+        'occupied' might indicate all occupied orbitals except possibly core orbitals,
+        'valence' might skip deep core and also skip high-lying virtual orbitals, etc.
         """,
     )
 
     core_threshold = Quantity(
         type=np.float64,
         description="""
-        the energy window for the first occupied MO to be localized (in a.u.).
+        The energy window or threshold defining which orbitals are considered core,
+        and thus excluded from the localization procedure. For example, an orbital with
+        an energy below -20.0 eV might be automatically treated as a frozen core orbital and
+        not localized.
         """,
     )
 
 
 class PNOSettings(NumericalSettings):
-    """Numerical settings that control pair natural orbitals (PNOs).
-    The nomenclature has been adapted from Molpro.
+    """
+    Numerical settings that control **Pair Natural Orbitals (PNOs)** in local correlation approaches.
+    PNOs are a compact representation of the virtual orbital space for each pair (or domain)
+    of occupied orbitals, improving efficiency in post-HF methods (e.g., local MP2, local CC).
+
+    The nomenclature for these thresholds is adapted from Molpro, ORCA, etc.
+    See also:
+      - H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz,
+        "Molpro: a general-purpose quantum chemistry program package," 
+        WIREs Comput. Mol. Sci. 2, 242 (2012).
+      - F. Neese, "The ORCA program system," 
+        WIREs Comput. Mol. Sci. 2, 73-78 (2012) for DLPNO expansions.
+      - G. Knizia, "Intrinsic Bond Orbitals," 
+        J. Chem. Theory Comput. 9, 4834 (2013) for orbital transformations in local correlation.
     """
 
     domain_connectivity = Quantity(

From b0eb159711a1f9cb66b568f06576ca89d0614977 Mon Sep 17 00:00:00 2001
From: EBB2675 <esmaboydas4@gmail.com>
Date: Tue, 28 Jan 2025 14:42:15 +0100
Subject: [PATCH 37/37] minor modifications

---
 .../schema_packages/atoms_state.py            |  30 ++--
 .../schema_packages/basis_set.py              |  30 ++--
 .../schema_packages/model_method.py           |  15 +-
 .../schema_packages/numerical_settings.py     | 129 +++++-------------
 4 files changed, 69 insertions(+), 135 deletions(-)

diff --git a/src/nomad_simulations/schema_packages/atoms_state.py b/src/nomad_simulations/schema_packages/atoms_state.py
index b2b97b06..790df1e7 100644
--- a/src/nomad_simulations/schema_packages/atoms_state.py
+++ b/src/nomad_simulations/schema_packages/atoms_state.py
@@ -660,20 +660,6 @@ class MolecularOrbitals(Entity):
 
     """
 
-    mo_type = Quantity(
-        type=str,
-        shape=['mo_num'],
-        description="""
-        Type of the molecular orbitals
-         e.g. 'canonical', 'localized'.
-        In case of CASSCF calculations, there will be orbital subspaces of different nature.
-        E.g. : 
-        Internal orbitals : canonical
-        Active orbitals   : natural
-        Virtual orbitals  : canonical
-        """,
-    )
-
     mo_num = Quantity(
         type=np.int32,
         description="""
@@ -689,6 +675,20 @@ class MolecularOrbitals(Entity):
         """,
     )
 
+    mo_type = Quantity(
+        type=str,
+        shape=['mo_num'],
+        description="""
+        Type of the molecular orbitals
+         e.g. 'canonical', 'localized'.
+        In case of CASSCF calculations, there will be orbital subspaces of different nature.
+        E.g. : 
+        Internal orbitals : canonical
+        Active orbitals   : natural
+        Virtual orbitals  : canonical
+        """,
+    )
+
     coefficient = Quantity(
         type=np.float64,
         shape=['mo_num', 'ao_num'],
@@ -743,5 +743,3 @@ class MolecularOrbitals(Entity):
         Typically 0 for alpha, 1 for beta.
         """,
     )
-
-
diff --git a/src/nomad_simulations/schema_packages/basis_set.py b/src/nomad_simulations/schema_packages/basis_set.py
index b66d2f50..536a39dd 100644
--- a/src/nomad_simulations/schema_packages/basis_set.py
+++ b/src/nomad_simulations/schema_packages/basis_set.py
@@ -187,18 +187,18 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 
 class AtomCenteredFunction(ArchiveSection):
     """
-    Specifies a single contracted basis function in an atom-centered basis set. 
-    
-    In many quantum-chemistry codes, an atom-centered basis set is composed of 
-    several "shells," each shell containing one or more basis functions of a certain 
-    angular momentum. For instance, a shell of p-type orbitals (L=1) typically 
+    Specifies a single contracted basis function in an atom-centered basis set.
+
+    In many quantum-chemistry codes, an atom-centered basis set is composed of
+    several "shells," each shell containing one or more basis functions of a certain
+    angular momentum. For instance, a shell of p-type orbitals (L=1) typically
     consists of 3 degenerate functions (p_x, p_y, p_z) if `harmonic_type='cartesian'`
     or 3 spherical harmonics if `harmonic_type='spherical'`.
-    
-    A single "atom-centered function" can be a linear combination of multiple 
-    primitive Gaussians (or Slater-type orbitals, STOs). 
-    In practice, these contract together to form the final basis function used by 
-    the SCF or post-SCF method. Often, each contraction is labeled by its 
+
+    A single "atom-centered function" can be a linear combination of multiple
+    primitive Gaussians (or Slater-type orbitals, STOs).
+    In practice, these contract together to form the final basis function used by
+    the SCF or post-SCF method. Often, each contraction is labeled by its
     angular momentum (e.g., s, p, d, f) and a set of exponents and coefficients.
 
     **References**:
@@ -353,16 +353,16 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 class AtomCenteredBasisSet(BasisSetComponent):
     """
     Defines an **atom-centered basis set** for quantum chemistry calculations.
-    Unlike plane-wave methods, these expansions are typically built around each atom's 
+    Unlike plane-wave methods, these expansions are typically built around each atom's
     position, using either:
       - Slater-type orbitals (STO)
       - Gaussian-type orbitals (GTO)
       - Numerical atomic orbitals (NAO)
       - Effective-core potentials or point-charges (PC, cECP, etc.)
-    
+
     This section references multiple `AtomCenteredFunction` objects, each describing
-    a single contracted function or shell. Additionally, one can label the overall 
-    basis set name (e.g., "cc-pVTZ", "def2-SVP", "6-31G**") and specify the high-level 
+    a single contracted function or shell. Additionally, one can label the overall
+    basis set name (e.g., "cc-pVTZ", "def2-SVP", "6-31G**") and specify the high-level
     role of the basis set in the calculation.
 
     **Common examples**:
@@ -411,7 +411,7 @@ class AtomCenteredBasisSet(BasisSetComponent):
             'orbital',
             'auxiliary_scf',
             'auxiliary_post_hf',
-            'cabs',  
+            'cabs',
         ),
         description="""
         The role of this basis set in the calculation:
diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index d3055ffc..c0e20aa1 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -1272,17 +1272,17 @@ class IntegralDecomposition(BaseModelMethod):
 
 class HartreeFock(ModelMethodElectronic):
     """
-    Defines a Hartree–Fock (HF) calculation using a specified reference determinant 
+    Defines a Hartree–Fock (HF) calculation using a specified reference determinant
     (RHF, UHF, or ROHF).
 
     In HF theory:
       - RHF  = Restricted Hartree–Fock, for closed-shell systems.
       - UHF  = Unrestricted Hartree–Fock, allows different orbitals for alpha/beta spin.
-      - ROHF = Restricted Open-Shell Hartree–Fock, a partially restricted approach for 
+      - ROHF = Restricted Open-Shell Hartree–Fock, a partially restricted approach for
                open-shell systems.
 
     **References**:
-      - Roothaan, C. C. J. (1951). "New Developments in Molecular Orbital Theory." 
+      - Roothaan, C. C. J. (1951). "New Developments in Molecular Orbital Theory."
         Rev. Mod. Phys. 23, 69.
       - Szabo, A., & Ostlund, N. S. (1989). *Modern Quantum Chemistry*. McGraw-Hill.
       - Jensen, F. (2007). *Introduction to Computational Chemistry*. 2nd ed., Wiley.
@@ -1306,7 +1306,7 @@ class HartreeFock(ModelMethodElectronic):
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         """
         Perform minimal consistency checks between the HF reference determinant
-        and the final molecular orbitals spin array (if available).
+        and the final molecular orbitals spin array (if available).NumericalIntegration
         """
         super().normalize(archive, logger)
 
@@ -1319,15 +1319,15 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
                     # For RHF, we typically expect spin=0 (alpha) only
                     if not np.all(mo_spin == 0):
                         logger.warning(
-                            "RHF reference used, but molecular_orbitals.spin contains non-zero spin indices."
+                            'RHF reference used, but molecular_orbitals.spin contains non-zero spin indices.'
                         )
                 elif self.reference_determinant == 'UHF':
                     # UHF often has spin=0 for alpha and spin=1 for beta
                     # If we only see spin=0, that's effectively no spin polarization
                     if len(unique_spins) == 1 and unique_spins[0] == 0:
                         logger.info(
-                            "UHF reference chosen, but only alpha spin found in MOs (spin=0). "
-                            "This might still be valid if spin polarization is zero."
+                            'UHF reference chosen, but only alpha spin found in MOs (spin=0). '
+                            'This might still be valid if spin polarization is zero.'
                         )
                 # For ROHF, spin indexing can vary across codes, so no strict check here.
 
@@ -1376,7 +1376,6 @@ class PerturbationMethod(ModelMethodElectronic):
     )
 
 
-
 class LocalCorrelation(ArchiveSection):
     """
     A base section used to define the parameters of a local correlation method for
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index eca33e56..b43a2177 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -955,7 +955,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
 class OrbitalLocalization(SelfConsistency):
     """
     Numerical settings that control orbital localization, typically applied
-    to transform canonical molecular orbitals into localized orbitals. 
+    to transform canonical molecular orbitals into localized orbitals.
     These localized orbitals can then be used for:
       - Local correlation methods (e.g., LMP2, local CC)
       - Interpretable chemical analysis (e.g., identifying bonds, lone pairs)
@@ -965,10 +965,10 @@ class OrbitalLocalization(SelfConsistency):
     iterative, requiring thresholds and iteration limits akin to an SCF loop.
 
     References:
-      - R. F. Boys, "Construction of Some Molecular Orbitals to Be Approximately Invariant 
+      - R. F. Boys, "Construction of Some Molecular Orbitals to Be Approximately Invariant
         for Changes in Molecular Conformation," Rev. Mod. Phys. 32, 296 (1960). [Boys]
       - J. Pipek, P. G. Mezey, J. Chem. Phys. 90, 4916 (1989). [PM method]
-      - J. L. Knizia, "Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory 
+      - J. L. Knizia, "Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory
         and Chemical Concepts," J. Chem. Theory Comput. 9, 4834 (2013). [IBO method]
       - P. Pulay, Chem. Phys. Lett. 100, 151 (1983). [Local correlation context]
     """
@@ -1010,113 +1010,50 @@ class OrbitalLocalization(SelfConsistency):
     )
 
 
-class PNOSettings(NumericalSettings):
+class PairApproximationSettings(NumericalSettings):
     """
-    Numerical settings that control **Pair Natural Orbitals (PNOs)** in local correlation approaches.
+    Numerical settings that control Pair Natural Orbitals (PNOs) in local correlation approaches.
     PNOs are a compact representation of the virtual orbital space for each pair (or domain)
     of occupied orbitals, improving efficiency in post-HF methods (e.g., local MP2, local CC).
 
+    The nomenclature for these thresholds is adapted from different programs (e.g., Molpro, ORCA).
+    This class is code-agnostic and allows a single 'approximation_level' to store
+    code-specific presets such as:
+    'tightPNO', 'normalPNO', 'loosePNO' (ORCA) or
+    'default', 'tight', 'vtight' (Molpro), etc.
+
+    For more fine-grained custom thresholds, the user can fill the other
+    domain/pair fields (occupancies, energy thresholds, etc.) on the parser side as needed.
+
     The nomenclature for these thresholds is adapted from Molpro, ORCA, etc.
     See also:
       - H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz,
-        "Molpro: a general-purpose quantum chemistry program package," 
+        "Molpro: a general-purpose quantum chemistry program package,"
         WIREs Comput. Mol. Sci. 2, 242 (2012).
-      - F. Neese, "The ORCA program system," 
+      - F. Neese, "The ORCA program system,"
         WIREs Comput. Mol. Sci. 2, 73-78 (2012) for DLPNO expansions.
-      - G. Knizia, "Intrinsic Bond Orbitals," 
+      - G. Knizia, "Intrinsic Bond Orbitals,"
         J. Chem. Theory Comput. 9, 4834 (2013) for orbital transformations in local correlation.
     """
 
-    domain_connectivity = Quantity(
-        type=int,
-        description="""
-        the connectivity parameter for domain extension.
-        """,
-    )
-
-    domain_radius = Quantity(
-        type=int,
-        description="""
-        the radius parameter for domain extension.
-        """,
-    )
-
-    t_domain_osv_occ = Quantity(
-        type=np.float32,
-        description="""
-        OSV domain occupation number threshold.
-        """,
-    )
-
-    t_occ_lmp2 = Quantity(
-        type=np.float32,
-        description="""
-        LMP2 PNO domains (occ. number threshold).
-        """,
-    )
-
-    t_energy_lmp2 = Quantity(
-        type=np.float32,
-        description="""
-        LMP2 PNO domains (energy threshold).
-        """,
-    )
-
-    t_occ_lccsd = Quantity(
-        type=np.float32,
-        description="""
-        LCCSD PNO domains (occ. number threshold).
-        """,
-    )
-
-    t_energy_lccsd = Quantity(
-        type=np.float32,
-        description="""
-        LCCSD PNO domains (energy threshold).
-        """,
-    )
+    # type = Quantity(
+    #     type=MEnum('PAO', 'OSV', 'PNO'),
+    #     description="""
+    #     PAO : pair atomic orbitals
+    #     OSV : orbital-specific virtuals
+    #     PNO : pair natural orbitals
+    #     """,
+    # )
 
-    t_close_pair = Quantity(
+    code_specific_approximation_tier = Quantity(
         type=str,
         description="""
-        close pair energy threshold.
-        """,
-    )
-
-    t_weak_pair = Quantity(
-        type=np.float32,
-        description="""
-        weak pair energy threshold.
-        """,
-    )
-
-    t_distant_pair = Quantity(
-        type=np.float32,
-        description="""
-        distant pair energy threshold
-        """,
-    )
-
-    t_verydistant_pair = Quantity(
-        type=np.float32,
-        description="""
-        very distant pair energy threshold
+        A single keyword or label indicating a preset PNO approximation level. Examples:
+         - In ORCA: 'loosePNO', 'normalPNO', 'tightPNO'
+         - In Molpro: 'default', 'tight', 'vtight'
+         - In Psi4: 'loose', 'normal', 'tight' for PNO_CONVERGENCE keyword
+         - Or any other custom label recognized by your code
+        
+        This field is purely descriptive to unify different codes' preset naming.
         """,
     )
-
-    t_triples_preselection = Quantity(
-        type=np.float32,
-        description="""
-        preselection of triples list.
-        """,
-    )
-
-    t_triples_iteration = Quantity(
-        type=np.float32,
-        description="""
-        selection of triples for iterations
-        """,
-    )
-
-    def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
-        super().normalize(archive, logger)