From 54e0778dc49b7e03a30ba7cb7b5bcf94f0746ab5 Mon Sep 17 00:00:00 2001
From: Toon Verstraelen <toon.verstraelen@UGent.be>
Date: Sun, 12 May 2019 19:59:59 +0200
Subject: [PATCH 1/3] Move MolecularOrbitals to its own module

---
 iodata/formats/cp2k.py    |  3 +-
 iodata/formats/fchk.py    |  3 +-
 iodata/formats/molden.py  |  3 +-
 iodata/formats/molekel.py |  3 +-
 iodata/formats/wfn.py     |  3 +-
 iodata/orbitals.py        | 83 +++++++++++++++++++++++++++++++++++++++
 iodata/utils.py           | 59 +---------------------------
 7 files changed, 95 insertions(+), 62 deletions(-)
 create mode 100644 iodata/orbitals.py

diff --git a/iodata/formats/cp2k.py b/iodata/formats/cp2k.py
index e732c357..067fb8ac 100644
--- a/iodata/formats/cp2k.py
+++ b/iodata/formats/cp2k.py
@@ -25,7 +25,8 @@
 from scipy.special import factorialk
 
 from ..basis import angmom_sti, MolecularBasis, Shell, HORTON2_CONVENTIONS
-from ..utils import LineIterator, MolecularOrbitals
+from ..orbitals import MolecularOrbitals
+from ..utils import LineIterator
 
 
 __all__ = []
diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py
index 54fb6064..2f70f402 100644
--- a/iodata/formats/fchk.py
+++ b/iodata/formats/fchk.py
@@ -25,7 +25,8 @@
 import numpy as np
 
 from ..basis import MolecularBasis, Shell, HORTON2_CONVENTIONS
-from ..utils import LineIterator, MolecularOrbitals, amu
+from ..orbitals import MolecularOrbitals
+from ..utils import LineIterator, amu
 
 
 __all__ = []
diff --git a/iodata/formats/molden.py b/iodata/formats/molden.py
index 39b5d19c..6e6418c8 100644
--- a/iodata/formats/molden.py
+++ b/iodata/formats/molden.py
@@ -28,8 +28,9 @@
                      convert_conventions, HORTON2_CONVENTIONS)
 from ..iodata import IOData
 from ..periodic import sym2num, num2sym
+from ..orbitals import MolecularOrbitals
 from ..overlap import compute_overlap, gob_cart_normalization
-from ..utils import angstrom, LineIterator, MolecularOrbitals
+from ..utils import angstrom, LineIterator
 
 
 __all__ = []
diff --git a/iodata/formats/molekel.py b/iodata/formats/molekel.py
index 5f6af966..3acd2bc4 100644
--- a/iodata/formats/molekel.py
+++ b/iodata/formats/molekel.py
@@ -25,7 +25,8 @@
 
 from .molden import CONVENTIONS, _fix_molden_from_buggy_codes
 from ..basis import angmom_sti, MolecularBasis, Shell
-from ..utils import angstrom, LineIterator, MolecularOrbitals
+from ..orbitals import MolecularOrbitals
+from ..utils import angstrom, LineIterator
 
 
 __all__ = []
diff --git a/iodata/formats/wfn.py b/iodata/formats/wfn.py
index 1e722210..2e75794a 100644
--- a/iodata/formats/wfn.py
+++ b/iodata/formats/wfn.py
@@ -26,7 +26,8 @@
 from ..basis import MolecularBasis, Shell
 from ..overlap import gob_cart_normalization
 from ..periodic import sym2num
-from ..utils import MolecularOrbitals, LineIterator
+from ..orbitals import MolecularOrbitals
+from ..utils import LineIterator
 
 
 __all__ = []
diff --git a/iodata/orbitals.py b/iodata/orbitals.py
new file mode 100644
index 00000000..6801c72f
--- /dev/null
+++ b/iodata/orbitals.py
@@ -0,0 +1,83 @@
+# IODATA is an input and output module for quantum chemistry.
+# Copyright (C) 2011-2019 The IODATA Development Team
+#
+# This file is part of IODATA.
+#
+# IODATA is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 3
+# of the License, or (at your option) any later version.
+#
+# IODATA is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, see <http://www.gnu.org/licenses/>
+# --
+"""Data structure for molecular orbitals."""
+
+
+from typing import NamedTuple
+
+import numpy as np
+
+
+__all__ = ['MolecularOrbitals']
+
+
+class MolecularOrbitals(NamedTuple):
+    """Molecular Orbitals Class.
+
+    Attributes
+    ----------
+    type : str
+        Molecular orbital type; choose from 'restricted', 'unrestricted', or 'generalized'.
+    norba : int
+        Number of alpha molecular orbitals. None in case of type=='generalized'.
+    norbb : int
+        Number of beta molecular orbitals. None in case of type=='generalized'.
+    occs : np.ndarray
+        Molecular orbital occupation numbers. The number of elements equals the
+        number of columns of coeffs.
+    coeffs : np.ndarray
+        Molecular orbital basis coefficients.
+        In case of restricted: shape = (nbasis, norb_a) = (nbasis, norb_b).
+        In case of unrestricted: shape = (nbasis, norb_a + norb_b).
+        In case of generalized: shape = (2*nbasis, norb), where norb is the
+        total number of orbitals (not defined by other attributes).
+    irreps : np.ndarray
+        Irreducible representation. The number of elements equals the
+        number of columns of coeffs.
+    energies : np.ndarray
+        Molecular orbital energies. The number of elements equals the
+        number of columns of coeffs.
+
+    """
+
+    type: str
+    norba: int
+    norbb: int
+    occs: np.ndarray
+    coeffs: np.ndarray
+    irreps: np.ndarray
+    energies: np.ndarray
+
+    @property
+    def nelec(self) -> float:
+        """Return the total number of electrons."""
+        return self.occs.sum()
+
+    @property
+    def spinpol(self) -> float:
+        """Return the spin multiplicity of the Slater determinant."""
+        if self.type == 'restricted':
+            nbeta = np.clip(self.occs, 0, 1).sum()
+            sq = self.nelec - 2 * nbeta
+        elif self.type == 'unrestricted':
+            sq = self.occs[:self.norba].sum() - self.occs[self.norba:].sum()
+        else:
+            # Not sure how to do this in a simply way.
+            raise NotImplementedError
+        return abs(sq)
diff --git a/iodata/utils.py b/iodata/utils.py
index 785b44fe..5bd34518 100644
--- a/iodata/utils.py
+++ b/iodata/utils.py
@@ -26,7 +26,8 @@
 from scipy.linalg import eigh
 
 
-__all__ = ['LineIterator', 'set_four_index_element', 'MolecularOrbitals']
+__all__ = ['LineIterator', 'Cube', 'set_four_index_element', 'volume',
+           'derive_naturals', 'check_dm']
 
 
 # The unit conversion factors below can be used as follows:
@@ -87,62 +88,6 @@ def back(self, line):
         self.lineno -= 1
 
 
