Fix "unsupported operand" issue when accessing IOData.charge

iodata/formats/qchemlog.py

@@ -37,19 +37,17 @@
 PATTERNS = ['*.qchemlog']
 
 
-@document_load_one("qchemlog", ['atcoords', 'atmasses', 'atnums', 'charge', 'energy', 'g_rot',
-                                'mo', 'lot', 'nelec', 'obasis_name', 'run_type', 'extra'],
+@document_load_one("qchemlog", ['atcoords', 'atmasses', 'atnums', 'energy', 'g_rot',
+                                'mo', 'lot', 'obasis_name', 'run_type', 'extra'],
                    ['athessian'])
 def load_one(lit: LineIterator) -> dict:
     """Do not edit this docstring. It will be overwritten."""
     data = load_qchemlog_low(lit)
 
     # add these labels if they are loaded
-    result_labels = ['atcoords', 'atmasses', 'atnums', 'charge', 'charge', 'energy', 'g_rot',
+    result_labels = ['atcoords', 'atmasses', 'atnums', 'energy', 'g_rot',
                      'run_type', 'athessian', 'lot', 'obasis_name']
     result = {label: data[label] for label in result_labels if data.get(label) is not None}
-    # add number of electrons
-    result['nelec'] = data['alpha_elec'] + data['beta_elec']
 
     # mulliken charges
     if data.get("mulliken_charges") is not None:
@@ -125,8 +123,7 @@ def load_qchemlog_low(lit: LineIterator) -> dict:
 
         # get the atomic information
         if line.startswith('$molecule'):
-            result = _helper_atoms(lit)
-            data['charge'], data['spin_multi'], data['atnums'], data['atcoords'] = result
+            data['atnums'], data['atcoords'] = _helper_atoms(lit)
             data['natom'] = len(data['atnums'])
         # job specifications
         elif line.startswith('$rem'):
@@ -168,8 +165,8 @@ def load_qchemlog_low(lit: LineIterator) -> dict:
 
 def _helper_atoms(lit: LineIterator) -> Tuple:
     """Load list of coordinates from an Q-Chem log output file format."""
-    # net charge and spin multiplicity
-    charge, spin_multi = [int(i) for i in next(lit).strip().split()]
+    # Skip line with net charge and spin multiplicity
+    next(lit)
     # atomic numbers and atomic coordinates (in Angstrom)
     atom_symbols = []
     atcoords = []
@@ -180,7 +177,7 @@ def _helper_atoms(lit: LineIterator) -> Tuple:
         atcoords.append([float(i) for i in line.strip().split()[1:]])
     atnums = np.array([sym2num[i] for i in atom_symbols])
     atcoords = np.array(atcoords)
-    return charge, spin_multi, atnums, atcoords
+    return atnums, atcoords
 
 
 def _helper_job(lit: LineIterator) -> Tuple:

iodata/iodata.py

@@ -114,7 +114,8 @@ class IOData:
         e.g. for a wire. Two vectors describe a 2D cell, e.g. for a membrane.
         Three vectors describe a 3D cell, e.g. a crystalline solid.
     charge
-        The net charge of the system.
+        The net charge of the system. When possible, this is derived from
+        atcorenums and nelec.
     core_energy
         The Hartree-Fock energy due to the core orbitals
     cube
@@ -247,36 +248,66 @@ class IOData:
     two_ints: dict = {}
     two_rdms: dict = {}
 
+    def __attrs_post_init__(self):
+        # Trigger setter to acchieve consistency in properties
+        # atcorenums, nelec, charge, spinmult. This is needed because the
+        # attr constructor bypasses the setters. See
+        # https://www.attrs.org/en/stable/init.html#private-attributes
+        if self._atcorenums is not None:
+            self.atcorenums = self._atcorenums
+        if self._charge is not None:
+            self.charge = self._charge
+        if self._nelec is not None:
+            self.nelec = self._nelec
+        if self._spinpol is not None:
+            self.spinpol = self._spinpol
+
     # Public interfaces to private attributes
 
     @property
     def atcorenums(self) -> np.ndarray:
         """Return effective core charges."""
-        result = self._atcorenums
-        if result is None and self.atnums is not None:
+        if self._atcorenums is None and self.atnums is not None:
             # Known bug in pylint. See
             # https://stackoverflow.com/questions/47972143/using-attr-with-pylint
             # https://github.com/PyCQA/pylint/issues/1694
             # pylint: disable=no-member
-            result = self.atnums.astype(float)
-            self._atcorenums = result
-        return result
+            self.atcorenums = self.atnums.astype(float)
+        return self._atcorenums
 
     @atcorenums.setter
     def atcorenums(self, atcorenums):
-        self._atcorenums = atcorenums
+        if atcorenums is None:
+            if self.nelec is not None and self._atcorenums is not None:
+                # Set _charge because charge can no longer be derived from
+                # atcorenums and nelec.
+                self._charge = self._atcorenums.sum() - self.nelec
+            self._atcorenums = None
+        else:
+            if self._charge is not None:
+                # _charge is treated as the dependent one, while atcorenums and
+                # nelec are treated as independent.
+                if self._nelec is None:
+                    # Switch to storing _nelec.
+                    self._nelec = atcorenums.sum() - self._charge
+                self._charge = None
+            self._atcorenums = np.asarray(atcorenums, dtype=float)
 
     @property
     def charge(self) -> float:
         """Return the net charge of the system."""
-        if self.atcorenums is None:
+        # The internal _charge is used only if it cannot be derived.
+        if self.atcorenums is None or self.nelec is None:
             return self._charge
         return self.atcorenums.sum() - self.nelec
 
     @charge.setter
     def charge(self, charge: float):
+        # The internal _charge is used only if atcorenums is None.
         if self.atcorenums is None:
             self._charge = charge
+        elif charge is None:
+            self.nelec = None
         else:
             self.nelec = self.atcorenums.sum() - charge
 
@@ -286,8 +317,8 @@ def natom(self) -> int:
         natom = None
         if self.atcoords is not None:
             natom = len(self.atcoords)
-        elif self.atcorenums is not None:
-            natom = len(self.atcorenums)
+        elif self._atcorenums is not None:
+            natom = len(self._atcorenums)
         elif self.atgradient is not None:
             natom = len(self.atgradient)
         elif self.atfrozen is not None:
@@ -301,9 +332,11 @@ def natom(self) -> int:
     @property
     def nelec(self) -> float:
         """Return the number of electrons."""
-        if self.mo is None:
-            return self._nelec
-        return self.mo.nelec
+        # When the molecular orbitals are present, they determine the number
+        # of electrons. Only when mo is absent, we use a stored value.
+        if self.mo is not None:
+            return self.mo.nelec
+        return self._nelec
 
     @nelec.setter
     def nelec(self, nelec: float):
@@ -320,12 +353,14 @@ def spinpol(self) -> float:
         number in ]0, 2[ implies spin polarizaiton, which may not always be a
         valid assumption.
         """
-        if self.mo is None:
-            return self._spinpol
-        return self.mo.spinpol
+        if self.mo is not None:
+            return self.mo.spinpol
+        return self._spinpol
 
     @spinpol.setter
     def spinpol(self, spinpol: float):
+        # When the molecular orbitals are present, they determine the spin
+        # polarization. Only when mo is absent, we use a stored value.
         if self.mo is None:
             self._spinpol = spinpol
         else: