Merge branch 'development_efrem'

cosymlib/file_io/__init__.py

@@ -119,10 +119,10 @@ def get_geometry_from_file_cor(file_name, read_multiple=False):
 
 # TODO: This function should handle open and close shell, mo energies and more. A lot to improve!
 def get_molecule_from_file_fchk(file_name, read_multiple=False):
-    key_list = ['Charge', 'Multiplicity', '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 Orbital Energies', 'Alpha MO coefficients', 'Beta MO coefficients', ]
+    key_list = ['Charge', 'Multiplicity', '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 Orbital Energies', 'Alpha MO coefficients', 'Beta MO coefficients']
     input_molecule = [[] for _ in range(len(key_list))]
     read = False
     with open(file_name, mode='r') as lines:
@@ -146,7 +146,7 @@ def get_molecule_from_file_fchk(file_name, read_multiple=False):
                 if key in line:
                     if n == len(key_list) - 1:
                         break
-                    if options and idn != 2:
+                    if options and idn != 4:
                         input_molecule[idn].append(int(line.split()[-1]))
                         n = idn
                         break
@@ -159,42 +159,44 @@ def get_molecule_from_file_fchk(file_name, read_multiple=False):
                     read = True
                     break
 
-        for n in range(2, len(input_molecule) - 1):
+        for n in range(4, len(input_molecule)):
             input_molecule[n] = reformat_input(input_molecule[n])
         bohr_to_angstrom = 0.529177249
-        coordinates = np.array(input_molecule[3], dtype=float).reshape(-1, 3) * bohr_to_angstrom
+        coordinates = np.array(input_molecule[5], dtype=float).reshape(-1, 3) * bohr_to_angstrom
 
-        atomic_number = [int(num) for num in input_molecule[2]]
+        atomic_number = [int(num) for num in input_molecule[4]]
         symbols = []
         for number in atomic_number:
             symbols.append(atomic_number_to_element(number))
 
         basis = basis_format(basis_set_name=basis_set,
                              atomic_numbers=atomic_number,
                              atomic_symbols=symbols,
-                             shell_type=input_molecule[4],
-                             n_primitives=input_molecule[5],
-                             atom_map=input_molecule[6],
-                             p_exponents=input_molecule[7],
-                             c_coefficients=input_molecule[8],
-                             p_c_coefficients=input_molecule[9])
+                             shell_type=input_molecule[6],
+                             n_primitives=input_molecule[7],
+                             atom_map=input_molecule[8],
+                             p_exponents=input_molecule[9],
+                             c_coefficients=input_molecule[10],
+                             p_c_coefficients=input_molecule[11])
 
         geometry = Geometry(symbols=symbols,
                             positions=coordinates,
                             name=name)
 
-        Ca = np.array(input_molecule[11], dtype=float).reshape(-1, int(np.sqrt(len(input_molecule[11]))))
-        energies_alpha = np.array(input_molecule[10], dtype=float).tolist()
-        if input_molecule[12]:
-            Cb = np.array(input_molecule[12], dtype=float).reshape(-1, int(np.sqrt(len(input_molecule[12]))))
+        Ca = np.array(input_molecule[13], dtype=float).reshape(-1, int(np.sqrt(len(input_molecule[13]))))
+        energies_alpha = np.array(input_molecule[12], dtype=float).tolist()
+        if input_molecule[14]:
+            Cb = np.array(input_molecule[14], dtype=float).reshape(-1, int(np.sqrt(len(input_molecule[14]))))
         else:
             Cb = []
 
         electronic_structure = ElectronicStructure(charge=input_molecule[0][0],
                                                    multiplicity=input_molecule[1][0],
                                                    basis=basis,
                                                    orbital_coefficients=[Ca, Cb],
-                                                   mo_energies=energies_alpha)
+                                                   mo_energies=energies_alpha,
+                                                   alpha_electrons=input_molecule[2],
+                                                   beta_electrons=input_molecule[3])
 
         if read_multiple:
             return [Molecule(geometry, electronic_structure)]

cosymlib/molecule/__init__.py

@@ -46,7 +46,8 @@ def electronic_structure(self):
                                                              orbital_coefficients=[eh.get_mo_coefficients(), []],
                                                              mo_energies=eh.get_mo_energies(),
                                                              multiplicity=eh.get_multiplicity(),
-                                                             valence_only=True)
+                                                             alpha_electrons=[eh.get_alpha_electrons()],
+                                                             beta_electrons=[eh.get_beta_electrons()])
             self.symmetry.set_electronic_structure(self._electronic_structure)
 
         return self._electronic_structure

cosymlib/molecule/electronic_structure/__init__.py

@@ -5,20 +5,27 @@ def __init__(self,
                  basis=None,
                  orbital_coefficients=None,
                  mo_energies=None,
-                 valence_only=False):
+                 alpha_electrons=None,
+                 beta_electrons=None):
 
         self._charge = charge
         self._multiplicity = multiplicity
         self._basis = basis
-        self._valence_only = valence_only
+        # self._valence_only = valence_only
         self._Ca = orbital_coefficients[0]