-class MolecularOrbitals(NamedTuple):
-    """Molecular Orbitals Class.
-
-    Attributes
-    ----------
-    type : str
-        Molecular orbital type; choose from 'restricted', 'unrestricted', or 'generalized'.
-    norba : int
-        Number of alpha molecular orbitals. None in case of type=='generalized'.
-    norbb : int
-        Number of beta molecular orbitals. None in case of type=='generalized'.
-    occs : np.ndarray
-        Molecular orbital occupation numbers. The number of elements equals the
-        number of columns of coeffs.
-    coeffs : np.ndarray
-        Molecular orbital basis coefficients.
-        In case of restricted: shape = (nbasis, norb_a) = (nbasis, norb_b).
-        In case of unrestricted: shape = (nbasis, norb_a + norb_b).
-        In case of generalized: shape = (2*nbasis, norb), where norb is the
-        total number of orbitals (not defined by other attributes).
-    irreps : np.ndarray
-        Irreducible representation. The number of elements equals the
-        number of columns of coeffs.
-    energies : np.ndarray
-        Molecular orbital energies. The number of elements equals the
-        number of columns of coeffs.
-
-    """
-
-    type: str
-    norba: int
-    norbb: int
-    occs: np.ndarray
-    coeffs: np.ndarray
-    irreps: np.ndarray
-    energies: np.ndarray
-
-    @property
-    def nelec(self):
-        """Return the total number of electrons."""
-        return self.occs.sum()
-
-    @property
-    def spinpol(self):
-        """Return the spin multiplicity of the Slater determinant."""
-        if self.type == 'restricted':
-            nbeta = np.clip(self.occs, 0, 1).sum()
-            sq = self.nelec - 2 * nbeta
-        elif self.type == 'unrestricted':
-            sq = self.occs[:self.norba].sum() - self.occs[self.norba:].sum()
-        else:
-            # Not sure how to do this in a simply way.
-            raise NotImplementedError
-        return abs(sq)
-
-
 class Cube(NamedTuple):
     """The volumetric data from a cube (or similar) file.
 

From 0f66aad5802c1e7c1e0fd11e2b7a7e9031ed9e89 Mon Sep 17 00:00:00 2001
From: Toon Verstraelen <toon.verstraelen@UGent.be>
Date: Sun, 12 May 2019 21:36:37 +0200
Subject: [PATCH 2/3] Convenient alpha and beta properties of mos

---
 iodata/formats/cp2k.py      |  22 +++---
 iodata/formats/fchk.py      |  13 ++--
 iodata/formats/molden.py    | 131 ++++++++++++++++++-----------------
 iodata/formats/molekel.py   |  64 +++++++++--------
 iodata/formats/wfn.py       |  31 ++++-----
 iodata/iodata.py            |   7 +-
 iodata/orbitals.py          | 133 ++++++++++++++++++++++++++++++------
 iodata/test/common.py       |   2 +-
 iodata/test/test_cp2k.py    | 103 ++++++++++++----------------
 iodata/test/test_fchk.py    |  77 +++++++++------------
 iodata/test/test_molden.py  |   7 +-
 iodata/test/test_molekel.py |   4 +-
 iodata/test/test_wfn.py     |  18 +++--
 13 files changed, 339 insertions(+), 273 deletions(-)

diff --git a/iodata/formats/cp2k.py b/iodata/formats/cp2k.py
index 067fb8ac..6e980cc5 100644
--- a/iodata/formats/cp2k.py
+++ b/iodata/formats/cp2k.py
@@ -448,20 +448,17 @@ def load_one(lit: LineIterator) -> dict:
 
     # Turn orbital data into a MolecularOrbitals object.
     if restricted:
-        mo_type = 'restricted'
         norb, nel = _get_norb_nel(oe_alpha)
-        norb_alpha = norb_beta = norb
         assert nel % 2 == 0
         orb_alpha_coeffs = np.zeros([obasis.nbasis, norb])
         orb_alpha_energies = np.zeros(norb)
         orb_alpha_occs = np.zeros(norb)
         _fill_orbitals(orb_alpha_coeffs, orb_alpha_energies, orb_alpha_occs,
                        oe_alpha, coeffs_alpha, obasis, restricted)
-        mo_occs = orb_alpha_occs
-        mo_coeffs = orb_alpha_coeffs
-        mo_energy = orb_alpha_energies
+        mo = MolecularOrbitals(
+            'restricted', norb, norb, 2 * orb_alpha_occs, orb_alpha_coeffs,
+            orb_alpha_energies, None)
     else:
-        mo_type = 'unrestricted'
         norb_alpha = _get_norb_nel(oe_alpha)[0]
         norb_beta = _get_norb_nel(oe_beta)[0]
         assert norb_alpha == norb_beta
@@ -476,12 +473,13 @@ def load_one(lit: LineIterator) -> dict:
         _fill_orbitals(orb_beta_coeffs, orb_beta_energies, orb_beta_occs,
                        oe_beta, coeffs_beta, obasis, restricted)
 
-        mo_occs = np.concatenate((orb_alpha_occs, orb_beta_occs), axis=0)
-        mo_energy = np.concatenate((orb_alpha_energies, orb_beta_energies), axis=0)
-        mo_coeffs = np.concatenate((orb_alpha_coeffs, orb_beta_coeffs), axis=1)
-
-    # create a MO namedtuple
-    mo = MolecularOrbitals(mo_type, norb_alpha, norb_beta, mo_occs, mo_coeffs, None, mo_energy)
+        mo = MolecularOrbitals(
+            'unrestricted', norb_alpha, norb_beta,
+            np.concatenate((orb_alpha_occs, orb_beta_occs), axis=0),
+            np.concatenate((orb_alpha_coeffs, orb_beta_coeffs), axis=1),
+            np.concatenate((orb_alpha_energies, orb_beta_energies), axis=0),
+            None,
+        )
 
     result = {
         'obasis': obasis,
diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py
index 2f70f402..6825291a 100644
--- a/iodata/formats/fchk.py
+++ b/iodata/formats/fchk.py
@@ -191,31 +191,28 @@ def load_one(lit: LineIterator) -> dict:
 
     norba = fchk['Alpha Orbital Energies'].shape[0]
     mo_coeffs = np.copy(fchk['Alpha MO coefficients'].reshape(nbasis_indep, nbasis).T)
-    mo_energy = np.copy(fchk['Alpha Orbital Energies'])
+    mo_energies = np.copy(fchk['Alpha Orbital Energies'])
 
     if 'Beta Orbital Energies' in fchk:
         # unrestricted
-        mo_type = 'unrestricted'
         norbb = fchk['Beta Orbital Energies'].shape[0]
         mo_coeffs_b = np.copy(fchk['Beta MO coefficients'].reshape(nbasis_indep, nbasis).T)
         mo_coeffs = np.concatenate((mo_coeffs, mo_coeffs_b), axis=1)
-        mo_energy = np.concatenate((mo_energy, np.copy(fchk['Beta Orbital Energies'])), axis=0)
+        mo_energies = np.concatenate((mo_energies, np.copy(fchk['Beta Orbital Energies'])), axis=0)
         mo_occs = np.zeros(2 * nbasis_indep)
         mo_occs[:nalpha] = 1.0
         mo_occs[nbasis_indep: nbasis_indep + nbeta] = 1.0
+        mo = MolecularOrbitals('unrestricted', norba, norbb, mo_occs, mo_coeffs, mo_energies, None)
     else:
         # restricted closed-shell and open-shell
-        mo_type = 'restricted'
-        norbb = norba
         mo_occs = np.zeros(nbasis_indep)
         mo_occs[:nalpha] = 1.0
         mo_occs[:nbeta] = 2.0
         if nalpha != nbeta:
             # delete dm_full_scf because it is known to be buggy
             result['one_rdms'].pop('scf')
-
-    # create a MO namedtuple
-    result['mo'] = MolecularOrbitals(mo_type, norba, norbb, mo_occs, mo_coeffs, None, mo_energy)
+        mo = MolecularOrbitals('restricted', norba, norba, mo_occs, mo_coeffs, mo_energies, None)
+    result['mo'] = mo
 
     # E) Load properties
     if 'Polarizability' in fchk:
diff --git a/iodata/formats/molden.py b/iodata/formats/molden.py
index 6e6418c8..a94b0844 100644
--- a/iodata/formats/molden.py
+++ b/iodata/formats/molden.py
@@ -104,12 +104,12 @@ def _load_low(lit: LineIterator) -> dict:
     atnums = None
     atcoords = None
     obasis = None
-    coeff_alpha = None
-    ener_alpha = None
-    occ_alpha = None
-    coeff_beta = None
-    ener_beta = None
-    occ_beta = None
+    coeffsa = None
+    energiesa = None
+    occsa = None
+    coeffsb = None
+    energiesb = None
+    occsb = None
     title = None
 
     line = next(lit)
@@ -154,8 +154,8 @@ def _load_low(lit: LineIterator) -> dict:
         # molecular-orbital coefficients.
         elif line == '[mo]':
             data_alpha, data_beta = _load_helper_coeffs(lit)
-            coeff_alpha, ener_alpha, occ_alpha = data_alpha
-            coeff_beta, ener_beta, occ_beta = data_beta
+            coeffsa, energiesa, occsa = data_alpha
+            coeffsb, energiesb, occsb = data_beta
 
     # Assign pure and Cartesian correctly. This needs to be done after reading
     # because the tags for pure functions may come after the basis set.
@@ -164,25 +164,21 @@ def _load_low(lit: LineIterator) -> dict:
         if shell.angmoms[0] in pure_angmoms:
             shell.kinds[0] = 'p'
 
-    if coeff_beta is None:
-        mo_type = 'restricted'
-        if coeff_alpha.shape[0] != obasis.nbasis:
+    if coeffsb is None:
+        if coeffsa.shape[0] != obasis.nbasis:
             lit.error("Number of alpha orbital coefficients does not match the size of the basis.")
-        norba = norbb = coeff_alpha.shape[1]
-        mo_occs = occ_alpha
-        mo_coeffs = coeff_alpha
-        mo_energy = ener_alpha
+        mo = MolecularOrbitals(
+            'restricted', coeffsa.shape[1], coeffsa.shape[1],
+            occsa, coeffsa, energiesa, None)
     else:
-        mo_type = 'unrestricted'
-        if coeff_beta.shape[0] != obasis.nbasis:
+        if coeffsb.shape[0] != obasis.nbasis:
             lit.error("Number of beta orbital coefficients does not match the size of the basis.")
-        norba = coeff_alpha.shape[1]
-        norbb = coeff_beta.shape[1]
-        mo_occs = np.concatenate((occ_alpha, occ_beta), axis=0)
-        mo_energy = np.concatenate((ener_alpha, ener_beta), axis=0)
-        mo_coeffs = np.concatenate((coeff_alpha, coeff_beta), axis=1)
-    # create a MO namedtuple
-    mo = MolecularOrbitals(mo_type, norba, norbb, mo_occs, mo_coeffs, None, mo_energy)
+        mo = MolecularOrbitals(
+            'unrestricted', coeffsa.shape[1], coeffsb.shape[1],
+            np.concatenate((occsa, occsb), axis=0),
+            np.concatenate((coeffsa, coeffsb), axis=1),
+            np.concatenate((energiesa, energiesb), axis=0),
+            None)
 
     result = {
         'atcoords': atcoords,
@@ -253,12 +249,12 @@ def _load_helper_obasis(lit: LineIterator) -> MolecularBasis:
 
 def _load_helper_coeffs(lit: LineIterator) -> Tuple:
     """Load the orbital coefficients."""
-    coeff_alpha = []
-    ener_alpha = []
-    occ_alpha = []
-    coeff_beta = []
-    ener_beta = []
-    occ_beta = []
+    coeffsa = []
+    energiesa = []
+    occsa = []
+    coeffsb = []
+    energiesb = []
+    occsb = []
 
     while True:
         try:
@@ -290,13 +286,13 @@ def _load_helper_coeffs(lit: LineIterator) -> Tuple:
         col = []
         # store column of coefficients, i.e. one orbital, energy and occ
         if info['spin'].strip().lower() == 'alpha':
-            coeff_alpha.append(col)
-            ener_alpha.append(energy)
-            occ_alpha.append(occ)
+            coeffsa.append(col)
+            energiesa.append(energy)
+            occsa.append(occ)
         else:
-            coeff_beta.append(col)
-            ener_beta.append(energy)
-            occ_beta.append(occ)
+            coeffsb.append(col)
+            energiesb.append(energy)
+            occsb.append(occ)
         for line in lit:
             words = line.split()
             if len(words) != 2 or not words[0].isdigit():
@@ -306,18 +302,18 @@ def _load_helper_coeffs(lit: LineIterator) -> Tuple:
                 break
             col.append(float(words[1]))
 
-    coeff_alpha = np.array(coeff_alpha).T
-    ener_alpha = np.array(ener_alpha)
-    occ_alpha = np.array(occ_alpha)
-    if not coeff_beta:
-        coeff_beta = None
-        ener_beta = None
-        occ_beta = None
+    coeffsa = np.array(coeffsa).T
+    energiesa = np.array(energiesa)
+    occsa = np.array(occsa)
+    if not coeffsb:
+        coeffsb = None
+        energiesb = None
+        occsb = None
     else:
-        coeff_beta = np.array(coeff_beta).T
-        ener_beta = np.array(ener_beta)
-        occ_beta = np.array(occ_beta)
-    return (coeff_alpha, ener_alpha, occ_alpha), (coeff_beta, ener_beta, occ_beta)
+        coeffsb = np.array(coeffsb).T
+        energiesb = np.array(energiesb)
+        occsb = np.array(occsb)
+    return (coeffsa, energiesa, occsa), (coeffsb, energiesb, occsb)
 
 
 def _is_normalized_properly(obasis: MolecularBasis, atcoords: np.ndarray,
@@ -523,15 +519,17 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
     """
     obasis = result['obasis']
     atcoords = result['atcoords']
