Rename and restructure MolecularSystem class

examples/coulomb_sturmian/find_optimal_k.py

@@ -28,8 +28,11 @@
 n_max = 5
 
 for atom, mult in [("He", 1), ("Be", 1), ("C", 3), ("Ne", 1)]:
+    system = molsturm.System(atom)
+    system.multiplicity = mult
+
     scfparams = molsturm.ScfParameters()
-    scfparams.system = molsturm.MolecularSystem(atom, multiplicity=mult)
+    scfparams.system = system
 
     k_guess = molsturm.sturmian.cs.empirical_kopt(scfparams.system)
     scfparams.basis = molsturm.construct_basis("sturmian/atomic/cs_reference_pc",

examples/posthf/be_ccpvdz_mp2.py

@@ -24,7 +24,7 @@
 import molsturm
 import molsturm.posthf
 
-sys = molsturm.MolecularSystem([4], [[0, 0, 0]])
+sys = molsturm.System([4], [[0, 0, 0]])
 res = molsturm.hartree_fock(sys, basis_type="gaussian", basis_set_name="cc-pvdz",
                             print_iterations=True, conv_tol=1e-10)
 molsturm.print_convergence_summary(res)

examples/posthf/be_cs_14_mp2.py

@@ -24,7 +24,7 @@
 import molsturm
 import molsturm.posthf
 
-sys = molsturm.MolecularSystem(["beryllium"], [[0, 0, 0]])
+sys = molsturm.System(["beryllium"], [[0, 0, 0]])
 bas = molsturm.construct_basis("sturmian/atomic", sys, k_exp=2.1, n_max=11, l_max=0,
                                backend="cs_reference_pc")
 res = molsturm.hartree_fock(sys, bas, conv_tol=1e-10, print_iterations=True)

examples/single_point/c_cs_py

@@ -24,8 +24,11 @@
 import molsturm
 
 def main():
+    carbon = molsturm.System.("c")
+    carbon.multiplicity = 3
+
     params = molsturm.ScfParameters()
-    params.system = molsturm.MolecularSystem("C", multiplicity=3)
+    params.system = carbon
     params.basis = molsturm.construct_basis("sturmian/atomic/cs_reference_pc",
                                             params.system, k_exp=3.3, n_max=4,
                                             l_max=3)

examples/state_interface/coupled_cluster_doubles.py

@@ -116,11 +116,11 @@ def ccd(state):
 
 
 if __name__ == "__main__":
-    sys = molsturm.MolecularSystem(
+    sys = molsturm.System.by_charge_multiplicity(
         atoms=["O", "O"],
         coords=[(0, 0, 0), (0, 0, 2.8535)],
-        multiplicity=3
     )
+    sys.multiplicity = 3
     state = molsturm.hartree_fock(sys, basis_type="gaussian",
                                   basis_set_name="6-31g", conv_tol=5e-7)
 

src/interface/python/molsturm/ScfParameters.py

@@ -25,7 +25,7 @@
 import gint.element
 from gint import split_basis_type
 from .ParameterMap import ParameterMap
-from .MolecularSystem import MolecularSystem
+from .System import System
 import gint.gaussian
 import gint.sturmian.atomic
 from ._iface_conversion import ParamSpecial
@@ -129,9 +129,9 @@ def from_args(cls, system, basis, conv_tol=None, max_iter=None, n_eigenpairs=Non
         """
 
         if isinstance(system, str):
-            system = MolecularSystem(system)
-        elif not isinstance(system, MolecularSystem):
-            raise TypeError("The first argument needs to be a MolecularSystem object or "
+            system = System(system)
+        elif not isinstance(system, System):
+            raise TypeError("The first argument needs to be a System object or "
                             "a string of exactly one atom.")
 
         # Build a parameter tree from the commandline arguments
@@ -193,7 +193,7 @@ def __normalise_system(self):
 
         system = self["system"]
 
-        # TODO This could be made much simpler if the MolecularSystem
+        # TODO This could be made much simpler if the System
         #      class had a from_params function which constructed
         #      (and checked) the class from such a subtree.
 
@@ -214,32 +214,26 @@ def __normalise_system(self):
         self.__normalise_numpy_array("system/coords", (n_atom, 3), dtype=float)
         self.__normalise_numpy_array("system/atom_numbers", (n_atom,), dtype=int)
 
-        if "multiplicity" in system:
-            if "n_alpha" in system or "n_beta" in system:
-                warnings.warn("Overriding system/n_alpha and system/n_beta in "
-                              "ScfParameters, since system/multiplicity is present.")
+        if "multiplicity" in system or "charge" in system or \
+           ("n_alpha" not in system and "n_beta" not in system):
+            sysobj = System(atoms=system["atom_numbers"].value,
+                            coords=system["coords"].value)
+
             try:
-                sysobj = MolecularSystem(
-                    atoms=system["atom_numbers"].value,
-                    coords=system["coords"].value,
-                    multiplicity=system.get("multiplicity", None),
-                    charge=system.get("charge", None)
-                )
-            except (ValueError, TypeError) as e:
+                sysobj.adjust_electrons(charge=system.pop("charge", None),
+                                        multiplicity=system.pop("multiplicity", None))
+            except ValueError as e:
                 raise ValueError("system/multiplicity or system/charge erroneous: " +
                                  str(e))
 
+            if "n_alpha" in system or "n_beta" in system:
+                warnings.warn("Overriding system/n_alpha and system/n_beta in "
+                              "ScfParameters, since system/multiplicity or "
+                              "system/charge is present.")
+
             system["n_alpha"] = np.uint64(sysobj.n_alpha)
             system["n_beta"] = np.uint64(sysobj.n_beta)
 
-            for k in ["charge", "multiplicity"]:
-                if k in system:
-                    del system[k]
-
-        if "charge" in system:
-            raise ValueError("If system/charge is present, system/multiplicity "
-                             "needs to be present as well.")
-
         if "n_alpha" not in system or "n_beta" not in system:
             raise KeyError("system/n_alpha and system/n_beta or alternatively "
                            "system/multiplicity need to be present.")
@@ -491,7 +485,7 @@ def system(self):
         object represented by the internal parameters.
         """
         self.__normalise_system()
-        return MolecularSystem(
+        return System(
             atoms=self["system/atom_numbers"].value,
             coords=self["system/coords"].value,
             electrons=(self["system/n_alpha"], self["system/n_beta"])