Merge branch 'master' into attribute_validation

iodata/basis.py

@@ -87,24 +87,25 @@ def angmom_its(angmom: Union[int, List[int]]) -> Union[str, List[str]]:
 
 @attr.s(auto_attribs=True, slots=True)
 class Shell:
-    """Describe a single shell in a molecular basis set.
+    """A shell in a molecular basis representing (generalized) contractions with the same exponents.
 
     Attributes
     ----------
     icenter
-        An integer referring to a row in the array atcoords in an IOData object.
+        An integer index specifying the row in the atcoords array of IOData object.
     angmoms
         An integer array of angular momentum quantum numbers, non-negative, with
         shape (ncon,).
     kinds
-        List of strings describing the kind of contraction: 'c' for Cartesian
+        List of strings describing the kind of contractions: 'c' for Cartesian
         and 'p' for pure. Pure functions are only allowed for angmom>1.
         The length equals the number of contractions: len(angmoms)=ncon.
     exponents
-        an array of exponents of primitives, with shape (nprim,).
+        The array containing the exponents of the primitives, with shape (nprim,).
     coeffs
-        an array with contraction coefficients, with shape (nprim, ncon). These
-        coefficients assume that the primitives are L2 (orbitals) or L1
+        The array containing the coefficients of the normalized primitives in each contraction;
+        shape = (nprim, ncon).
+        These coefficients assume that the primitives are L2 (orbitals) or L1
         (densities) normalized, but contractions are not necessarily normalized.
         (This depends on the code which generated the contractions.)
 
@@ -117,26 +118,26 @@ class Shell:
     coeffs: np.ndarray = attr.ib(validator=validate_shape(("exponents", 0), ("kinds", 0)))
 
     @property
-    def nbasis(self) -> int:
-        """Return the number of basis functions."""
+    def nbasis(self) -> int:   # noqa: D401
+        """Number of basis functions (e.g. 3 for a P shell and 4 for an SP shell)."""
         result = 0
         for angmom, kind in zip(self.angmoms, self.kinds):
             if kind == 'c':  # Cartesian
                 result += ((angmom + 1) * (angmom + 2)) // 2
             elif kind == 'p' and angmom >= 2:
                 result += 2 * angmom + 1
             else:
-                raise TypeError('Unknown shell kind \'{}\'.'.format(kind))
+                raise TypeError('Unknown shell kind \'{}\'; expected \'c\' or \'p\'.'.format(kind))
         return result
 
     @property
-    def nprim(self) -> int:
-        """Return the number of primitives. Also known as the contraction length."""
+    def nprim(self) -> int:   # noqa: D401
+        """Number of primitives, also known as the contraction length."""
         return len(self.exponents)
 
     @property
-    def ncon(self) -> int:
-        """Return the number of contractions."""
+    def ncon(self) -> int:   # noqa: D401
+        """Number of contractions. This is usually 1; e.g., it would be 2 for an SP shell."""
         return len(self.angmoms)
 
     def evolve(self, **changes):
@@ -145,15 +146,17 @@ def evolve(self, **changes):
 
 
 class MolecularBasis(NamedTuple):
-    """Describe a complete molecular orbital or density basis set.
+    """A complete molecular orbital or density basis set.
 
     Attributes
     ----------
     shells
-        a list of objects of the type Shell
+        A list of objects of type Shell which can support generalized contractions.
     conventions
-        a dictionary with as key a typle of angular momentum integer and kind
-        character, and as value a list of basis function strings, e.g.
+        A dictionary specifying the ordered basis functions for a given angular momentum and kind.
+        The key is a tuple of angular momentum integer and kind character ('c' for Cartesian
+        and 'p' for pure/spherical) and the value is a list of basis function strings.
+        For example,
 
         .. code-block:: python
 
@@ -175,7 +178,7 @@ class MolecularBasis(NamedTuple):
                 (2, 'p'): ['s2', 's1', 'c0', 'c1', 'c2'],
                 # Different quantum-chemistry codes may use incompatible
                 # orderings and sign conventions. E.g. Molden files written