-        if not orbital_coefficients[1]:
+        if len(orbital_coefficients[1]) == 0:
             self._Cb = orbital_coefficients[0]
         else:
             self._Cb = orbital_coefficients[1]
 
         self._mo_energies = mo_energies
 
+        self._alpha_occupancy = [1] * int(alpha_electrons[0])
+        self._beta_occupancy = [1] * int(beta_electrons[0])
+        if self._multiplicity > 1:
+            for i in range(self._multiplicity-1):
+                self._beta_occupancy.append(0)
+
     @property
     def charge(self):
         return self._charge
@@ -43,6 +50,22 @@ def coefficients_b(self):
     def energies(self):
         return self._mo_energies
 
+    # @property
+    # def valence_only(self):
+    #     return self._valence_only
+
+    @property
+    def alpha_occupancy(self):
+        return self._alpha_occupancy
+
+    @property
+    def beta_occupancy(self):
+        return self._beta_occupancy
+
+    @property
+    def alpha_electrons(self):
+        return sum(self._alpha_occupancy)
+
     @property
-    def valence_only(self):
-        return self._valence_only
+    def beta_electrons(self):
+        return sum(self._beta_occupancy)

cosymlib/simulation/__init__.py

@@ -5,6 +5,9 @@ class ExtendedHuckel:
 
     def __init__(self, geometry):
         self._EH = huckelpy.ExtendedHuckel(geometry.get_positions(), geometry.get_symbols())
+        self._alpha_electrons = None
+        self._beta_electrons = None
+        self._total_electrons = self._EH.get_number_of_electrons()
 
     def get_mo_coefficients(self):
         return self._EH.get_eigenvectors()
@@ -18,6 +21,14 @@ def get_mo_energies(self):
     def get_multiplicity(self):
         return self._EH.get_multiplicity()
 
+    def get_alpha_electrons(self):
+        if self._alpha_electrons is None:
+            self._alpha_electrons = self._total_electrons // 2 + self.get_multiplicity() - 1
+        return self._alpha_electrons
+
+    def get_beta_electrons(self):
+        return self._total_electrons - self._alpha_electrons
+
     def build_fchk_file(self, name):
         txt_fchk = huckelpy.file_io.build_fchk(self._EH)
         open(name + '.fchk', 'w').write(txt_fchk)

cosymlib/symmetry/__init__.py

@@ -115,7 +115,8 @@ def _get_wfnsym_results(self, group):
                                           charge=self._electronic_structure.charge,
                                           multiplicity=self._electronic_structure.multiplicity,
                                           group=group.upper(),
-                                          valence_only=self._electronic_structure.valence_only)
+                                          alpha_occupancy=self._electronic_structure.alpha_occupancy,
+                                          beta_occupancy=self._electronic_structure.beta_occupancy)
         return self._results[key]
 
     ##########################################

examples/normal_modes_api.py

@@ -0,0 +1,59 @@
+import numpy as np
+from cosymlib.file_io import read_generic_structure_file, get_file_xyz_txt
+from cosymlib.molecule.geometry import Geometry
+
+
+def read_normal_modes_gaussian(file_name):
+    read = False
+    frequencies = []
+    normal_modes_matrices = []
+    with open(file_name, mode='r') as lines:
+        for line in lines:
+            if read:
+                if len(line.split()) == 3 or not line.strip():
+                    read = False
+                else:
+                    for i in range(int(len(line.split()[2:])/3)):
+                        atom_vibration = line.split()[2 + 3*i:3 + 2 + 3*i]
+                        normal_modes_matrices[-3+i].append([float(x)for x in atom_vibration])
+
+            if 'Atom  AN' in line:
+                read = True
+                for _ in range(line.split()[2:].count('X')):
+                    normal_modes_matrices.append([])
+            if 'Frequencies --' in line:
+                frequencies.append(line.split()[2:])
+    frequencies = [float(frequency) for sublist in frequencies for frequency in sublist]
+
+    return frequencies, normal_modes_matrices
+
+
+def get_nm_vibration_path(geometry, normal_mode, points=10, backward=False, k=1.5):
+
+    nm_np_matrix = np.array(normal_mode)
+    geom_np = np.array(geometry)
+    if backward:
+        x = np.linspace(-k, k, points)
+    else:
+        x = np.linspace(k/points, k, points)
+
+    structures_path = []
+    for dx in x:
+        structures_path.append(geom_np + dx*nm_np_matrix)
+
+    return structures_path
+
+
+molecule = read_generic_structure_file('sf6.fchk')
+freq, nm_martices = read_normal_modes_gaussian('sf6_freq.out')
+k_points = 1.5
+n_freq = 0
+path = get_nm_vibration_path(molecule.geometry.get_positions(), nm_martices[0], k=k_points)
+
+x = np.linspace(k_points/10, k_points, 10)
+output = open('sf6_freq{}.xyz'.format(n_freq+1), 'w')
+for ids, structure in enumerate(path):
+    xi = '{:.2f}'.format(x[ids]).replace('.', '')
+    txt = get_file_xyz_txt(Geometry(structure, symbols=molecule.geometry.get_symbols(), name='freq : ' + str(freq[0]) +
+                                                                                             ' ' + str(xi)))
+    output.write(txt)