Merge branch 'development_efrem'

cosymlib/__init__.py

@@ -1,4 +1,4 @@
-__version__ = '0.9.1'
+__version__ = '0.9.2'
 
 from cosymlib.molecule import Molecule
 from cosymlib.molecule.geometry import Geometry
@@ -10,6 +10,7 @@
 import matplotlib.pyplot as plt
 
 import sys
+import os
 import numpy as np
 
 
@@ -59,7 +60,8 @@ class Cosymlib:
     """
     Main cosymlib class
 
-    :param structures: List of :class:`~cosymlib.molecule.geometry.Geometry` or :class:`~cosymlib.molecule.Molecule` type objects
+    :param structures: List of :class:`~cosymlib.molecule.geometry.Geometry` or :class:`~cosymlib.molecule.Molecule`
+    type objects
     :type structures: list
     :param ignore_atoms_labels: Ignore atomic element labels is symmetry calculations
     :type ignore_atoms_labels: bool
@@ -156,6 +158,8 @@ def print_shape_measure(self, shape_reference, central_atom=0, fix_permutation=F
         :type shape_reference: str
         :param central_atom: Position of the central atom
         :type central_atom: int
+        :param fix_permutation: Do not permute atoms during shape calculations
+        :type fix_permutation: bool
         :param output: Display hook
         :type output: hook
         """
@@ -212,13 +216,15 @@ def print_shape_structure(self, shape_reference, central_atom=0, fix_permutation
 
         shape_results_structures = []
         references = []
+        sym_labels = []
         for reference in reference_list:
             if isinstance(reference, str):
                 references.append(reference)
-
+                sym_labels.append(get_sym_from_label(reference))
             else:
                 references.append(reference.name)
-                reference = reference.get_positions()
+                sym_labels.append('')
+                # reference = reference.get_positions()
 
             shape_results_structures.append(self.get_shape_measure(reference,
                                                                    'structure',
@@ -229,7 +235,7 @@ def print_shape_structure(self, shape_reference, central_atom=0, fix_permutation
             geometries.append(molecule.geometry)
             for idl, reference in enumerate(references):
                 shape_results_structures[idl][idm].set_name(molecule.name + ' ' + reference + ' ' +
-                                                            get_sym_from_label(reference))
+                                                            sym_labels[idl])
                 geometries.append(shape_results_structures[idl][idm])
 
         print("\nOriginal structures vs reference polyhedra in file {}\n".format(output.name))
@@ -267,7 +273,7 @@ def print_geometric_measure_info(self, label, multi=1, central_atom=0, center=No
             for idn, array in enumerate(molecule.geometry.get_positions()):
                 array = array - cm
                 txt += '{:2} {:12.8f} {:12.8f} {:12.8f}\n'.format(molecule.geometry.get_symbols()[idn],
-                                                                      array[0], array[1], array[2])
+                                                                  array[0], array[1], array[2])
             txt += sep_line
 
             txt += 'Optimal permutation\n'
@@ -278,7 +284,7 @@ def print_geometric_measure_info(self, label, multi=1, central_atom=0, center=No
             txt += 'Inverted structure\n'
             for idn, axis in enumerate(molecule.geometry._symmetry.nearest_structure(label)):
                 txt += '{:2} {:12.8f} {:12.8f} {:12.8f}\n'.format(molecule.geometry.get_symbols()[idn],
-                                                                      axis[0], axis[1], axis[2])
+                                                                  axis[0], axis[1], axis[2])
             txt += '\n'
 
             txt += 'Reference axis\n'
@@ -373,6 +379,14 @@ def print_electronic_symmetry_measure(self, group, axis=None, axis2=None, center
             txt += '{:<9} '.format(molecule.name) + '  '.join(['{:7.3f}'.format(s) for s in wf_measure['csm']])
             txt += '\n'
             first = False
+        axes_information = molecule.get_symmetry_axes(group, axis=axis, axis2=axis2, center=center)
+        txt += '\nSymmetry parameters\n'
+        txt += 'center: ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['center']])
+        txt += '\n'
+        txt += 'axis  : ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['axis']])
+        txt += '\n'
+        txt += 'axis2 : ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['axis2']])
+        txt += '\n'
 
         output.write(txt)
 
@@ -402,7 +416,8 @@ def print_electronic_density_measure(self, group, axis=None, axis2=None, center=
             first = False
             axes_information = molecule.get_symmetry_axes(group, axis=axis, axis2=axis2, center=center)
             txt2 += '{:<9} '.format(molecule.name) + '{:7.3f}\n'.format(dens_measure['csm'])
-            txt2 += '\ncenter: ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['center']])
+            txt2 += '\nSymmetry parameters\n'
+            txt2 += 'center: ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['center']])
             txt2 += '\n'
             txt2 += 'axis  : ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['axis']])
             txt2 += '\n'
@@ -432,7 +447,7 @@ def print_wnfsym_sym_matrices(self, group, axis=None, axis2=None, center=None, o
             txt += sep_line
             for idn, array in enumerate(geometry.get_positions()):
                 txt += '{:2} {:11.6f} {:11.6f} {:11.6f}\n'.format(geometry.get_symbols()[idn],
-                                                                      array[0], array[1], array[2])
+                                                                  array[0], array[1], array[2])
             txt += sep_line
             for i, group in enumerate(sym_lables):
                 txt += '\n'
@@ -446,7 +461,7 @@ def print_wnfsym_sym_matrices(self, group, axis=None, axis2=None, center=None, o
                 for idn, array in enumerate(geometry.get_positions()):
                     array2 = np.dot(array, sym_matrices[i].T)
                     txt += '{:2} {:11.6f} {:11.6f} {:11.6f}\n'.format(geometry.get_symbols()[idn],
-                                                                          array2[0], array2[1], array2[2])
+                                                                      array2[0], array2[1], array2[2])
         output.write(txt)
 
     def print_wnfsym_sym_ovelap(self, group, axis=None, axis2=None, center=None, output=sys.stdout):
@@ -559,6 +574,15 @@ def print_wnfsym_irreducible_repr(self, group, axis=None, axis2=None, center=Non
             txt += 'WFN ' + '  '.join(['{:7.3f}'.format(s) for s in data_wf['total']])
             txt += '\n'
 
+        axes_information = molecule.get_symmetry_axes(group, axis=axis, axis2=axis2, center=center)
+        txt += '\nSymmetry parameters\n'
+        txt += 'center: ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['center']])
+        txt += '\n'
+        txt += 'axis  : ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['axis']])
+        txt += '\n'
+        txt += 'axis2 : ' + '  '.join(['{:12.8f}'.format(s) for s in axes_information['axis2']])
+        txt += '\n'
+
         output.write(txt)
 
     def plot_mo_diagram(self, group, axis=None, axis2=None, center=None):
@@ -733,23 +757,26 @@ def print_minimum_distortion_path_shape(self, shape_label1, shape_label2, centra
         :type max_gco: float
         :param num_points: Number of points
         :type num_points: int
-        :param output: Display hook
-        :type output: hook
+        :param output1: Display hook
+        :type output1: hook
         """
 
         if output is not None:
-            output = open(output + '_pth.csv', 'w')
-            output2 = open(output + '_pth.xyz', 'w')
-            output3 = open(output, 'w')
+            output_name, file_extension = os.path.splitext(output.name)
+            output1 = open(output_name + '_pth.csv', 'w')
+            output2 = open(output_name + '_pth.xyz', 'w')
+            output3 = output
         else:
-            output = sys.stdout
+            output1 = sys.stdout
             output2 = sys.stdout
             output3 = sys.stdout
 
         csm, devpath, gen_coord = self.get_path_parameters(shape_label1, shape_label2, central_atom=central_atom)
 
-        txt_params = 'Deviation threshold to calculate Path deviation function: {:2.1f}% - {:2.1f}%\n'.format(min_dev, max_dev)
-        txt_params += 'Deviation threshold to calculate Generalized Coordinate: {:2.1f}% - {:2.1f}%\n'.format(min_gco, max_gco)
+        txt_params = 'Deviation threshold to calculate Path deviation function: ' \
+                     '{:2.1f}% - {:2.1f}%\n'.format(min_dev, max_dev)
+        txt_params += 'Deviation threshold to calculate Generalized Coordinate: ' \
+                      '{:2.1f}% - {:2.1f}%\n'.format(min_gco, max_gco)
         txt_params += '\n'
         txt_params += '{:9} '.format('structure'.upper())
         for csm_label in list(csm.keys()):
@@ -769,6 +796,7 @@ def print_minimum_distortion_path_shape(self, shape_label1, shape_label2, centra
             txt_params += '{:^8.1f} {:^8.1f}'.format(devpath[idx], gen_coord[idx])
             txt_params += '\n'
 
+        # self.print_shape_measure()
         txt_params += 'skipped {} structure/s\n\n'.format(filter_mask.count(False))
         output3.write(txt_params)
 
@@ -788,15 +816,15 @@ def print_minimum_distortion_path_shape(self, shape_label1, shape_label2, centra
             txt_path += '{:6.3f}  {:6.3f}'.format(path[0][idx], path[1][idx])
             txt_path += '\n'
         txt_path += '\n'
-        output.write(txt_path)
+        output1.write(txt_path)
 
         test_structures = []
         for ids, structure in enumerate(path[2]):
             test_structures.append(Geometry(symbols=['' for _ in range(len(structure))],
                                             positions=structure, name='map_structure{}'.format(ids)))
         output2.write(file_io.get_file_xyz_txt(test_structures))
 
-        if output is sys.stdout:
+        if output1 is sys.stdout:
             import matplotlib.pyplot as plt
             plt.plot(path[0], path[1], 'k', linewidth=2.0)
             plt.scatter(np.array(csm[label1_name])[filter_mask],

cosymlib/file_io/__init__.py

@@ -227,18 +227,18 @@ def get_molecule_from_file_fchk(file_name, read_multiple=False):
         else:
             Cb = []
 
-        electronic_structure = ElectronicStructure(charge=input_molecule[0][0],
-                                                   multiplicity=input_molecule[1][0],
-                                                   basis=basis,
+        electronic_structure = ElectronicStructure(basis=basis,
                                                    orbital_coefficients=[Ca, Cb],
+                                                   charge=input_molecule[0][0],
+                                                   multiplicity=input_molecule[1][0],
                                                    mo_energies=energies_alpha,
-                                                   alpha_electrons=input_molecule[2],
-                                                   beta_electrons=input_molecule[3])
+                                                   alpha_occupancy=[1]*int(input_molecule[2][0]),
+                                                   beta_occupancy=[1]*int(input_molecule[3][0]))
 
         if read_multiple:
             return [Molecule(geometry, electronic_structure)]
-        else:
-            return Molecule(geometry, electronic_structure)
+
+        return Molecule(geometry, electronic_structure)
 
 
 def get_molecule_from_file_molden(file_name, read_multiple=False):
@@ -299,7 +299,7 @@ def get_molecule_from_file_molden(file_name, read_multiple=False):
                             spin = re.split('[= ]', line)[-1]
                             input_molecule[spin + ' MO coefficients'].append([])
                         elif 'Occup' in line:
-                            occupation[spin].append(float(line.split('=')[-1]))
+                            occupation[spin].append(int(float(line.split('=')[-1])))
                         else:
                             input_molecule[spin+' MO coefficients'][-1].append(line.split()[1])
 
@@ -327,21 +327,22 @@ def get_molecule_from_file_molden(file_name, read_multiple=False):
         symbols = []
         for atom_number in input_molecule['Atomic numbers']:
             total_n_electrons += int(atom_number)
-            symbols.append(tools.atomic_number_to_element(atom_number))
+            symbols.append(atomic_number_to_element(atom_number))
 
-        occupation['Alpha'] = np.array(occupation['Alpha'], dtype=float)
-        occupation['Beta'] = np.array(occupation['Beta'], dtype=float)
+        new_occupation = []
         if len(occupation['Beta']) == 0:
-            if 2.0 not in occupation['Alpha']:
-                occupation['Alpha'] = 2*occupation['Alpha']
-            n_electrons = sum(occupation['Alpha'])
-            input_molecule['Multiplicity'].append(sum(map(lambda x: x % 2 == 1, occupation['Alpha'])) + 1)
-        else:
-            if 2.0 not in occupation['Beta']:
-                occupation['Beta'] = 2*occupation['Beta']
-            n_electrons = sum(occupation['Alpha'] + occupation['Beta'])
-            input_molecule['Multiplicity'].append(sum(occupation['Alpha'] - occupation['Beta']) + 1)
-
+            if 2 not in occupation['Alpha']:
+                occupation['Beta'] = occupation['Alpha']
+            else:
+                for occup in occupation['Alpha']:
+                    new_occupation.append(occup//2)
+                    if occup == 2:
+                        occupation['Beta'].append(1)
+                    else:
+                        occupation['Beta'].append(0)
+                occupation['Alpha'] = new_occupation
+        n_electrons = sum(occupation['Alpha'] + occupation['Beta'])
+        input_molecule['Multiplicity'].append(sum(occupation['Alpha']) - sum(occupation['Beta']) + 1)
         input_molecule['Charge'].append(int(total_n_electrons - n_electrons))
 
         basis = basis_format(basis_set_name='UNKNOWN',
@@ -364,16 +365,18 @@ def get_molecule_from_file_molden(file_name, read_multiple=False):
         else:
             Cb = []
 
-        ee = ElectronicStructure(charge=input_molecule['Charge'][0],
-                                 multiplicity=input_molecule['Multiplicity'][0],
-                                 basis=basis,
+        ee = ElectronicStructure(basis=basis,
                                  orbital_coefficients=[Ca, Cb],
-                                 mo_energies=input_molecule['MO Energies'])
+                                 charge=input_molecule['Charge'][0],
+                                 multiplicity=input_molecule['Multiplicity'][0],
+                                 mo_energies=input_molecule['MO Energies'],
+                                 alpha_occupancy=occupation['Alpha'],
+                                 beta_occupancy=occupation['Beta'])
 
         if read_multiple:
             return [Molecule(geometry, ee)]
-        else:
-            Molecule(geometry, ee)
+
+        return Molecule(geometry, ee)
 
 
 def get_geometry_from_file_ref(file_name, read_multiple=False):

cosymlib/molecule/__init__.py

@@ -2,6 +2,7 @@
 from cosymlib.molecule.electronic_structure import ElectronicStructure
 from cosymlib.simulation import ExtendedHuckel
 from warnings import warn
+from copy import deepcopy
 
 
 # Gets the parameters defined in the arguments of the function and sets them to Symmetry instance
@@ -40,6 +41,12 @@ def __init__(self, geometry,
         self._shape = geometry._shape
         self._symmetry.set_electronic_structure(electronic_structure)
 
+    def molecule_copy(self):
+        return deepcopy(self)
+
+    def set_name(self, name):
+        self._name = name
+
     @property
     def name(self):
         return self._name
@@ -69,8 +76,8 @@ def electronic_structure(self):
                                                              orbital_coefficients=[eh.get_mo_coefficients(), []],
                                                              mo_energies=eh.get_mo_energies(),
                                                              multiplicity=eh.get_multiplicity(),
-                                                             alpha_electrons=[eh.get_alpha_electrons()],
-                                                             beta_electrons=[eh.get_beta_electrons()])
+                                                             alpha_occupancy=[1]*eh.get_alpha_electrons(),
+                                                             beta_occupancy=[1]*eh.get_beta_electrons())
             self.symmetry.set_electronic_structure(self._electronic_structure)
 
         return self._electronic_structure