Merge remote-tracking branch 'origin/master'

# Conflicts:
#	cosymlib/shape/tools.py

cosymlib/__init__.py

@@ -6,7 +6,7 @@
 from cosymlib import tools
 from cosymlib.utils import swap_vectors
 from cosymlib.utils import get_shape_path, plot_molecular_orbital_diagram, plot_symmetry_energy_evolution
-from cosymlib.shape.tools import get_structure_references
+from cosymlib.shape.tools import get_structure_references, get_sym_from_label
 import matplotlib.pyplot as plt
 
 import sys
@@ -23,6 +23,38 @@ def _get_symgroup_arguments(locals):
     return kwargs
 
 
+def get_table_format(labels, molecules_names, data):
+
+    txt = 'Structure'
+    max_len_name = 12
+    max_name = 0
+    names = []
+    for name in molecules_names:
+        if len(name) > max_len_name:
+            max_name = 12
+            names.append(name[:(max_len_name - 1)])
+        elif len(name) > max_name:
+            max_name = len(name)
+            names.append(name)
+
+    n = max_len_name - 3
+    for label in labels:
+        n += len(label)
+        txt += '{}'.format(label.rjust(n))
+        n = 11 - len(label)
+    txt += '\n\n'
+
+    for idx, name in enumerate(names):
+        txt += '{:{width}}'.format(name + ',', width=max_len_name)
+        n = 11
+        for idn, label in enumerate(labels):
+            txt += '{:>{width}.{prec}f},'.format(data[idn][idx], width=n, prec=3)
+            n = 10
+        txt += '\n'
+    txt += '\n'
+    return txt
+
+
 class Cosymlib:
     """
     Main cosymlib class
@@ -152,31 +184,7 @@ def print_shape_measure(self, shape_reference, central_atom=0, fix_permutation=F
             else:
                 references_names.append(reference)
 
-        txt_shape = '{}'.format('Structure')
-        max_name = len(max(molecules_names, key=len))
-        if max_name < 9:
-            n = 5
-        else:
-            n = max_name - 4
-        for label in references_names:
-            n += len(label)
-            txt_shape += '{}'.format(label.rjust(n))
-            n = 12 - len(label)
-        txt_shape += '\n\n'
-
-        for idx, name in enumerate(molecules_names):
-            max_name = len(max(molecules_names, key=len))
-            txt_shape += '{}'.format(name)
-            if max_name < 9:
-                n = 18 - len(name)
-            else:
-                n = 9 + max_name - len(name)
-            for idn, label in enumerate(references_names):
-                txt_shape += ', {:{width}.{prec}f}'.format(measure_list[idn][idx], width=n, prec=3)
-                n = 11
-            txt_shape += '\n'
-        txt_shape += '\n'
-
+        txt_shape = get_table_format(references_names, molecules_names, measure_list)
         output.write(txt_shape)
 
     def print_shape_structure(self, shape_reference, central_atom=0, fix_permutation=False, output=sys.stdout):
@@ -220,9 +228,11 @@ def print_shape_structure(self, shape_reference, central_atom=0, fix_permutation
         for idm, molecule in enumerate(self._molecules):
             geometries.append(molecule.geometry)
             for idl, reference in enumerate(references):
-                 shape_results_structures[idl][idm].set_name(molecule.name + '_' + reference)
-                 geometries.append(shape_results_structures[idl][idm])
+                shape_results_structures[idl][idm].set_name(molecule.name + ' ' + reference + ' ' +
+                                                            get_sym_from_label(reference))
+                geometries.append(shape_results_structures[idl][idm])
 
+        print("\nOriginal structures vs reference polyhedra in file {}\n".format(output.name))
         for geometry in geometries:
             output.write(file_io.get_file_xyz_txt(geometry))
 
@@ -327,6 +337,18 @@ def print_symmetry_nearest_structure(self, label, central_atom=0, center=None, o
             geometry = molecule.geometry.get_symmetry_nearest_structure(**kwargs)
             output.write(file_io.get_file_xyz_txt(geometry))
 
+    def print_chirality_measure(self, labels, central_atom=0, center=None, output=sys.stdout):
+
+        csm_list = []
+        for label in labels:
+            csm_list.append([geometry.get_symmetry_measure(label, central_atom=central_atom, center=center)
+                             for geometry in self.get_geometries()])
+
+        molecules_names = [molecule.name for molecule in self._molecules]
+        txt = get_table_format(labels, molecules_names, csm_list)
+
+        output.write(txt)
+
     # This should be substituted by calling methods within this class
     def OLD_print_wnfsym_measure_verbose(self, group, axis=None, axis2=None, center=None, output=sys.stdout):
         self.print_wnfsym_sym_matrices(group, axis=axis, axis2=axis2, center=center, output=output)
@@ -644,7 +666,7 @@ def get_shape_measure(self, label, kind, central_atom=0, fix_permutation=False):
         """
         get_measure = 'get_shape_' + kind
 