-                # by Orca use the following convention for pure f functions:
+                # by ORCA use the following convention for pure f functions:
                 (3, 'p'): ['c0', 'c1', 's1', 'c2', 's2', '-c3', '-s3'],
                 # Note that the minus sign in the last two basis functions
                 # denotes that the signs of these harmonics have been changed.
@@ -185,7 +188,8 @@ class MolecularBasis(NamedTuple):
         in :ref:`basis_conventions`.
 
     primitive_normalization
-        Either 'L1' or 'L2'.
+        The normalization convention of primitives, which can be 'L2' (orbitals) or 'L1'
+        (densities) normalized.
 
     """
 
@@ -194,8 +198,8 @@ class MolecularBasis(NamedTuple):
     primitive_normalization: str
 
     @property
-    def nbasis(self) -> int:
-        """Return the number of basis functions."""
+    def nbasis(self) -> int:   # noqa: D401
+        """Number of basis functions."""
         return sum(shell.nbasis for shell in self.shells)
 
     def get_segmented(self):

iodata/formats/molden.py

@@ -19,7 +19,7 @@
 """Molden file format.
 
 Many QC codes can write out Molden files, e.g. `Molpro <https://www.molpro.net/>`_,
-`Orca <https://orcaforum.cec.mpg.de/>`_, `PSI4 <http://www.psicode.org/>`_,
+`Orca <https://sites.google.com/site/orcainputlibrary/>`_, `PSI4 <http://www.psicode.org/>`_,
 `Molden <http://www.cmbi.ru.nl/molden/>`_, `Turbomole <http://www.turbomole.com/>`_. Keep
 in mind that several of these write incorrect versions of the file format, but these
 errors are corrected when loading them with IOData.

iodata/formats/molekel.py

@@ -20,17 +20,19 @@
 
 This format is used by two programs:
 `Molekel <http://ugovaretto.github.io/molekel/wiki/pmwiki.php/Main/HomePage.html>`_ and