-    if result['mo'].type == 'restricted':
-        coeffs_a = result['mo'].coeffs
-        coeffs_b = None
-    elif result['mo'].type == 'unrestricted':
-        coeffs_a = result['mo'].coeffs[:, :result['mo'].norba]
-        coeffs_b = result['mo'].coeffs[:, result['mo'].norba:]
+    if result['mo'].kind == 'restricted':
+        coeffsa = result['mo'].coeffs
+        # Skip testing coeffsb if it is the same array as coeffsa.
+        coeffsb = None
+    elif result['mo'].kind == 'unrestricted':
+        coeffsa = result['mo'].coeffsa
+        coeffsb = result['mo'].coeffsb
     else:
-        raise ValueError('Molecular orbital type={0} not recognized'.format(result['mo'].type))
-    if _is_normalized_properly(obasis, atcoords, coeffs_a, coeffs_b):
+        raise ValueError('Molecular orbital kind={0} not recognized'.format(result['mo'].kind))
+
+    if _is_normalized_properly(obasis, atcoords, coeffsa, coeffsb):
         # The file is good. No need to change obasis.
         return
     print('5:Detected incorrect normalization of orbitals loaded from a Molden or MKL file.')
@@ -539,7 +537,7 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
     # --- ORCA
     print('5:Trying to fix it as if it was a file generated by ORCA.')
     orca_obasis = _fix_obasis_orca(obasis)
-    if _is_normalized_properly(orca_obasis, atcoords, coeffs_a, coeffs_b):
+    if _is_normalized_properly(orca_obasis, atcoords, coeffsa, coeffsb):
         print('5:Detected typical ORCA errors in file. Fixing them...')
         result['obasis'] = orca_obasis
         return
@@ -548,7 +546,7 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
     print('5:Trying to fix it as if it was a file generated by PSI4 (pre 1.0).')
     psi4_obasis = _fix_obasis_psi4(obasis)
     if psi4_obasis is not None and \
-       _is_normalized_properly(psi4_obasis, atcoords, coeffs_a, coeffs_b):
+       _is_normalized_properly(psi4_obasis, atcoords, coeffsa, coeffsb):
         print('5:Detected typical PSI4 errors in file. Fixing them...')
         result['obasis'] = psi4_obasis
         return
@@ -557,7 +555,7 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
     print('5:Trying to fix it as if it was a file generated by Turbomole.')
     turbom_obasis = _fix_obasis_turbomole(obasis)
     if turbom_obasis is not None and \
-       _is_normalized_properly(turbom_obasis, atcoords, coeffs_a, coeffs_b):
+       _is_normalized_properly(turbom_obasis, atcoords, coeffsa, coeffsb):
         print('5:Detected typical Turbomole errors in file. Fixing them...')
         result['obasis'] = turbom_obasis
         return
@@ -565,7 +563,7 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
     # --- Renormalized contractions
     print('5:Last resort: trying by renormalizing all contractions')
     normed_obasis = _fix_obasis_normalize_contractions(obasis)
-    if _is_normalized_properly(normed_obasis, atcoords, coeffs_a, coeffs_b):
+    if _is_normalized_properly(normed_obasis, atcoords, coeffsa, coeffsb):
         print('5:Detected unnormalized contractions in file. Fixing them...')
         result['obasis'] = normed_obasis
         return
@@ -664,17 +662,18 @@ def dump_one(f: TextIO, data: IOData):
     permutation, signs = convert_conventions(obasis, CONVENTIONS)
 
     # Print the mean-field orbitals
-    if data.mo.type == 'unrestricted':
+    if data.mo.kind == 'unrestricted':
         f.write('[MO]\n')
-        norba = data.mo.norba
-        _dump_helper_orb(f, 'Alpha', data.mo.energies[:norba], data.mo.occs[:norba],
-                         data.mo.coeffs[:, :norba][permutation] * signs.reshape(-1, 1))
-        _dump_helper_orb(f, 'Beta', data.mo.energies[norba:], data.mo.occs[norba:],
-                         data.mo.coeffs[:, norba:][permutation] * signs.reshape(-1, 1))
-    else:
+        _dump_helper_orb(f, 'Alpha', data.mo.energiesa, data.mo.occsa,
+                         data.mo.coeffsa[permutation] * signs.reshape(-1, 1))
+        _dump_helper_orb(f, 'Beta', data.mo.energiesb, data.mo.occsb,
+                         data.mo.coeffsb[permutation] * signs.reshape(-1, 1))
+    elif data.mo.kind == 'restricted':
         f.write('[MO]\n')
         _dump_helper_orb(f, 'Alpha', data.mo.energies, data.mo.occs,
                          data.mo.coeffs[permutation] * signs.reshape(-1, 1))
+    else:
+        raise NotImplementedError
 
 
 def _dump_helper_orb(f, spin, orb_energies, orb_occs, orb_coeffs):
diff --git a/iodata/formats/molekel.py b/iodata/formats/molekel.py
index 3acd2bc4..2fdb59d8 100644
--- a/iodata/formats/molekel.py
+++ b/iodata/formats/molekel.py
@@ -162,12 +162,12 @@ def load_one(lit: LineIterator) -> dict:
     atnums = None
     atcoords = None
     obasis = None
-    coeff_alpha = None
-    ener_alpha = None
-    occ_alpha = None
-    coeff_beta = None
-    ener_beta = None
-    occ_beta = None
+    coeffsa = None
+    energiesa = None
+    occsa = None
+    coeffsb = None
+    energiesb = None
+    occsb = None
     # Using a loop because we're not entirely sure if sections in an MKL file
     # have a fixed order.
     while True:
@@ -184,13 +184,13 @@ def load_one(lit: LineIterator) -> dict:
         elif line == '$BASIS':
             obasis = _load_helper_obasis(lit)
         elif line == '$COEFF_ALPHA':
-            coeff_alpha, ener_alpha = _load_helper_coeffs(lit, obasis.nbasis)
+            coeffsa, energiesa = _load_helper_coeffs(lit, obasis.nbasis)
         elif line == '$OCC_ALPHA':
-            occ_alpha = _load_helper_occ(lit)
+            occsa = _load_helper_occ(lit)
         elif line == '$COEFF_BETA':