-        return [geometry.get_shape_measure(label, central_atom=central_atom, fix_permutation=fix_permutation)
+        return [getattr(geometry, get_measure)(label, central_atom=central_atom, fix_permutation=fix_permutation)
                 for geometry in self.get_geometries()]
 
     def get_molecule_path_deviation(self, shape_label1, shape_label2, central_atom=0):

cosymlib/file_io/__init__.py

@@ -3,7 +3,7 @@
 from cosymlib.molecule.electronic_structure import ElectronicStructure
 from cosymlib.file_io.tools import extract_geometries
 from cosymlib.tools import atomic_number_to_element
-from cosymlib.file_io import custom_errors
+from cosymlib.file_io import errors
 
 import numpy as np
 import os, re
@@ -22,14 +22,12 @@ def _non_blank_lines(f):
 
 def check_geometries_vertices(geometries, file_name):
     n_atoms = geometries[0].get_n_atoms()
-    idl = 0
     for idg, geometry in enumerate(geometries):
         if geometry.get_n_atoms() != n_atoms:
-            warnings.warn('Error in structure {}: Line {} of file {}\n '
-                          'Number of vertices does not match with other structures\n'.format(idg + 1, idl + 1,
+            warnings.warn('Structure {} with name {} of file {}\n '
+                          'Number of vertices does not match with other structures\n'.format(idg + 1, geometry.name,
                                                                                              file_name),
-                          custom_errors.DifferentVerteciesWarning)
-        idl += geometry.get_n_atoms() + 2
+                          errors.DifferentVerteciesWarning)
 
 
 # Read INPUT files
@@ -62,13 +60,13 @@ def get_geometry_from_file_xyz(file_name, read_multiple=False):
     input_molecule = [[], []]
     geometries = []
     n_atoms = None
+    n_atoms_line = 1
     with open(file_name, mode='r') as lines:
         for idl, line in enumerate(lines):
-            # if '$' in line or '#' in line or '!' in line:
             if line.lstrip().startswith(comment_line):
                 pass
             elif line.strip() == '':
-                warnings.warn('Line {} is empty'.format(idl + 1), custom_errors.EmptyLineWarning)
+                warnings.warn('Line {} is empty'.format(idl + 1), errors.EmptyLineWarning)
             else:
                 try:
                     float(line.split()[1])
@@ -77,9 +75,9 @@ def get_geometry_from_file_xyz(file_name, read_multiple=False):
                 except IndexError:
                     if input_molecule[0]:
                         if len(input_molecule[0]) != n_atoms:
-                            warnings.warn('Number of atoms around line {} and number '
-                                          'of atoms provided are not equal'.format(idl - len(input_molecule[0]) - 1),
-                                          custom_errors.MissingLineWarning)
+                            warnings.warn('Number of atoms in line {} and number '
+                                          'of atoms provided are not equal'.format(n_atoms_line),
+                                          errors.MissingAtomWarning)
                         geometries.append(Geometry(symbols=input_molecule[0],
                                                    positions=input_molecule[1],
                                                    name=name))
@@ -88,13 +86,14 @@ def get_geometry_from_file_xyz(file_name, read_multiple=False):
                     input_molecule = [[], []]
                     if n_atoms is None:
                         n_atoms = int(line.split()[0])
+                        n_atoms_line = idl + 1
                     else:
                         name = line.split()[0]
 
         if len(input_molecule[0]) != n_atoms:
             warnings.warn('Number of atoms in line {} and number '
-                          'of atoms provided are not equal'.format(idl - len(input_molecule[0]) - 1),
-                          custom_errors.MissingLineWarning)
+                          'of atoms provided are not equal'.format(n_atoms_line),
+                          errors.MissingAtomWarning)
 
         geometries.append(Geometry(symbols=input_molecule[0],
                                    positions=input_molecule[1],
@@ -123,7 +122,7 @@ def get_geometry_from_file_cor(file_name, read_multiple=False):
             if line.lstrip().startswith(comment_line):
                 pass
             elif line.strip() == '':
-                warnings.warn('Line {} is empty'.format(idl + 1), custom_errors.EmptyLineWarning)
+                warnings.warn('Line {} is empty'.format(idl + 1), errors.EmptyLineWarning)
             else:
                 try:
                     float(line.split()[1])

cosymlib/file_io/classic_inputs.py

@@ -1,5 +1,5 @@
 from cosymlib.file_io import get_geometry_from_file_cor
-from cosymlib.file_io import custom_errors
+from cosymlib.file_io import errors
 
 import os
 import tempfile
@@ -55,7 +55,7 @@ def read_old_input(file_name):
                 idl += 1
                 for line in tmp_lines:
                     if line.strip() == '':
-                        warnings.warn('Line {} is empty'.format(idl + 1), custom_errors.EmptyLineWarning)
+                        warnings.warn('Line {} is empty'.format(idl + 1), errors.EmptyLineWarning)
                     else:
                         tmp.write(line)
                     idl += 1

cosymlib/file_io/custom_errors.py → cosymlib/file_io/errors.py

@@ -8,7 +8,7 @@ def warning_on_one_line(message, category, filename, lineno, file=None, line=Non
 warnings.formatwarning = warning_on_one_line
 
 
-class MissingLineWarning(Warning):
+class MissingAtomWarning(Warning):
     pass
 
 

cosymlib/shape/tools.py

@@ -5,6 +5,22 @@
 ideal_structures = None
 
 
+def distance(p1, p2):
+    return np.sqrt((p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2)
+
+
+def resize_structure(structure, center=(0, 0, 0), resize_distance=2.5):
+    resized_structure = []
+    max_distance = max([distance(coordinate, center) for coordinate in structure])
+    if max_distance >= resize_distance:
+        resized_structure = structure
+    else:
+        for coordinate in structure:
+            new_coordinate = resize_distance*(np.array(coordinate))/max_distance
+            resized_structure.append(list(new_coordinate))
+    return resized_structure
+
+
 def get_test_structure(label, central_atom=0):
     global ideal_structures
     from cosymlib.molecule.geometry import Geometry
@@ -14,22 +30,21 @@ def get_test_structure(label, central_atom=0):
         with open(file_path, 'r') as stream:
             ideal_structures = yaml.load(stream, Loader=yaml.FullLoader)
 
-    try:
-        if central_atom == 0:
-            coordinates = ideal_structures[label][:-1]
-            return Geometry(positions=coordinates,
-                            name=label,
-                            symbols=['L'] * len(coordinates),
-                            connectivity=[])
+    if central_atom == 0:
+        coordinates = ideal_structures[label][:-1]
+        # coordinates = resize_structure(ideal_structures[label][:-1])
+        return Geometry(positions=coordinates,
+                        name=label,
+                        symbols=['L'] * len(coordinates),
+                        connectivity=[])
 
-        else:
-            coordinates = ideal_structures[label]
-            return Geometry(positions=coordinates,
-                            name=label,
-                            symbols=['L'] * (len(coordinates)-1) + ['M'],
-                            connectivity=[])
-    except KeyError as e:
-        raise Exception('Shape reference label "{}" not found!'.format(e.args[0]))
+    else:
+        coordinates = ideal_structures[label]
+        # coordinates = resize_structure(ideal_structures[label], center=ideal_structures[label][-1])
+        return Geometry(positions=coordinates,
+                        name=label,
+                        symbols=['L'] * (len(coordinates)-1) + ['M'],
+                        connectivity=[])
 
 
 def get_structure_references(vertices):
@@ -55,6 +70,13 @@ def get_shape_label(code, vertices):
             return label[0]
 
 
+def get_sym_from_label(label):
+    vertices = label[-1]
+    for labels in shape_structure_references['{} Vertices'.format(vertices)]:
+        if labels[0] == label:
+            return labels[2]
+
+
 def get_shape_label_info(vertices, old=False, with_central_atom=False):
     if with_central_atom:
         vertices = vertices-1
@@ -161,4 +183,4 @@ def get_shape_label_info(vertices, old=False, with_central_atom=False):
                               '24 Vertices': [['TCU-24', 1, 'Oh', 'Truncated cube'],
                                               ['TOC-24', 2, 'Oh', 'Truncated octahedron']],
                               '48 Vertices': [['TCOC-48', 1, 'Oh', 'Truncated cuboctahedron']],
-                              '60 Vertices': [['TRIC-60', 1, 'Ih', 'Truncated icosahedron (fullerene)']]}
+                              '60 Vertices': [['TRIC-60', 1, 'Ih', 'Truncated icosahedron (fullerene)']]}

cosymlib/symmetry/tools.py

@@ -9,6 +9,7 @@ def print_symmetry_labels():
     for k, label in operations.items():
         print('{:6} {}'.format(k, label))
 
+
 def get_group_num_from_label(label):
 
     operations = {'Cn':  '',

examples/coord.xyz

@@ -12,7 +12,7 @@
    H      0.2051      0.8240     -0.6786
    H      0.3345     -0.9314     -0.4496
    H     -1.0685     -0.0537      0.1921
-    5
+ 5
  Ammonium2
    N      0.0000      0.0000      0.1000
    H      0.5288      0.1610      0.9359

scripts/chirality

@@ -23,8 +23,9 @@ parser.add_argument(type=str,
 parser.add_argument('-m', '--measure',
                     dest='measure',
                     metavar='SG',
-                    default=False,
-                    help='compute the SG symmetry measure of the input structures')
+                    default=0,
+                    type=int,
+                    help='compute the chiral symmetry measure of the input structures')
 parser.add_argument('-o', '--output',
                     dest='output_name',
                     metavar='filename',
@@ -76,33 +77,36 @@ if args.version:
 common_output = open(args.output_name, 'w') if args.output_name is not None else sys.stdout
 print_header(common_output)
 
-try:
-    structures = read_generic_structure_file(args.input_file, read_multiple=True)
-except Exception as e:
-    sys.exit('No input file selected! An existing file must be provide')
-
+if args.input_file is None:
+    parser.error('No input file selected! An existing file must be provide')
 
+structures = read_generic_structure_file(args.input_file, read_multiple=True)
 structure_set = Cosymlib(structures)
 
 if args.info:
     print_input_info(structure_set.get_geometries(), output=common_output)
     exit()
 
-# if not args.measure:
-#     reference = ['ci', 'cs']
-# elif args.measure not in ['ci', 'cs']:
-#     reference = ['ci', 'cs', args.measure]
-# else:
-#     reference = args.measure
-if not args.measure:
-    sys.exit('User should provide a symmetry group measure')
-reference = args.measure.lower()
-
-# Symgroup commands
+if args.measure == 0:
+    reference = ['Ci', 'Cs']
+elif args.measure == 1:
+    reference = ['Cs']
+elif args.measure == 2:
+    reference = ['Ci']
+elif int(args.measure) % 2 != 0:
+    raise ValueError('For Sn chirality measure choose an even value between 2-n')
+else:
+    for n in range(4, args.measure + 1):
+        if n%2 == 0:
+            reference = ['Ci', 'Cs', 'S'+str(n)]
+
+reference = [r.lower() for r in reference]
+
+# Chirality commands
 if args.measure:
-    structure_set.print_geometric_symmetry_measure(reference,
-                                                   central_atom=args.central_atom,
-                                                   center=args.center,
-                                                   output=common_output)
+    structure_set.print_chirality_measure(reference,
+                                          central_atom=args.central_atom,
+                                          center=args.center,
+                                          output=common_output)
 
 print_footer(common_output)