reorganized some code

set alpha/beta occupancy priority in the total electron calculation

python/unittest/test_ammonia.py

@@ -18,14 +18,14 @@ def setUp(self):
                                   multiplicity=1,
                                   group='C3v')
 
+    def test_csm_dens(self):
+        np.testing.assert_allclose(0.14761255, self.structure.csm_dens, rtol=1e-03)
+
     def test_csm_dens_coef(self):
         csm_coef_test = [1.000000, 0.999995, 0.997539, 0.997542, 0.997544]
 
         np.testing.assert_allclose(csm_coef_test, self.structure.csm_dens_coef, rtol=1e-06)
 
-    def test_csm_dens(self):
-        np.testing.assert_allclose(0.14761255, self.structure.csm_dens, rtol=1e-03)
-
     def test_dens_occupancy(self):
         self.structure = WfnSympy(coordinates=self.data['coordinates'],
                                   symbols=self.data['symbols'],
@@ -37,7 +37,7 @@ def test_dens_occupancy(self):
                                   multiplicity=1,
                                   group='C3v',
                                   alpha_occupancy=[0, 1, 1, 1, 1])
-        np.testing.assert_allclose(4.439216, self.structure.csm_dens, rtol=1e-03)
+        np.testing.assert_allclose(4.25514, self.structure.csm_dens, rtol=1e-03)
 
     def test_precision(self):
         self.structure = WfnSympy(coordinates=self.data['coordinates'],

python/wfnsympy/__init__.py

@@ -301,34 +301,47 @@ def __init__(self,
                 for pc in shell['p_con_coefficients']:
                     p_c_coefficients.append(pc)
 
+        self._n_shell = np.unique(atom_map, return_counts=True)[1]
+
         # convert from Angstroms to Bohr
         self._coordinates = np.array(coordinates)
 
         # get atomic numbers
         self._atomic_numbers = [symbol_map[i] for i in symbols]
 
+        # Create MO coefficients in contiguous list
+        self._n_mo = len(alpha_mo_coeff)
+        self._ca = np.ascontiguousarray(alpha_mo_coeff).flatten().tolist()
+        if beta_mo_coeff is not None:
+            self._cb = np.ascontiguousarray(beta_mo_coeff).flatten().tolist()
+        else:
+            self._cb = self._ca
+
+        if self._ca == self._cb:
+            self._unrestricted = False
+        else:
+            self._unrestricted = True
+
         # get valence electrons
         self._valence_electrons = get_valence_electrons(self._atomic_numbers, self._charge)
 
         # get total number of electrons
-        if valence_only:
-            self._total_electrons = self._valence_electrons
+        if self._alpha_occupancy is not None:
+            if self._beta_occupancy is None:
+                self._total_electrons = 2*np.sum(self._alpha_occupancy)
+            else:
+                self._total_electrons = np.sum(self._alpha_occupancy) + np.sum(self._beta_occupancy)
         else:
-            self._total_electrons = np.sum(self._atomic_numbers) - self._charge
+            if valence_only:
+                self._total_electrons = self._valence_electrons
+            else:
+                self._total_electrons = np.sum(self._atomic_numbers) - self._charge
 
         # Check total_electrons compatible with multiplicity
         if (np.remainder(self._total_electrons, 2) == np.remainder(self._multiplicity, 2) or
             self._total_electrons < self._multiplicity):
             raise MultiplicityError(self._multiplicity, self._total_electrons)
 
-        # Create MO coefficients in contiguous list
-        self._n_mo = len(alpha_mo_coeff)
-        self._ca = np.ascontiguousarray(alpha_mo_coeff).flatten().tolist()
-        if beta_mo_coeff is not None:
-            self._cb = np.ascontiguousarray(beta_mo_coeff).flatten().tolist()
-        else:
-            self._cb = self._ca
-
         # Check if electrons fit in provided MO
         if (self._total_electrons + self._multiplicity - 1)/2 > self._n_mo:
             self._total_electrons = self._n_mo * 2
@@ -337,14 +350,12 @@ def __init__(self,
         if self._valence_electrons >= self._total_electrons:
             self._valence_electrons = 0
 
-        if self._ca == self._cb:
-            self._unrestricted = False
-        else:
-            self._unrestricted = True
-
         if self._alpha_occupancy is None:
             self._alpha_occupancy = [0]*int(self._n_mo)
             self._alpha_occupancy[:int(self._total_electrons//2)] = [1]*int(self._total_electrons//2)
+            if self._multiplicity > 1:
+                for i in range(self._multiplicity - 1):
+                    self._alpha_occupancy[int(self._total_electrons // 2) + i] = 1
         else:
             if len(self._alpha_occupancy) != self._n_mo:
                 for _ in range(int(self._n_mo - len(self._alpha_occupancy))):
@@ -358,11 +369,6 @@ def __init__(self,
                 for _ in range(int(self._n_mo - len(self._beta_occupancy))):
                     self._beta_occupancy.append(0)
 
-        if self._multiplicity > 1:
-            for i in range(self._multiplicity-1):
-                self._alpha_occupancy[int(self._total_electrons//2)+i] = 1
-        self._n_shell = np.unique(atom_map, return_counts=True)[1]
-
         # Transform symbols type to correct Fortran char*2 type
         self._symbols = np.array([list('{:<2}'.format(char)) for char in symbols], dtype='S')