-            coeff_beta, ener_beta = _load_helper_coeffs(lit, obasis.nbasis)
+            coeffsb, energiesb = _load_helper_coeffs(lit, obasis.nbasis)
         elif line == '$OCC_BETA':
-            occ_beta = _load_helper_occ(lit)
+            occsb = _load_helper_occ(lit)
 
     if charge is None:
         lit.error('Charge and spin polarization not found.')
@@ -198,40 +198,38 @@ def load_one(lit: LineIterator) -> dict:
         lit.error('Coordinates not found.')
     if obasis is None:
         lit.error('Orbital basis not found.')
-    if coeff_alpha is None:
+    if coeffsa is None:
         lit.error('Alpha orbitals not found.')
-    if occ_alpha is None:
+    if occsa is None:
         lit.error('Alpha occupation numbers not found.')
 
     nelec = atnums.sum() - charge
-    if coeff_beta is None:
+    if coeffsb is None:
         # restricted closed-shell
-        mo_type = 'restricted'
         assert nelec % 2 == 0
-        assert abs(occ_alpha.sum() - nelec) < 1e-7
-        norba = norbb = coeff_alpha.shape[1]
-        mo_occs = occ_alpha
-        mo_coeffs = coeff_alpha
-        mo_energy = ener_alpha
+        assert abs(occsa.sum() - nelec) < 1e-7
+        mo = MolecularOrbitals(
+            'restricted', coeffsa.shape[1], coeffsa.shape[1],
+            occsa, coeffsa, energiesa, None)
     else:
-        mo_type = 'unrestricted'
-        if occ_beta is None:
+        if occsb is None:
             lit.error('Beta occupation numbers not found in mkl file while '
                       'beta orbitals were present.')
-        nalpha = int(np.round(occ_alpha.sum()))
-        nbeta = int(np.round(occ_beta.sum()))
+        nalpha = int(np.round(occsa.sum()))
+        nbeta = int(np.round(occsb.sum()))
         assert abs(spinpol - abs(nalpha - nbeta)) < 1e-7
         assert nelec == nalpha + nbeta
