charge and multiplicity variables useless in class Wfnsympy

README.md

@@ -76,8 +76,8 @@ wf_results = WfnSympy(coordinates=[[ 0.00000000, 0.00000000, -0.04280085],
                       symbols=['O', 'H', 'H'],
                       basis=basis,
                       alpha_mo_coeff=mo_coefficients,
-                      charge=0,
-                      multiplicity=1,
+                      # charge=0,
+                      # multiplicity=1,
                       group='C2v')
 
 wf_results.print_CSM()

examples/python_api.py

@@ -10,10 +10,12 @@
                      basis=data['basis'],
                      alpha_mo_coeff=data['mo_coefficients']['alpha'],
                      # center=[0., 0., 0.],
-                     # VAxis=[0., 0., 1.],
-                     # VAxis2=[0, 1, 0],
-                     charge=0,
-                     multiplicity=1,
+                     # axis=[0., 0., 1.],
+                     # axis2=[0, 1, 0],
+                     alpha_occupancy=data['alpha_occupancy'],
+                     beta_occupancy=data['beta_occupancy'],
+                     # charge=0,
+                     # multiplicity=1,
                      group='Td')
 
 

python/WFNSYMLIB.pyf

@@ -3,13 +3,14 @@
 
 python module WFNSYMLIB ! in 
     interface  ! in :WFNSYMLIB
-        subroutine mainlib(Etot, NBas, NMo, Norb, NAt, NCShell, iZAt, AtLab, Alph, COrb, NShell, RAt, n_prim, shell_type, &
-        igroup, ngroup, Ca, Cb, RCread, VAxis, VAxis2, iMult, DoOp, OutDim, OutGrim, OutCMS, OutSymLab, OutSDiagA, &
-        OutSDiagB, OutWDiagA, OutWDiagB, OutTbl, OutIRLab, OutgIRA, OutgIRB, OutgIRwfA, OutgIRwfB, OutgIRwf, OutSymMat) ! in :WFNSYMLIB:lib_main.F
+        subroutine mainlib(AOccup, BOccup, NBas, NMo, Norb, NAt, NCShell, iZAt, AtLab, Alph, COrb, NShell, RAt, n_prim, &
+        shell_type, igroup, ngroup, Ca, Cb, RCread, VAxis, VAxis2, nAOccup, nBOccup, DoOp, &
+        OutDim, OutGrim, OutCMS, OutSymLab, OutSDiagA, OutSDiagB, OutWDiagA, OutWDiagB, OutTbl, OutIRLab, OutgIRA, &
+        OutgIRB, OutgIRwfA, OutgIRwfB, OutgIRwf, OutSymMat) ! in :WFNSYMLIB:lib_main.F
 
 !           Inputs
             logical :: DoOp
-            integer :: Etot, NBas, Norb, NAt, NCShell, igroup, ngroup, iMult, NMo
+            integer :: NBas, Norb, NAt, NCShell, igroup, ngroup, NMo
 
             ! integer, check(len(Ca)/NBas>=NMo),depend(Ca), intent(hide) :: NMo=len(Ca)/NBas
             ! integer, intent(hide) :: NOrb_w
@@ -21,6 +22,10 @@ python module WFNSYMLIB ! in
             integer dimension(NAt),intent(in),depend(NAt) :: iZAt, NShell
             integer dimension(NCShell),intent(in),depend(NCShell) :: n_prim, shell_type
             character dimension(NAt),intent(in),depend(NAt) :: AtLab*2
+            real*8 dimension(nAOccup), intent(in) :: AOccup
+            real*8 dimension(nBOccup), intent(in) :: BOccup
+            integer, optional,check(len(AOccup)>=nAOccup),depend(AOccup) :: nAOccup=len(AOccup)
+            integer, optional,check(len(BOccup)>=nBOccup),depend(BOccup) :: nBOccup=len(BOccup)
 
 !           Outputs
             integer dimension(3),intent(out) :: OutDim

python/unittest/test_ammonia.py

@@ -11,12 +11,15 @@ def setUp(self):
         self.structure = WfnSympy(coordinates=self.data['coordinates'],
                                   symbols=self.data['symbols'],
                                   basis=self.data['basis'],
-                                  VAxis=[0.55944897, 0.69527034, -0.4512383 ],
-                                  VAxis2=[0., -0.4512383, -0.69527034],
+                                  axis=[0.55944897, 0.69527034, -0.4512383],
+                                  axis2=[0., -0.4512383, -0.69527034],
                                   alpha_mo_coeff=self.data['mo_coefficients']['alpha'],
-                                  charge=0,
-                                  multiplicity=1,
-                                  group='C3v')
+                                  # charge=0,
+                                  # multiplicity=1,
+                                  group='C3v',
+                                  # alpha_occupancy=self.data['alpha_occupancy'],
+                                  # beta_occupancy=self.data['beta_occupancy']
+                                  )
 
     def test_csm_dens(self):
         np.testing.assert_allclose(0.148, self.structure.csm_dens, rtol=1e-03)