-`Orca <https://orcaforum.cec.mpg.de/>`_.
+`Orca <https://sites.google.com/site/orcainputlibrary/>`_.
 """
 
+import warnings
 
-from typing import Tuple, List
+from typing import Tuple, List, TextIO
 
 import numpy as np
 
 from .molden import CONVENTIONS, _fix_molden_from_buggy_codes
-from ..basis import angmom_sti, MolecularBasis, Shell
-from ..docstrings import document_load_one
+from ..basis import angmom_sti, angmom_its, convert_conventions, MolecularBasis, Shell
+from ..docstrings import document_load_one, document_dump_one
+from ..iodata import IOData
 from ..orbitals import MolecularOrbitals
 from ..utils import angstrom, LineIterator
 
@@ -47,6 +49,17 @@ def _load_helper_charge_spinpol(lit: LineIterator) -> List[int]:
     return charge, spinpol
 
 
+def _load_helper_charges(lit: LineIterator) -> dict:
+    atcharges = []
+    for line in lit:
+        line = line.strip()
+        if line == '$END':
+            break
+        atcharges.append(float(line))
+
+    return {'mulliken': np.array(atcharges)}
+
+
 def _load_helper_atoms(lit: LineIterator) -> Tuple[np.ndarray, np.ndarray]:
     atnums = []
     atcoords = []
@@ -100,6 +113,7 @@ def _load_helper_obasis(lit: LineIterator) -> MolecularBasis:
 def _load_helper_coeffs(lit: LineIterator, nbasis: int) -> Tuple[np.ndarray, np.ndarray]:
     coeffs = []
     energies = []
+    irreps = []
 
     in_orb = 0
     for line in lit:
@@ -112,7 +126,7 @@ def _load_helper_coeffs(lit: LineIterator, nbasis: int) -> Tuple[np.ndarray, np.
             ncol = len(words)
             assert ncol > 0
             for word in words:
-                assert word == 'a1g'
+                irreps.append(word)
             cols = [np.zeros((nbasis, 1), float) for _ in range(ncol)]
             in_orb = 1
         elif in_orb == 1:
@@ -134,7 +148,7 @@ def _load_helper_coeffs(lit: LineIterator, nbasis: int) -> Tuple[np.ndarray, np.
                 in_orb = 0
                 coeffs.extend(cols)
 
-    return np.hstack(coeffs), np.array(energies)
+    return np.hstack(coeffs), np.array(energies), irreps
 
 
 def _load_helper_occ(lit: LineIterator) -> np.ndarray:
@@ -149,7 +163,7 @@ def _load_helper_occ(lit: LineIterator) -> np.ndarray:
 
 
 # pylint: disable=too-many-branches,too-many-statements
-@document_load_one("Molekel", ['atcoords', 'atnums', 'mo', 'obasis'])
+@document_load_one("Molekel", ['atcoords', 'atnums', 'mo', 'obasis'], ['atcharges'])
 def load_one(lit: LineIterator) -> dict:
     """Do not edit this docstring. It will be overwritten."""
     charge = None
@@ -162,6 +176,9 @@ def load_one(lit: LineIterator) -> dict:
     coeffsb = None
     energiesb = None
     occsb = None
+    atcharges = None
+    irrepsa = None
+    irrepsb = None
     # Using a loop because we're not entirely sure if sections in an MKL file
     # have a fixed order.
     while True:
@@ -173,16 +190,18 @@ def load_one(lit: LineIterator) -> dict:
             break
         if line == '$CHAR_MULT':
             charge, spinpol = _load_helper_charge_spinpol(lit)
+        elif line == '$CHARGES':
+            atcharges = _load_helper_charges(lit)
         elif line == '$COORD':
             atnums, atcoords = _load_helper_atoms(lit)
         elif line == '$BASIS':
             obasis = _load_helper_obasis(lit)
         elif line == '$COEFF_ALPHA':
-            coeffsa, energiesa = _load_helper_coeffs(lit, obasis.nbasis)
+            coeffsa, energiesa, irrepsa = _load_helper_coeffs(lit, obasis.nbasis)
         elif line == '$OCC_ALPHA':
             occsa = _load_helper_occ(lit)
         elif line == '$COEFF_BETA':
-            coeffsb, energiesb = _load_helper_coeffs(lit, obasis.nbasis)
+            coeffsb, energiesb, irrepsb = _load_helper_coeffs(lit, obasis.nbasis)
         elif line == '$OCC_BETA':
             occsb = _load_helper_occ(lit)
 
@@ -204,7 +223,7 @@ def load_one(lit: LineIterator) -> dict:
         assert abs(occsa.sum() - nelec) < 1e-7
         mo = MolecularOrbitals(
             'restricted', coeffsa.shape[1], coeffsa.shape[1],
-            occsa, coeffsa, energiesa, None)
+            occsa, coeffsa, energiesa, irrepsa)
     else:
         if occsb is None:
             lit.error('Beta occupation numbers not found in mkl file while '
@@ -223,13 +242,152 @@ def load_one(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,
         'atnums': atnums,
         'obasis': obasis,
         'mo': mo,
+        'atcharges': atcharges,
     }
     _fix_molden_from_buggy_codes(result, lit)
     return result
+
+
+@document_dump_one("Molekel", ['atcoords', 'atnums', 'mo', 'obasis'], ['atcharges'])
+def dump_one(f: TextIO, data: IOData):
+    """Do not edit this docstring. It will be overwritten."""
+    # Header
+    f.write('$MKL\n')
+    f.write('#\n')
+    f.write('# MKL format file produced by IOData\n')
+    f.write('#\n')
+
+    # CHAR_MUL
+    f.write('$CHAR_MULT\n')
+    f.write('  {:.0f} {:.0f}\n'.format(data.charge, data.spinpol + 1))
+    f.write('$END\n')
+    f.write('\n')
+
+    # COORD
+    atcoords = data.atcoords / angstrom
+    f.write('$COORD\n')
+    for n, coord in zip(data.atnums, atcoords):
+        f.write('   {:d}   {: ,.6f}  {: ,.6f}  {: ,.6f}\n'.format(n, coord[0], coord[1], coord[2]))
+    f.write('$END\n')
+    f.write('\n')
+
+    # CHARGES
+    if 'mulliken' in data.atcharges:
+        charges = data.atcharges['mulliken']
+    else:
+        charges = {}
+        warnings.warn('Skip writing Mulliken charges, because they are not stored in IOData.')
+    f.write('$CHARGES\n')
+    for charge in charges:
+        f.write('  {: ,.6f}\n'.format(charge))
+    f.write('$END\n')
+    f.write('\n')
+
+    # BASIS
+    f.write('$BASIS\n')
+    iatom_last = 0
+    for shell in data.obasis.shells:
+        iatom_new = shell.icenter
+        if iatom_new != iatom_last:
+            f.write('$$\n')
+        for iangmom, (angmom, kind) in enumerate(zip(shell.angmoms, shell.kinds)):
+            iatom_last = shell.icenter
+            nbasis = len(CONVENTIONS[(angmom, kind)])
+            f.write(' {} {:1s} 1.00\n'.format(nbasis, angmom_its(angmom).capitalize()))
+            for exponent, coeff in zip(shell.exponents, shell.coeffs[:, iangmom]):
+                f.write('{:20.10f} {:17.10f}\n'.format(exponent, coeff))
+    f.write('\n')
+    f.write('$END\n')
+    f.write('\n')
+
+    if data.mo.kind == 'restricted':
+        # COEFF_ALPHA
+        f.write('$COEFF_ALPHA\n')
+        _dump_helper_coeffs(f, data, spin='a')
+
+        # OCC_ALPHA
+        f.write('$OCC_ALPHA\n')
+        _dump_helper_occ(f, data, spin='ab')
+
+    # Not taking into account generalized.
+    elif data.mo.kind == 'unrestricted':
+        # COEFF_ALPHA
+        f.write('$COEFF_ALPHA\n')
+        _dump_helper_coeffs(f, data, spin='a')
+
+        # OCC_ALPHA
+        f.write('$OCC_ALPHA\n')
+        _dump_helper_occ(f, data, spin='a')
+        f.write('\n')
+
+        # COEFF_BETA
+        f.write('$COEFF_BETA\n')
+        _dump_helper_coeffs(f, data, spin='b')
+
+        # OCC_BETA
+        f.write('$OCC_BETA\n')
+        _dump_helper_occ(f, data, spin='b')
+
+    else:
+        raise ValueError(f"The MKL format does not support {data.mo.kind} orbitals.")
+
+
+# Defining help dumping functions
+def _dump_helper_coeffs(f, data, spin=None):
+    permutation, signs = convert_conventions(data.obasis, CONVENTIONS)
+    if spin == 'a':
+        norb = data.mo.norba
+        coeff = data.mo.coeffsa[permutation] * signs.reshape(-1, 1)
+        ener = data.mo.energiesa
+        if data.mo.irreps is not None:
+            irreps = data.mo.irreps[:norb]
+        else:
+            irreps = ['a1g'] * norb
+    elif spin == 'b':
+        norb = data.mo.norbb
+        coeff = data.mo.coeffsb[permutation] * signs.reshape(-1, 1)
+        ener = data.mo.energiesb
+        if data.mo.irreps is not None:
+            irreps = data.mo.irreps[norb:]
+        else:
+            irreps = ['a1g'] * norb
+    else:
+        raise IOError('A spin must be specified')
+
+    for j in range(0, norb, 5):
+        en = ' '.join(['   {: ,.12f}'.format(e) for e in ener[j:j + 5]])
+        irre = ' '.join(['{}'.format(irr) for irr in irreps[j:j + 5]])
+        f.write(irre + '\n')
+        f.write(en + '\n')
+        for orb in coeff[:, j:j + 5]:
+            coeffs = ' '.join(['  {: ,.12f}'.format(c) for c in orb])
+            f.write(coeffs + '\n')
+
+    f.write(' $END\n')
+    f.write('\n')
+
+
+def _dump_helper_occ(f, data, spin=None):
+    if spin == 'a':
+        norb = data.mo.norba
+        occ = data.mo.occsa
+    elif spin == 'b':
+        norb = data.mo.norbb
+        occ = data.mo.occsb
+    elif spin == 'ab':
+        norb = data.mo.norba
+        occ = data.mo.occs
+    else:
+        raise IOError('A spin must be specified')
+
+    for j in range(0, norb, 5):
+        occs = ' '.join(['  {: ,.7f}'.format(o) for o in occ[j:j + 5]])
+        f.write(occs + '\n')
+    f.write(' $END\n')