-        assert coeff_alpha.shape == coeff_beta.shape
-        assert ener_alpha.shape == ener_beta.shape
-        assert occ_alpha.shape == occ_beta.shape
-        norba = coeff_alpha.shape[1]
-        norbb = coeff_beta.shape[1]
-        mo_occs = np.concatenate((occ_alpha, occ_beta), axis=0)
-        mo_energy = np.concatenate((ener_alpha, ener_beta), axis=0)
-        mo_coeffs = np.concatenate((coeff_alpha, coeff_beta), axis=1)
-    # create a MO namedtuple
-    mo = MolecularOrbitals(mo_type, norba, norbb, mo_occs, mo_coeffs, None, mo_energy)
+        assert coeffsa.shape == coeffsb.shape
+        assert energiesa.shape == energiesb.shape
+        assert occsa.shape == occsb.shape
+        mo = MolecularOrbitals(
+            'unrestricted',
+            coeffsa.shape[1],
+            coeffsb.shape[1],
+            np.concatenate((occsa, occsb), axis=0),
+            np.concatenate((coeffsa, coeffsb), axis=1),
+            np.concatenate((energiesa, energiesb), axis=0),
+            None)
 
     result = {
         'atcoords': atcoords,
diff --git a/iodata/formats/wfn.py b/iodata/formats/wfn.py
index 2e75794a..4c658e00 100644
--- a/iodata/formats/wfn.py
+++ b/iodata/formats/wfn.py
@@ -326,25 +326,24 @@ def load_one(lit: LineIterator) -> dict:
             scales.append(gob_cart_normalization(shell.exponents[0], np.array([nx, ny, nz])))
     scales = np.array(scales)
     mo_coefficients /= scales.reshape(-1, 1)
+    norb = mo_coefficients.shape[1]
     # make the wavefunction
     if mo_occ.max() > 1.0:
-        # close shell system
-        mo_type = 'restricted'
-        na_orb = len(mo_occ)
-        nb_orb = len(mo_occ)
+        # closed-shell system
+        mo = MolecularOrbitals(
+            'restricted', norb, norb,
+            mo_occ, mo_coefficients, mo_energy, None)
     else:
-        # open shell system
-        mo_type = 'unrestricted'
-        # counting the number of alpha and beta orbitals
-        n = 1
-        while (n < mo_coefficients.shape[1]
-               and mo_energy[n] >= mo_energy[n - 1]
-               and mo_count[n] == mo_count[n - 1] + 1):
-            n += 1
-        na_orb = n
-        nb_orb = len(mo_occ) - n
-    # create a MO namedtuple
-    mo = MolecularOrbitals(mo_type, na_orb, nb_orb, mo_occ, mo_coefficients, None, mo_energy)
+        # open-shell system
+        # counting the number of alpha orbitals
+        norba = 1
+        while (norba < mo_coefficients.shape[1]
+               and mo_energy[norba] >= mo_energy[norba - 1]
+               and mo_count[norba] == mo_count[norba - 1] + 1):
+            norba += 1
+        mo = MolecularOrbitals(
+            'unrestricted', norba, norb - norba,
+            mo_occ, mo_coefficients, mo_energy, None)
 
     result = {
         'title': title,
diff --git a/iodata/iodata.py b/iodata/iodata.py
index 5ee3fd0c..f4da0cd3 100644
--- a/iodata/iodata.py
+++ b/iodata/iodata.py
@@ -344,7 +344,12 @@ def charge(self, charge: float):
 
     @property
     def spinpol(self) -> float:
-        """Return the spin multiplicity."""
+        """Return the spin polarization.
+
+        Warning: for restricted wavefunctions, it is assumed that an occupation
+        number in ]0, 2[ implies spin polarizaiton, which may not always be a
+        valid assumption.
+        """
         mo = getattr(self, 'mo', None)
         if mo is None:
             return self._spinpol
diff --git a/iodata/orbitals.py b/iodata/orbitals.py
index 6801c72f..e60baf12 100644
--- a/iodata/orbitals.py
+++ b/iodata/orbitals.py
@@ -28,56 +28,147 @@
 
 
 class MolecularOrbitals(NamedTuple):
-    """Molecular Orbitals Class.
+    """Molecular Orbitals base Class.
 
     Attributes
     ----------
-    type : str
-        Molecular orbital type; choose from 'restricted', 'unrestricted', or 'generalized'.
+    kind
+        'restricted', 'unrestricted', 'generalized'
     norba : int
         Number of alpha molecular orbitals. None in case of type=='generalized'.
     norbb : int
         Number of beta molecular orbitals. None in case of type=='generalized'.
-    occs : np.ndarray
+    occs
         Molecular orbital occupation numbers. The number of elements equals the
         number of columns of coeffs.
-    coeffs : np.ndarray
+    coeffs
         Molecular orbital basis coefficients.
         In case of restricted: shape = (nbasis, norb_a) = (nbasis, norb_b).
         In case of unrestricted: shape = (nbasis, norb_a + norb_b).
         In case of generalized: shape = (2*nbasis, norb), where norb is the
-        total number of orbitals (not defined by other attributes).
-    irreps : np.ndarray
-        Irreducible representation. The number of elements equals the
-        number of columns of coeffs.
-    energies : np.ndarray
+        total number of orbitals.
+    energies
         Molecular orbital energies. The number of elements equals the
         number of columns of coeffs.
+    irreps
+        Irreducible representation. The number of elements equals the
+        number of columns of coeffs.
+
+    Warning: the interpretation of the occupation numbers may only be suitable
+    for single-reference orbitals (not fractionally occupied natural orbitals.)
+    When an occupation number is in ]0, 1], it is assumed that an alpha orbital
+    is (fractionally) occupied. When an occupation number is in ]1, 2], it is
+    assumed that the alpha orbital is fully occupied and the beta orbital is
+    (fractionally) occupied.
 
     """
 
-    type: str
+    kind: str
     norba: int
     norbb: int
     occs: np.ndarray
     coeffs: np.ndarray
-    irreps: np.ndarray
     energies: np.ndarray
+    irreps: np.ndarray
 
     @property
     def nelec(self) -> float:
         """Return the total number of electrons."""
         return self.occs.sum()
 
+    @property
+    def nbasis(self):
+        """Return the number of spatial basis functions."""
+        if self.kind == 'generalized':
+            return self.coeffs.shape[0] // 2
+        return self.coeffs.shape[0]
+
+    @property
+    def norb(self):
+        """Return the number of orbitals."""
+        return self.coeffs.shape[1]
+
     @property
     def spinpol(self) -> float:
-        """Return the spin multiplicity of the Slater determinant."""
-        if self.type == 'restricted':
+        """Return the spin polarization of the Slater determinant."""
+        if self.kind == 'restricted':
             nbeta = np.clip(self.occs, 0, 1).sum()
-            sq = self.nelec - 2 * nbeta
-        elif self.type == 'unrestricted':
-            sq = self.occs[:self.norba].sum() - self.occs[self.norba:].sum()
-        else:
-            # Not sure how to do this in a simply way.
-            raise NotImplementedError
-        return abs(sq)
+            return abs(self.nelec - 2 * nbeta)
+        if self.kind == 'unrestricted':
+            return abs(self.occsa.sum() - self.occsb.sum())
+        raise NotImplementedError
+
+    @property
+    def occsa(self):
+        """Return alpha occupation numbers."""
+        if self.kind == 'restricted':
+            return np.clip(self.occs, 0, 1)
+        if self.kind == 'unrestricted':
+            return self.occs[:self.norba]
+        raise NotImplementedError
+
+    @property
+    def occsb(self):
+        """Return beta occupation numbers."""
+        if self.kind == 'restricted':
+            return self.occs - np.clip(self.occs, 0, 1)
+        if self.kind == 'unrestricted':
+            return self.occs[self.norba:]
+        raise NotImplementedError
+
+    @property
+    def coeffsa(self):
+        """Return alpha orbital coefficients."""
+        if self.kind == 'restricted':
+            return self.coeffs
+        if self.kind == 'unrestricted':
+            return self.coeffs[:, :self.norba]
+        raise NotImplementedError
+
+    @property
+    def coeffsb(self):
+        """Return beta orbital coefficients."""
+        if self.kind == 'restricted':
+            return self.coeffs
+        if self.kind == 'unrestricted':
+            return self.coeffs[:, self.norba:]
+        raise NotImplementedError
+
+    @property
+    def irrepsa(self):
+        """Return alpha irreps."""
+        if self.kind == 'restricted':
+            return self.irreps
+        if self.kind == 'unrestricted':
+            return self.irreps[:self.norba]
+        raise NotImplementedError
+
+    @property
+    def irrepsb(self):
+        """Return beta irreps."""
+        if self.kind == 'restricted':
+            return self.irreps
+        if self.kind == 'unrestricted':
+            return self.irreps[self.norba:]
+        raise NotImplementedError
+
+    @property
+    def energiesa(self):
+        """Return alpha orbital energies."""
+        if self.kind == 'restricted':
+            return self.energies
+        if self.kind == 'unrestricted':
+            return self.energies[:self.norba]
+        raise NotImplementedError
+
+    @property
+    def energiesb(self):
+        """Return beta orbital energies."""
+        if self.kind == 'restricted':
+            return self.energies
+        if self.kind == 'unrestricted':
+            return self.energies[self.norba:]
+        raise NotImplementedError
+
+
+MolecularOrbitals.__defaults__ = (None,)
diff --git a/iodata/test/common.py b/iodata/test/common.py
index 1f29965f..5e7f05e8 100644
--- a/iodata/test/common.py
+++ b/iodata/test/common.py
@@ -109,7 +109,7 @@ def compare_mols(mol1, mol2):
         assert mol2.obasis is None
     # wfn
     permutation, signs = convert_conventions(mol1.obasis, mol2.obasis.conventions)
-    assert mol1.mo.type == mol2.mo.type
+    assert mol1.mo.kind == mol2.mo.kind
     assert_allclose(mol1.mo.occs, mol2.mo.occs)
     assert_allclose(mol1.mo.energies, mol2.mo.energies)
     assert_allclose(mol1.mo.coeffs[permutation] * signs.reshape(-1, 1), mol2.mo.coeffs, atol=1e-8)
diff --git a/iodata/test/test_cp2k.py b/iodata/test/test_cp2k.py
index eed4a2e2..f11c850e 100644
--- a/iodata/test/test_cp2k.py
+++ b/iodata/test/test_cp2k.py
@@ -34,20 +34,17 @@
     from importlib.resources import path
 
 
-# TODO: add more obasis tests?
-
-
 def test_atom_si_uks():
     with path('iodata.test.data', 'atom_si.cp2k.out') as fn_out:
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [14])
     assert_equal(mol.atcorenums, [4])
-    assert mol.mo.type == 'unrestricted'
-    assert_equal(mol.mo.occs[:mol.mo.norba], [1, 2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0])
-    assert_equal(mol.mo.occs[mol.mo.norba:], [1, 0, 0, 0])
-    assert_allclose(mol.mo.energies[:mol.mo.norba],
+    assert mol.mo.kind == 'unrestricted'
+    assert_equal(mol.mo.occsa, [1, 2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0])
+    assert_equal(mol.mo.occsb, [1, 0, 0, 0])
+    assert_allclose(mol.mo.energiesa,
                     [-0.398761, -0.154896, -0.154896, -0.154896], atol=1.e-4)
-    assert_allclose(mol.mo.energies[mol.mo.norba:],
+    assert_allclose(mol.mo.energiesb,
                     [-0.334567, -0.092237, -0.092237, -0.092237], atol=1.e-4)
     assert_allclose(mol.energy, -3.761587698067, atol=1.e-10)
     assert len(mol.obasis.shells) == 3
@@ -58,8 +55,8 @@ def test_atom_si_uks():
     assert mol.obasis.shells[2].kinds == ['p']
     # check mo normalization
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp)
-    check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp)
+    check_orthonormal(mol.mo.coeffsa, olp)
+    check_orthonormal(mol.mo.coeffsb, olp)
 
 
 def test_atom_o_rks():
@@ -67,8 +64,8 @@ def test_atom_o_rks():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [8])
     assert_equal(mol.atcorenums, [6])
-    assert mol.mo.type == 'restricted'
-    assert_equal(mol.mo.occs, [1, 1, 1, 1])
+    assert mol.mo.kind == 'restricted'
+    assert_equal(mol.mo.occs, [2, 2, 2, 2])
     assert_allclose(mol.mo.energies, [0.102709, 0.606458, 0.606458, 0.606458], atol=1.e-4)
     assert_allclose(mol.energy, -15.464982778766, atol=1.e-10)
     assert_equal(mol.obasis.shells[0].angmoms, [0, 0])
@@ -88,19 +85,16 @@ def test_carbon_gs_ae_contracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [6])
-    assert mol.mo.type == 'unrestricted'
-    assert_allclose(mol.mo.occs[:mol.mo.norba],
-                    [1, 1, 2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0])
-    assert_allclose(mol.mo.energies[:mol.mo.norba],
-                    [-10.058194, -0.526244, -0.214978, -0.214978, -0.214978])
-    assert_allclose(mol.mo.occs[mol.mo.norba:], [1, 1, 0, 0, 0])
-    assert_allclose(mol.mo.energies[mol.mo.norba:],
-                    [-10.029898, -0.434300, -0.133323, -0.133323, -0.133323])
+    assert mol.mo.kind == 'unrestricted'
+    assert_allclose(mol.mo.occsa, [1, 1, 2 / 3, 2 / 3, 2 / 3])
+    assert_allclose(mol.mo.energiesa, [-10.058194, -0.526244, -0.214978, -0.214978, -0.214978])
+    assert_allclose(mol.mo.occsb, [1, 1, 0, 0, 0])
+    assert_allclose(mol.mo.energiesb, [-10.029898, -0.434300, -0.133323, -0.133323, -0.133323])
     assert_allclose(mol.energy, -37.836423363057)
     # check mo normalization
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp)
-    check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp)
+    check_orthonormal(mol.mo.coeffsa, olp)
+    check_orthonormal(mol.mo.coeffsb, olp)
 
 
 def test_carbon_gs_ae_uncontracted():
@@ -108,19 +102,16 @@ def test_carbon_gs_ae_uncontracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [6])
-    assert mol.mo.type == 'unrestricted'
-    assert_allclose(mol.mo.occs[:mol.mo.norba],
-                    [1, 1, 2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0])
-    assert_allclose(mol.mo.energies[:mol.mo.norba],
-                    [-10.050076, -0.528162, -0.217626, -0.217626, -0.217626])
-    assert_allclose(mol.mo.occs[mol.mo.norba:], [1, 1, 0, 0, 0])
-    assert_allclose(mol.mo.energies[mol.mo.norba:],
-                    [-10.022715, -0.436340, -0.137135, -0.137135, -0.137135])
+    assert mol.mo.kind == 'unrestricted'
+    assert_allclose(mol.mo.occsa, [1, 1, 2 / 3, 2 / 3, 2 / 3])
+    assert_allclose(mol.mo.energiesa, [-10.050076, -0.528162, -0.217626, -0.217626, -0.217626])
+    assert_allclose(mol.mo.occsb, [1, 1, 0, 0, 0])
+    assert_allclose(mol.mo.energiesb, [-10.022715, -0.436340, -0.137135, -0.137135, -0.137135])
     assert_allclose(mol.energy, -37.842552743398)
     # check mo normalization
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp)
-    check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp)
+    check_orthonormal(mol.mo.coeffsa, olp)
+    check_orthonormal(mol.mo.coeffsb, olp)
 
 
 def test_carbon_gs_pp_contracted():
@@ -128,18 +119,16 @@ def test_carbon_gs_pp_contracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [4])
-    assert mol.mo.type == 'unrestricted'
-    assert_allclose(mol.mo.occs[:mol.mo.norba], [1, 2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0])
-    assert_allclose(mol.mo.energies[:mol.mo.norba],
-                    [-0.528007, -0.219974, -0.219974, -0.219974])
-    assert_allclose(mol.mo.occs[mol.mo.norba:], [1, 0, 0, 0])
-    assert_allclose(mol.mo.energies[mol.mo.norba:],
-                    [-0.429657, -0.127060, -0.127060, -0.127060])
+    assert mol.mo.kind == 'unrestricted'
+    assert_allclose(mol.mo.occsa, [1, 2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0])
+    assert_allclose(mol.mo.energiesa, [-0.528007, -0.219974, -0.219974, -0.219974])
+    assert_allclose(mol.mo.occsb, [1, 0, 0, 0])
+    assert_allclose(mol.mo.energiesb, [-0.429657, -0.127060, -0.127060, -0.127060])
     assert_allclose(mol.energy, -5.399938535844)
     # check mo normalization
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp)
-    check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp)
+    check_orthonormal(mol.mo.coeffsa, olp)
+    check_orthonormal(mol.mo.coeffsb, olp)
 
 
 def test_carbon_gs_pp_uncontracted():
@@ -147,18 +136,16 @@ def test_carbon_gs_pp_uncontracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [4])
-    assert mol.mo.type == 'unrestricted'
-    assert_allclose(mol.mo.occs[:mol.mo.norba], [1, 2.0 / 3.0, 2.0 / 3.0, 2.0 / 3.0])
-    assert_allclose(mol.mo.energies[:mol.mo.norba],
-                    [-0.528146, -0.219803, -0.219803, -0.219803])
-    assert_allclose(mol.mo.occs[mol.mo.norba:], [1, 0, 0, 0])
-    assert_allclose(mol.mo.energies[mol.mo.norba:],
-                    [-0.429358, -0.126411, -0.126411, -0.126411])
+    assert mol.mo.kind == 'unrestricted'
+    assert_allclose(mol.mo.occsa, [1, 2 / 3, 2 / 3, 2 / 3])
+    assert_allclose(mol.mo.energiesa, [-0.528146, -0.219803, -0.219803, -0.219803])
+    assert_allclose(mol.mo.occsb, [1, 0, 0, 0])
+    assert_allclose(mol.mo.energiesb, [-0.429358, -0.126411, -0.126411, -0.126411])
     assert_allclose(mol.energy, -5.402288849332)
     # check mo normalization
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp)
-    check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp)
+    check_orthonormal(mol.mo.coeffsa, olp)
+    check_orthonormal(mol.mo.coeffsb, olp)
 
 
 def test_carbon_sc_ae_contracted():
@@ -166,8 +153,8 @@ def test_carbon_sc_ae_contracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [6])
-    assert mol.mo.type == 'restricted'
-    assert_allclose(mol.mo.occs, [1, 1, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0])
+    assert mol.mo.kind == 'restricted'
+    assert_allclose(mol.mo.occs, [2, 2, 2 / 3, 2 / 3, 2 / 3])
     assert_allclose(mol.mo.energies, [-10.067251, -0.495823, -0.187878, -0.187878, -0.187878])
     assert_allclose(mol.energy, -37.793939631890)
     # check mo normalization
@@ -180,8 +167,8 @@ def test_carbon_sc_ae_uncontracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [6])
-    assert mol.mo.type == 'restricted'
-    assert_allclose(mol.mo.occs, [1, 1, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0])
+    assert mol.mo.kind == 'restricted'
+    assert_allclose(mol.mo.occs, [2, 2, 2 / 3, 2 / 3, 2 / 3])
     assert_allclose(mol.mo.energies, [-10.062206, -0.499716, -0.192580, -0.192580, -0.192580])
     assert_allclose(mol.energy, -37.800453482378)
     # check mo normalization
@@ -194,8 +181,8 @@ def test_carbon_sc_pp_contracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [4])
-    assert mol.mo.type == 'restricted'
-    assert_allclose(mol.mo.occs, [1, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0])
+    assert mol.mo.kind == 'restricted'
+    assert_allclose(mol.mo.occs, [2, 2 / 3, 2 / 3, 2 / 3])
     assert_allclose(mol.mo.energies, [-0.500732, -0.193138, -0.193138, -0.193138])
     assert_allclose(mol.energy, -5.350765755382)
     # check mo normalization
@@ -208,8 +195,8 @@ def test_carbon_sc_pp_uncontracted():
         mol = load_one(str(fn_out))
     assert_equal(mol.atnums, [6])
     assert_equal(mol.atcorenums, [4])
-    assert mol.mo.type == 'restricted'
-    assert_allclose(mol.mo.occs, [1, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0])
+    assert mol.mo.kind == 'restricted'
+    assert_allclose(mol.mo.occs, [2, 2 / 3, 2 / 3, 2 / 3])
     assert_allclose(mol.mo.energies, [-0.500238, -0.192365, -0.192365, -0.192365])
     assert_allclose(mol.energy, -5.352864672201)
     # check mo normalization
diff --git a/iodata/test/test_fchk.py b/iodata/test/test_fchk.py
index bbeb3c3e..3201d65a 100644
--- a/iodata/test/test_fchk.py
+++ b/iodata/test/test_fchk.py
@@ -61,7 +61,7 @@ def test_load_fchk_hf_sto3g_num():
     assert mol.run_type == 'energy'
     assert mol.lot == 'rhf'
     assert mol.obasis_name == 'sto-3g'
-    assert mol.mo.type == 'restricted'
+    assert mol.mo.kind == 'restricted'
     assert mol.spinpol == 0
     assert mol.obasis.nbasis == 6
     assert len(mol.obasis.shells) == 3
@@ -157,46 +157,37 @@ def test_load_fchk_water_sto3g_hf():
 
 
 def test_load_fchk_lih_321g_hf():
-    fields = load_fchk_helper_internal('li_h_3-21G_hf_g09.fchk')
-    assert len(fields['obasis'].shells) == 5
-    assert fields['obasis'].nbasis == 11
-    assert len(fields['atcoords']) == len(fields['atnums'])
-    assert fields['atcoords'].shape[1] == 3
-    assert len(fields['atnums']) == 2
-
-    orb_alpha_coeffs = fields['mo'].coeffs[:, :fields['mo'].norba]
-    orb_alpha_energies = fields['mo'].energies[:fields['mo'].norba]
-    orb_alpha_occs = fields['mo'].occs[:fields['mo'].norba]
-
-    assert_allclose(orb_alpha_energies[0], (-2.76117), atol=1.e-4)
-    assert_allclose(orb_alpha_energies[-1], 0.97089, atol=1.e-4)
-    assert_allclose(orb_alpha_coeffs[0, 0], 0.99105, atol=1.e-4)
-    assert_allclose(orb_alpha_coeffs[1, 0], 0.06311, atol=1.e-4)
-    assert orb_alpha_coeffs[3, 2] < 1.e-4
-    assert_allclose(orb_alpha_coeffs[-1, 9], 0.13666, atol=1.e-4)
-    assert_allclose(orb_alpha_coeffs[4, -1], 0.17828, atol=1.e-4)
-    assert_equal(orb_alpha_occs.sum(), 2)
-    assert_equal(orb_alpha_occs.min(), 0.0)
-    assert_equal(orb_alpha_occs.max(), 1.0)
-
-    orb_beta_coeffs = fields['mo'].coeffs[:, fields['mo'].norba:]
-    orb_beta_energies = fields['mo'].energies[fields['mo'].norba:]
-    orb_beta_occs = fields['mo'].occs[fields['mo'].norba:]
-
-    assert_allclose(orb_beta_energies[0], -2.76031, atol=1.e-4)
-    assert_allclose(orb_beta_energies[-1], 1.13197, atol=1.e-4)
-    assert_allclose(orb_beta_coeffs[0, 0], 0.99108, atol=1.e-4)
-    assert_allclose(orb_beta_coeffs[1, 0], 0.06295, atol=1.e-4)
-    assert abs(orb_beta_coeffs[3, 2]) < 1e-4
-    assert_allclose(orb_beta_coeffs[-1, 9], 0.80875, atol=1.e-4)
-    assert_allclose(orb_beta_coeffs[4, -1], -0.15503, atol=1.e-4)
-    assert_equal(orb_beta_occs.sum(), 1)
-    assert_equal(orb_beta_occs.min(), 0.0)
-    assert_equal(orb_beta_occs.max(), 1.0)
-    assert_equal(orb_alpha_occs.shape[0], orb_alpha_coeffs.shape[0])
-    assert_equal(orb_beta_occs.shape[0], orb_beta_coeffs.shape[0])
-    energy = fields['energy']
-    assert_allclose(energy, -7.687331212191968E+00)
+    mol = load_fchk_helper('li_h_3-21G_hf_g09.fchk')
+    assert len(mol.obasis.shells) == 5
+    assert mol.obasis.nbasis == 11
+    assert len(mol.atcoords) == len(mol.atnums)
+    assert mol.atcoords.shape[1] == 3
+    assert len(mol.atnums) == 2
+    assert_allclose(mol.energy, -7.687331212191968E+00)
+
+    assert_allclose(mol.mo.energiesa[0], (-2.76117), atol=1.e-4)
+    assert_allclose(mol.mo.energiesa[-1], 0.97089, atol=1.e-4)
+    assert_allclose(mol.mo.coeffsa[0, 0], 0.99105, atol=1.e-4)
+    assert_allclose(mol.mo.coeffsa[1, 0], 0.06311, atol=1.e-4)
+    assert mol.mo.coeffsa[3, 2] < 1.e-4
+    assert_allclose(mol.mo.coeffsa[-1, 9], 0.13666, atol=1.e-4)
+    assert_allclose(mol.mo.coeffsa[4, -1], 0.17828, atol=1.e-4)
+    assert_equal(mol.mo.occsa.sum(), 2)
+    assert_equal(mol.mo.occsa.min(), 0.0)
+    assert_equal(mol.mo.occsa.max(), 1.0)
+
+    assert_allclose(mol.mo.energiesb[0], -2.76031, atol=1.e-4)
+    assert_allclose(mol.mo.energiesb[-1], 1.13197, atol=1.e-4)
+    assert_allclose(mol.mo.coeffsb[0, 0], 0.99108, atol=1.e-4)
+    assert_allclose(mol.mo.coeffsb[1, 0], 0.06295, atol=1.e-4)
+    assert abs(mol.mo.coeffsb[3, 2]) < 1e-4
+    assert_allclose(mol.mo.coeffsb[-1, 9], 0.80875, atol=1.e-4)
+    assert_allclose(mol.mo.coeffsb[4, -1], -0.15503, atol=1.e-4)
+    assert_equal(mol.mo.occsb.sum(), 1)
+    assert_equal(mol.mo.occsb.min(), 0.0)
+    assert_equal(mol.mo.occsb.max(), 1.0)
+    assert_equal(mol.mo.occsa.shape[0], mol.mo.coeffsa.shape[0])
+    assert_equal(mol.mo.occsb.shape[0], mol.mo.coeffsb.shape[0])
 
 
 def test_load_fchk_ghost_atoms():
@@ -456,10 +447,10 @@ def test_atmasses():
 
 def test_spinpol():
     mol1 = load_fchk_helper('ch3_rohf_sto3g_g03.fchk')
-    assert mol1.mo.type == 'restricted'
+    assert mol1.mo.kind == 'restricted'
     assert mol1.spinpol == 1
     mol2 = load_fchk_helper('li_h_3-21G_hf_g09.fchk')
-    assert mol2.mo.type == 'unrestricted'
+    assert mol2.mo.kind == 'unrestricted'
     assert mol2.spinpol == 1
     with pytest.raises(TypeError):
         mol2.spinpol = 2
diff --git a/iodata/test/test_molden.py b/iodata/test/test_molden.py
index 97c76dd4..9fa1439d 100644
--- a/iodata/test/test_molden.py
+++ b/iodata/test/test_molden.py
@@ -46,12 +46,14 @@ def test_load_molden_li2_orca():
 
     # Check geometry
     assert_equal(mol.atnums, [3, 3])
+    assert_allclose(mol.mo.occsa[:4], [1, 1, 1, 0])
+    assert_allclose(mol.mo.occsb[:4], [1, 1, 0, 0])
     assert_allclose(mol.atcoords[1], [5.2912331750, 0.0, 0.0])
 
     # Check normalization
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp, 1e-5)
-    check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp, 1e-5)
+    check_orthonormal(mol.mo.coeffsa, olp, 1e-5)
+    check_orthonormal(mol.mo.coeffsb, olp, 1e-5)
 
     # Check Mulliken charges
     charges = compute_mulliken_charges(mol)
@@ -68,6 +70,7 @@ def test_load_molden_h2o_orca():
 
     # Check geometry
     assert_equal(mol.atnums, [8, 1, 1])
+    assert_allclose(mol.mo.occs[:6], [2, 2, 2, 2, 2, 0])
     assert_allclose(mol.atcoords[2], [0.0, -0.1808833432, 1.9123825806])
 
     # Check normalization
diff --git a/iodata/test/test_molekel.py b/iodata/test/test_molekel.py
index c52ff440..f9e22cfe 100644
--- a/iodata/test/test_molekel.py
+++ b/iodata/test/test_molekel.py
@@ -72,8 +72,8 @@ def test_load_mkl_li2():
         mol = load_one(str(fn_mkl))
     # check mo normalization
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp, 1e-5)
-    check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp, 1e-5)
+    check_orthonormal(mol.mo.coeffsa, olp, 1e-5)
+    check_orthonormal(mol.mo.coeffsb, olp, 1e-5)
 
 
 def test_load_mkl_h2():
diff --git a/iodata/test/test_wfn.py b/iodata/test/test_wfn.py
index 4548d1f3..d6d817c4 100644
--- a/iodata/test/test_wfn.py
+++ b/iodata/test/test_wfn.py
@@ -134,11 +134,9 @@ def check_wfn(fn_wfn, nbasis, energy, charges_mulliken):
     assert mol.obasis.nbasis == nbasis
     # check orthonormal mo
     olp = compute_overlap(mol.obasis, mol.atcoords)
-    if mol.mo.type == 'restricted':
-        check_orthonormal(mol.mo.coeffs, olp, 1.e-5)
-    elif mol.mo.type == 'unrestricted':
-        check_orthonormal(mol.mo.coeffs[:, :mol.mo.norba], olp, 1.e-5)
-        check_orthonormal(mol.mo.coeffs[:, mol.mo.norba:], olp, 1.e-5)
+    check_orthonormal(mol.mo.coeffsa, olp, 1.e-5)
+    if mol.mo.kind == 'unrestricted':
+        check_orthonormal(mol.mo.coeffsb, olp, 1.e-5)
     # check energy & atomic charges
     if energy is not None:
         assert_allclose(mol.energy, energy, rtol=0., atol=1.e-5)
@@ -202,15 +200,15 @@ def test_load_wfn_o2_virtual():
                     -149.664140769678, np.array([0.0, 0.0]))
     # check MO occupation
     assert_equal(mol.mo.occs.shape, (88,))
-    assert_allclose(mol.mo.occs[:mol.mo.norba], [1.] * 9 + [0.] * 35)
-    assert_allclose(mol.mo.occs[mol.mo.norba:], [1.] * 7 + [0.] * 37)
+    assert_allclose(mol.mo.occsa, [1.] * 9 + [0.] * 35)
+    assert_allclose(mol.mo.occsb, [1.] * 7 + [0.] * 37)
     # check MO energies
     assert_equal(mol.mo.energies.shape, (88,))
-    mo_energies_a = mol.mo.energies[:mol.mo.norba]
+    mo_energies_a = mol.mo.energiesa
     assert_allclose(mo_energies_a[0], -20.752000, rtol=0, atol=1.e-6)
     assert_allclose(mo_energies_a[10], 0.179578, rtol=0, atol=1.e-6)
     assert_allclose(mo_energies_a[-1], 51.503193, rtol=0, atol=1.e-6)
-    mo_energies_b = mol.mo.energies[mol.mo.norba:]
+    mo_energies_b = mol.mo.energiesb
     assert_allclose(mo_energies_b[0], -20.697027, rtol=0, atol=1.e-6)
     assert_allclose(mo_energies_b[15], 0.322590, rtol=0, atol=1.e-6)
     assert_allclose(mo_energies_b[-1], 51.535258, rtol=0, atol=1.e-6)
@@ -239,4 +237,4 @@ def test_load_wfn_lih_cation_fci():
     assert_equal(mol.atnums, [3, 1])
     assert_equal(mol.mo.occs.shape, (11,))
     assert_allclose(mol.mo.occs.sum(), 3., rtol=0., atol=1.e-6)
-    # assert abs(mol.mo.occs[:mol.mo.norba].sum() - 1.5) < 1.e-6
+    # assert abs(mol.mo.occsa.sum() - 1.5) < 1.e-6

From b2766acfbfaf0a496d10b070c472dfe66dcb94d3 Mon Sep 17 00:00:00 2001
From: Toon Verstraelen <toon.verstraelen@UGent.be>
Date: Sun, 12 May 2019 22:05:19 +0200
Subject: [PATCH 3/3] Load irreps from molden file

---
 iodata/formats/molden.py   | 56 ++++++++++++++++++++++----------------
 iodata/test/common.py      |  3 +-
 iodata/test/test_molden.py | 13 +++++++++
 3 files changed, 47 insertions(+), 25 deletions(-)

diff --git a/iodata/formats/molden.py b/iodata/formats/molden.py
index a94b0844..9e934002 100644
--- a/iodata/formats/molden.py
+++ b/iodata/formats/molden.py
@@ -154,8 +154,8 @@ def _load_low(lit: LineIterator) -> dict:
         # molecular-orbital coefficients.
         elif line == '[mo]':
             data_alpha, data_beta = _load_helper_coeffs(lit)
-            coeffsa, energiesa, occsa = data_alpha
-            coeffsb, energiesb, occsb = data_beta
+            occsa, coeffsa, energiesa, irrepsa = data_alpha
+            occsb, coeffsb, energiesb, irrepsb = data_beta
 
     # Assign pure and Cartesian correctly. This needs to be done after reading
     # because the tags for pure functions may come after the basis set.
@@ -169,7 +169,7 @@ def _load_low(lit: LineIterator) -> dict:
             lit.error("Number of alpha orbital coefficients does not match the size of the basis.")
         mo = MolecularOrbitals(
             'restricted', coeffsa.shape[1], coeffsa.shape[1],
-            occsa, coeffsa, energiesa, None)
+            occsa, coeffsa, energiesa, irrepsa)
     else:
         if coeffsb.shape[0] != obasis.nbasis:
             lit.error("Number of beta orbital coefficients does not match the size of the basis.")
@@ -178,7 +178,7 @@ def _load_low(lit: LineIterator) -> dict:
             np.concatenate((occsa, occsb), axis=0),
             np.concatenate((coeffsa, coeffsb), axis=1),
             np.concatenate((energiesa, energiesb), axis=0),
-            None)
+            irrepsa + irrepsb)
 
     result = {
         'atcoords': atcoords,
@@ -249,16 +249,18 @@ def _load_helper_obasis(lit: LineIterator) -> MolecularBasis:
 
 def _load_helper_coeffs(lit: LineIterator) -> Tuple:
     """Load the orbital coefficients."""
+    occsa = []
     coeffsa = []
     energiesa = []
-    occsa = []
+    irrepsa = []
+    occsb = []
     coeffsb = []
     energiesb = []
-    occsb = []
+    irrepsb = []
 
     while True:
         try:
-            line = next(lit).lower().strip()
+            line = next(lit).strip()
         except StopIteration:
             # We have no proper way to check when a Molden file has ended, so
             # we must anticipate for it here.
@@ -281,18 +283,21 @@ def _load_helper_coeffs(lit: LineIterator) -> Tuple:
                 break
             key, value = line.split('=')
             info[key.strip().lower()] = value
-        energy = float(info['ene'])
         occ = float(info['occup'])
         col = []
+        energy = float(info['ene'])
+        irrep = info.get('sym', '??').strip()
         # store column of coefficients, i.e. one orbital, energy and occ
         if info['spin'].strip().lower() == 'alpha':
+            occsa.append(occ)
             coeffsa.append(col)
             energiesa.append(energy)
-            occsa.append(occ)
+            irrepsa.append(irrep)
         else:
+            occsb.append(occ)
             coeffsb.append(col)
             energiesb.append(energy)
-            occsb.append(occ)
+            irrepsb.append(irrep)
         for line in lit:
             words = line.split()
             if len(words) != 2 or not words[0].isdigit():
@@ -313,7 +318,7 @@ def _load_helper_coeffs(lit: LineIterator) -> Tuple:
         coeffsb = np.array(coeffsb).T
         energiesb = np.array(energiesb)
         occsb = np.array(occsb)
-    return (coeffsa, energiesa, occsa), (coeffsb, energiesb, occsb)
+    return (occsa, coeffsa, energiesa, irrepsa), (occsb, coeffsb, energiesb, irrepsb)
 
 
 def _is_normalized_properly(obasis: MolecularBasis, atcoords: np.ndarray,
@@ -664,27 +669,30 @@ def dump_one(f: TextIO, data: IOData):
     # Print the mean-field orbitals
     if data.mo.kind == 'unrestricted':
         f.write('[MO]\n')
-        _dump_helper_orb(f, 'Alpha', data.mo.energiesa, data.mo.occsa,
-                         data.mo.coeffsa[permutation] * signs.reshape(-1, 1))
-        _dump_helper_orb(f, 'Beta', data.mo.energiesb, data.mo.occsb,
-                         data.mo.coeffsb[permutation] * signs.reshape(-1, 1))
+        _dump_helper_orb(f, 'Alpha', data.mo.occsa,
+                         data.mo.coeffsa[permutation] * signs.reshape(-1, 1),
+                         data.mo.energiesa, data.mo.irrepsa)
+        _dump_helper_orb(f, 'Beta', data.mo.occsb,
+                         data.mo.coeffsb[permutation] * signs.reshape(-1, 1),
+                         data.mo.energiesb, data.mo.irrepsb)
     elif data.mo.kind == 'restricted':
         f.write('[MO]\n')
-        _dump_helper_orb(f, 'Alpha', data.mo.energies, data.mo.occs,
-                         data.mo.coeffs[permutation] * signs.reshape(-1, 1))
+        _dump_helper_orb(f, 'Alpha', data.mo.occs,
+                         data.mo.coeffs[permutation] * signs.reshape(-1, 1),
+                         data.mo.energies, data.mo.irreps)
     else:
         raise NotImplementedError
 
 
-def _dump_helper_orb(f, spin, orb_energies, orb_occs, orb_coeffs):
-    for ifn in range(orb_coeffs.shape[1]):
-        f.write(f' Ene= {orb_energies[ifn]:.17e}\n')
-        f.write(' Sym=     1a\n')
+def _dump_helper_orb(f, spin, occs, coeffs, energies, irreps):
+    for ifn in range(coeffs.shape[1]):
+        f.write(f' Ene= {energies[ifn]:.17e}\n')
+        f.write(f' Sym= {irreps[ifn]}\n')
         f.write(f' Spin= {spin}\n')
-        f.write(f' Occup= {orb_occs[ifn]:.17e}\n')
-        for ibasis in range(orb_coeffs.shape[0]):
+        f.write(f' Occup= {occs[ifn]:.17e}\n')
+        for ibasis in range(coeffs.shape[0]):
             # The original molden floating-point formatting is too low
             # precision. Molden also reads high-precision, so we use this
             # instead.
             # f.write('{:4d} {:10.6f}\n'.format(ibasis + 1, orb_coeffs[ibasis, ifn]))
-            f.write('{:4d} {:.17e}\n'.format(ibasis + 1, orb_coeffs[ibasis, ifn]))
+            f.write('{:4d} {:.17e}\n'.format(ibasis + 1, coeffs[ibasis, ifn]))
diff --git a/iodata/test/common.py b/iodata/test/common.py
index 5e7f05e8..c1fcee7a 100644
--- a/iodata/test/common.py
+++ b/iodata/test/common.py
@@ -111,8 +111,9 @@ def compare_mols(mol1, mol2):
     permutation, signs = convert_conventions(mol1.obasis, mol2.obasis.conventions)
     assert mol1.mo.kind == mol2.mo.kind
     assert_allclose(mol1.mo.occs, mol2.mo.occs)
-    assert_allclose(mol1.mo.energies, mol2.mo.energies)
     assert_allclose(mol1.mo.coeffs[permutation] * signs.reshape(-1, 1), mol2.mo.coeffs, atol=1e-8)
+    assert_allclose(mol1.mo.energies, mol2.mo.energies)
+    assert_equal(mol1.mo.irreps, mol2.mo.irreps)
     # operators and density matrices
     cases = [
         ('one_ints', ['olp', 'kin_ao', 'na_ao']),
diff --git a/iodata/test/test_molden.py b/iodata/test/test_molden.py
index 9fa1439d..adc6ecca 100644
--- a/iodata/test/test_molden.py
+++ b/iodata/test/test_molden.py
@@ -48,6 +48,9 @@ def test_load_molden_li2_orca():
     assert_equal(mol.atnums, [3, 3])
     assert_allclose(mol.mo.occsa[:4], [1, 1, 1, 0])
     assert_allclose(mol.mo.occsb[:4], [1, 1, 0, 0])
+    assert_equal(mol.mo.irreps, ['1a'] * mol.mo.norb)
+    assert_equal(mol.mo.irrepsa, ['1a'] * mol.mo.norba)
+    assert_equal(mol.mo.irrepsb, ['1a'] * mol.mo.norbb)
     assert_allclose(mol.atcoords[1], [5.2912331750, 0.0, 0.0])
 
     # Check normalization
@@ -71,6 +74,7 @@ def test_load_molden_h2o_orca():
     # Check geometry
     assert_equal(mol.atnums, [8, 1, 1])
     assert_allclose(mol.mo.occs[:6], [2, 2, 2, 2, 2, 0])
+    assert_equal(mol.mo.irreps, ['1a'] * mol.mo.norb)
     assert_allclose(mol.atcoords[2], [0.0, -0.1808833432, 1.9123825806])
 
     # Check normalization
@@ -255,6 +259,15 @@ def test_load_molden_nh3_turbomole():
     assert_allclose(charges, molden_charges, atol=1.e-3)
 
 
+def test_load_molden_f():
+    with path('iodata.test.data', 'F.molden') as fn_molden:
+        mol = load_one(str(fn_molden))
+    assert_allclose(mol.mo.occsa[:6], [1, 1, 1, 1, 1, 0])
+    assert_allclose(mol.mo.occsb[:6], [1, 1, 1, 1, 0, 0])
+    assert_equal(mol.mo.irrepsa[:6], ['Ag', 'Ag', 'B3u', 'B2u', 'B1u', 'B3u'])
+    assert_equal(mol.mo.irrepsb[:6], ['Ag', 'Ag', 'B3u', 'B2u', 'B1u', 'B3u'])
+
+
 def check_load_dump_consistency(fn, tmpdir):
     """Check if data is preserved after dumping and loading a Molden file.