@@ -30,11 +33,11 @@ def test_dens_occupancy(self):
         self.structure = WfnSympy(coordinates=self.data['coordinates'],
                                   symbols=self.data['symbols'],
                                   basis=self.data['basis'],
-                                  VAxis=[0.55944897, 0.69527034, -0.4512383 ],
-                                  VAxis2=[0., -0.4512383, -0.69527034],
+                                  axis=[0.55944897, 0.69527034, -0.4512383],
+                                  axis2=[0., -0.4512383, -0.69527034],
                                   alpha_mo_coeff=self.data['mo_coefficients']['alpha'],
-                                  charge=0,
-                                  multiplicity=1,
+                                  # charge=0,
+                                  # multiplicity=1,
                                   group='C3v',
                                   alpha_occupancy=[0, 1, 1, 1, 1])
         np.testing.assert_allclose(4.507, self.structure.csm_dens, rtol=1e-03)
@@ -43,11 +46,13 @@ def test_precision(self):
         self.structure = WfnSympy(coordinates=self.data['coordinates'],
                                   symbols=self.data['symbols'],
                                   basis=self.data['basis'],
-                                  VAxis=[0.55944897, 0.69527034, -0.4512383],
-                                  VAxis2=[0., -0.4512383, -0.69527034],
+                                  axis=[0.55944897, 0.69527034, -0.4512383],
+                                  axis2=[0., -0.4512383, -0.69527034],
                                   alpha_mo_coeff=self.data['mo_coefficients']['alpha'],
-                                  charge=0,
-                                  multiplicity=1,
+                                  # charge=0,
+                                  # multiplicity=1,
                                   group='C3v',
-                                  tolerance=1e-04)
+                                  tolerance=1e-04,
+                                  alpha_occupancy=self.data['alpha_occupancy'],
+                                  beta_occupancy=self.data['beta_occupancy'])
         np.testing.assert_allclose(0.148, self.structure.csm_dens, rtol=1e-03)

python/unittest/test_pirrol.py

@@ -12,11 +12,13 @@ def setUp(self):
                                symbols=data['symbols'],
                                basis=data['basis'],
                                center=[0., 0., 0.],
-                               VAxis=[1., 0., 0.],
-                               # VAxis2=[0., 0., 1.],
+                               axis=[1., 0., 0.],
+                               # axis2=[0., 0., 1.],
                                alpha_mo_coeff=data['mo_coefficients']['alpha'][:18],
-                               charge=0,
-                               multiplicity=1,
+                               alpha_occupancy=data['alpha_occupancy'],
+                               beta_occupancy=data['beta_occupancy'],
+                               # charge=0,
+                               # multiplicity=1,
                                group='C6v')
 
     def test_symlab(self):

python/unittest/test_water.py

@@ -43,8 +43,8 @@ def setUp(self):
                                    symbols=['O', 'H', 'H'],
                                    basis=basis,
                                    alpha_mo_coeff=mo_coefficients[:5],
-                                   charge=0,
-                                   multiplicity=1,
+                                   # charge=0,
+                                   # multiplicity=1,
                                    group='C2v')
 
         self.wf_results.print_alpha_mo_IRD()

python/wfnsympy/__init__.py

@@ -41,22 +41,22 @@ def __exit__(self, exc_type, exc_value, traceback):
         self.F.close()
 
 
-def get_valence_electrons(atomic_numbers, charge):
-    valence_electrons = 0
-    for number in atomic_numbers:
-        if 2 >= number > 0:
-            valence_electrons += np.mod(number, 2)
-        if 18 >= number > 2:
-            valence_electrons += np.mod(number - 2, 8)
-        if 54 >= number > 18:
-            valence_electrons += np.mod(number - 18, 18)
-        if 118 >= number > 54:
-            valence_electrons += np.mod(number - 54, 32)
-        if number > 118:
-            raise Exception('Atomic number size not implemented')
-
-    valence_electrons -= charge
-    return valence_electrons
+# def get_valence_electrons(atomic_numbers, charge):
+#     valence_electrons = 0
+#     for number in atomic_numbers:
+#         if 2 >= number > 0:
+#             valence_electrons += np.mod(number, 2)
+#         if 18 >= number > 2:
+#             valence_electrons += np.mod(number - 2, 8)
+#         if 54 >= number > 18:
+#             valence_electrons += np.mod(number - 18, 18)
+#         if 118 >= number > 54:
+#             valence_electrons += np.mod(number - 54, 32)
+#         if number > 118:
+#             raise Exception('Atomic number size not implemented')
+#
+#     valence_electrons -= charge
+#     return valence_electrons
 
 
 def _get_group_num_from_label(label):
@@ -280,14 +280,13 @@ def __init__(self,
                  basis,  # basis dictionary
                  alpha_mo_coeff,  # Nbas x Nbas
                  center=None,  # in Angstrom
-                 VAxis=None,
-                 VAxis2=None,
-                 charge=0,
-                 multiplicity=1,
+                 axis=None,
+                 axis2=None,
+                 # charge=0,
+                 # multiplicity=1,
                  beta_mo_coeff=None,  # Nbas x Nbas
                  group=None,
                  do_operation=False,
-                 # valence_only=False,
                  alpha_occupancy=None,
                  beta_occupancy=None,
                  tolerance=1e-8):
@@ -302,10 +301,10 @@ def __init__(self,
 
         self._do_operation = do_operation
         self._center = center
-        self._axis = VAxis
-        self._axis2 = VAxis2
-        self._charge = charge
-        self._multiplicity = multiplicity
+        self._axis = axis
+        self._axis2 = axis2
+        # self._charge = charge
+        # self._multiplicity = multiplicity
         self._alpha_occupancy = alpha_occupancy
         self._beta_occupancy = beta_occupancy
         self._toldens = tolerance
@@ -358,39 +357,35 @@ def __init__(self,
         else:
             self._unrestricted = True
 
-        # get valence electrons
-        # self._valence_electrons = get_valence_electrons(self._atomic_numbers, self._charge)
-
         # get total number of 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:
-        #     if valence_only:
-        #         self._total_electrons = self._valence_electrons
-        #     else:
-            self._total_electrons = np.sum(self._atomic_numbers) - self._charge
+            # self._total_electrons = np.sum(self._atomic_numbers) - self._charge
+            self._total_electrons = np.sum(self._atomic_numbers)
             # 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)
+            # 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)
 
         # Check if electrons fit in provided MO
-        if (self._total_electrons + self._multiplicity - 1)/2 > self._n_mo:
+        # if (self._total_electrons + self._multiplicity - 1)/2 > self._n_mo:
+        if self._total_electrons/2 > self._n_mo:
             self._total_electrons = self._n_mo * 2
-            self._multiplicity = 1
-
-        # if self._valence_electrons >= self._total_electrons:
-        #     self._valence_electrons = 0
+            # self._multiplicity = 1
 
+        # self._total_electrons += 1
         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
+            if self._total_electrons%2 != 0:
+                self._alpha_occupancy[int(self._total_electrons//2)] = 1
+            # 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))):
@@ -518,12 +513,11 @@ def target_function(gamma, VAxis):
         # Start calculation
         with _captured_stdout() as E:
             coordinates_bohr = np.array(self._coordinates) / _bohr_to_angstrom
-            out_data = mainlib(self._total_electrons, self._n_bas, self._n_mo,
+            out_data = mainlib(self._alpha_occupancy, self._beta_occupancy, self._n_bas, self._n_mo,
                                self._n_uncontr_orbitals, self._n_atoms, self._ntot_shell, self._atomic_numbers,
                                self._symbols, self._alpha, self._uncontracted_coefficients, self._n_shell,
                                coordinates_bohr, self._n_primitives, self._shell_type, self._igroup, self._ngroup,
-                               self._ca, self._cb, self._center, self._axis, self._axis2,
-                               self._multiplicity, self._do_operation)
+                               self._ca, self._cb, self._center, self._axis, self._axis2, self._do_operation)
             E.seek(0)
             capture = E.read()
         capture = capture.decode().split('\n')

python/wfnsympy/file_io.py

@@ -71,7 +71,7 @@ def basis_format(basis_set_name,
 
 
 def get_data_from_file_fchk(file_name):
-    key_list = ['Charge', 'Multiplicity', 'Atomic numbers', 'Current cartesian coordinates',
+    key_list = ['Number of alpha electrons', 'Number of beta electrons', 'Atomic numbers', 'Current cartesian coordinates',
                 'Shell type', 'Number of primitives per shell', 'Shell to atom map', 'Primitive exponents',
                 'Contraction coefficients', 'P(S=P) Contraction coefficients',
                 'Alpha MO coefficients', 'Beta MO coefficients']
@@ -126,6 +126,12 @@ def get_symbols_from_atomic_numbers(atomic_numbers):
         atomic_numbers = [int(num) for num in input_molecule[2]]
         symbols = get_symbols_from_atomic_numbers(atomic_numbers)
 
+        alpha_occupancy = [1] * int(input_molecule[0][0])
+        beta_occupancy = [1] * int(input_molecule[1][0])
+        if len(alpha_occupancy) != len(beta_occupancy):
+            for i in range(abs(len(alpha_occupancy) - len(beta_occupancy))):
+                beta_occupancy.append(0)
+
         basis = basis_format(basis_set_name=basis_set,
                              atomic_numbers=atomic_numbers,
                              atomic_symbols=symbols,
@@ -145,7 +151,9 @@ def get_symbols_from_atomic_numbers(atomic_numbers):
         return {'coordinates': coordinates,
                 'symbols': symbols,
                 'basis': basis,
-                'mo_coefficients': coefficients}
+                'mo_coefficients': coefficients,
+                'alpha_occupancy': alpha_occupancy,
+                'beta_occupancy': beta_occupancy}
 
 
 if __name__ == '__main__':

python/wfnsympy/optimize.py

@@ -127,11 +127,11 @@ def target_function(alpha, beta):
                           symbols=data['symbols'],
                           basis=data['basis'],
                           center=c_geom,
-                          VAxis=ax1,
-                          # VAxis2=ax2,
+                          axis=ax1,
+                          # axis2=ax2,
                           alpha_mo_coeff=data['mo_coefficients']['alpha'],
-                          charge=0,
-                          multiplicity=1,
+                          # charge=0,
+                          # multiplicity=1,
                           group='C5')
         res = np.sum(np.sqrt(pirrol.csm_coef))
         return res

src/lib_main.F

@@ -1,10 +1,10 @@
 #include "wfnsym.h"
 
 c      PROGRAM Prog_Main
-      SUBROUTINE mainlib(NEtot, NBas, NMo, Norb,
+      SUBROUTINE mainlib(AOccup, BOccup, NBas, NMo, Norb,
      & NAt_2, NTotShell, iZAt, AtLab, Alph, COrb, NShell,
      & RAt, NPrim, shell_type, igroup_2, ngroup_2, Ca, Cb,
-     & RCread_2, VAxis_2, VAxis2_2, iMult, DoOp,
+     & RCread_2, VAxis_2, VAxis2_2, nAOccup, nBOccup, DoOp,
      & OutDim, OutGrim, OutCSM, OutSymLab,
      & OutSDiagA, OutSDiagB, OutWDiagA, OutWDiagB, OutTbl,
      & OutIRLab, OutgIRA, OutgIRB,OutgIRwfA,OutgIRwfB,
@@ -52,11 +52,12 @@ SUBROUTINE mainlib(NEtot, NBas, NMo, Norb,
 c					 True : 5 decimals
 c--------------------------------------------------------------
       IMPLICIT REAL*8 (a-h,o-z)
-      INTEGER :: iZAt(NAt_2), NPrim(NTotShell)
+      INTEGER :: iZAt(NAt_2), NPrim(NTotShell), nAOccup, nBOccup
       CHARACTER :: AtLab(NAt_2)*2
       REAL*8 :: Alph(Norb),  COrb(Norb), RAt(NAt_2, 3)
       REAL*8 :: CA(NMo*NBas),  CB(NMo*NBas)
       REAL*8 :: RCread_2(3), VAxis_2(3), VAxis2_2(3)
+      REAL*8 :: AOccup(nAOccup), BOccup(nBOccup)
       LOGICAL :: DoOp
 c     Output
       INTEGER :: OutDim(3)
@@ -137,16 +138,19 @@ SUBROUTINE mainlib(NEtot, NBas, NMo, Norb,
 
 
 c Set up some variables
-      NOa = (NEtot + iMult - 1)/2
-      NOb = NOa - iMult + 1
+
+      NOa = SUM(AOccup)
+      NOb = SUM(BOccup)
+c      NOa = (NEtot + iMult - 1)/2
+c      NOb = NOa - iMult + 1
 c      Ncore = (NEtot-NEval)/2
 c      NOva = NOa - Ncore
 c      NOvb = NOb - Ncore
       NBas2 = NBas*NBas
-      IF(NOa+NOb .NE. NEtot) THEN
-        WRITE(iout, *) 'ERROR. Charge/Multiplicity inconsistency'
-        STOP
-      ENDIF
+c      IF(NOa+NOb .NE. NEtot) THEN
+c        WRITE(iout, *) 'ERROR. Charge/Multiplicity inconsistency'
+c        STOP
+c      ENDIF
 
 
 c Where do we center the molecule?