commit 1d8c374e394c53254524ef7e0ec1e7ee002141bf
Author: Carles Martí <carles.marti2@gmail.com>
Date:   Fri Dec 4 14:08:58 2020 +0100

    Simplified the reading of coordinates and energies and the selection of conformers.

diff --git a/modules/calculation.py b/modules/calculation.py
index 50d77fb..9be7269 100644
--- a/modules/calculation.py
+++ b/modules/calculation.py
@@ -17,17 +17,15 @@ def check_finished_calcs(run_type, code):
     """
     from glob import glob
     import ase.io
-    from pycp2k import CP2K
-    from modules.utilities import tail
+    from modules.utilities import tail, is_binary
 
     finished_calcs = []
     unfinished_calcs = []
     for conf_dir in sorted(os.listdir(run_type), key=_human_key):
-        conf_path = f'{run_type}/{conf_dir}'
+        conf_path = f'{run_type}/{conf_dir}/'
         if not os.path.isdir(conf_path) or 'conf_' not in conf_dir:
             continue
         if code == 'cp2k':
-            cp2k = CP2K()
             restart_file_list = glob(f"{conf_path}/*-1.restart")
             if len(restart_file_list) == 0:
                 logger.warning(f"No *-1.restart file found on {conf_path}.")
@@ -42,7 +40,9 @@ def check_finished_calcs(run_type, code):
                 continue
             out_files = []
             for file in os.listdir(conf_path):
-                with open(conf_path+"/"+file, "rb") as out_fh:
+                if is_binary(conf_path+file):
+                    continue
+                with open(conf_path+file, "rb") as out_fh:
                     tail_out_str = tail(out_fh)
                 if tail_out_str.count("PROGRAM STOPPED IN") == 1:
                     out_files.append(file)
diff --git a/modules/formats.py b/modules/formats.py
index 1d9ecde..c652a45 100644
--- a/modules/formats.py
+++ b/modules/formats.py
@@ -166,7 +166,8 @@ def read_coords_cp2k(file, spec_atoms=tuple()):
     force_eval = cp2k.CP2K_INPUT.FORCE_EVAL_list[0]
     raw_coords = force_eval.SUBSYS.COORD.Default_keyword
     symbols = [atom.split()[0] for atom in raw_coords]
-    positions = np.array([atom.split()[1:] for atom in raw_coords])
+    positions = np.array([[float(coord) for coord in atom.split()[1:]]
+                          for atom in raw_coords])
     if len(spec_atoms) > 0:
         add_special_atoms(spec_atoms)
     return Atoms(symbols=symbols, positions=positions)
@@ -190,45 +191,11 @@ def read_coords_vasp(file, spec_atoms=tuple()):
     return ase.io.read(file)
 
 
-def collect_coords(conf_list, code, run_type, spec_atms=tuple()):
-    """Directs the reading of coordinates on a set of subdirectories.
-
-    Given a dockonsurf directory hierarchy: project/run_type/conf_X
-    (run_type = ['isolated', 'screening' or 'refinement']) with finished
-    calculations produced by a given code, stored in every conf_X subdirectory,
-    it collects the coordinates of the specified conf_X subdirectories in a
-    single run_type by calling the adequate function (depending on the code) and
-    returns a list of ase.Atoms objects.
-
-    @param conf_list: List of directories where to read the coords from.
-    @param code: the code that produced the calculation results files.
-    @param run_type: the type of calculation (and also the name of the folder)
-                     containing the calculation subdirectories.
-    @param spec_atms: List of tuples containing pairs of new/traditional
-        chemical symbols.
-    @return: list of ase.Atoms objects.
-    """
-    from glob import glob
-    atoms_list = []
-    for conf in conf_list:
-        if code == 'cp2k':
-            atoms_list.append(read_coords_cp2k(glob(f"{run_type}/{conf}/*-1"
-                                                    f".restart")[0], spec_atms))
-        elif code == 'vasp':
-            atoms_list.append(read_coords_vasp(f"{run_type}/{conf}/OUTCAR",
-                                               spec_atms))
-        else:
-            err_msg = f"Collect coordinates not implemented for '{code}'."
-            logger.error(err_msg)
-            raise NotImplementedError(err_msg)
-    return atoms_list
-
-
 def read_energy_cp2k(file):
-    """Reads the CP2K out file and returns its final energy.
+    """Reads the CP2K output file and returns its final energy.
 
     @param file: The file from which the energy should be read.
-    @return: The last energy on the out file.
+    @return: The last energy on the out file in eV.
     """
     out_fh = open(file, 'r')
     energy = None
@@ -239,38 +206,50 @@ def read_energy_cp2k(file):
     return energy
 
 
-def collect_energies(conf_list, code, run_type):
-    """Directs the reading of energies on a set of subdirectories.
+def collect_confs(dir_list, code, run_type, spec_atms=tuple()):
+    """Reads the coordinates and energies of a list of finished calculations.
 
     Given a dockonsurf directory hierarchy: project/run_type/conf_X
-    (run_type = ['isolated', 'screening' or 'refinement']) with finished
-    calculations produced by a given code, stored in every conf_X subdirectory,
-    it collects the energies of the specified conf_X subdirectories in a
-    single run_type by calling the adequate function (depending on the code) and
-    returns a list of energies.
-
-    @param conf_list: List of directories where to read the energy.
-    @param code: The code that produced the calculation output files.
-    @param run_type: The type of calculation (and also the name of the folder)
-                     containing the calculation subdirectories.
-    @return: List of energies
+    (run_type = ['isolated', 'screening' or 'refinement']) it reads the
+    coordinates of each conf_X, it assigns its total energy from the calculation
+    and assigns the conf_X label to track its origin. Finally it returns the
+    ase.Atoms object.
+
+    @param dir_list: List of directories where to read the coords from.
+    @param code: the code that produced the calculation results files.
+    @param run_type: the type of calculation (and also the name of the folder)
+        containing the calculation subdirectories.
+    @param spec_atms: List of tuples containing pairs of new/traditional
+        chemical symbols.
+    @return: list of ase.Atoms objects.
     """
     from glob import glob
-    import numpy as np
-
-    energies = []
-    for conf in conf_list:
+    import os
+    from modules.utilities import is_binary
+    atoms_list = []
+    for conf_dir in dir_list:
+        conf_path = f"{run_type}/{conf_dir}/"
         if code == 'cp2k':
-            energies.append(read_energy_cp2k(
-                glob(f"{run_type}/{conf}/*.out")[0]))
+            ase_atms = read_coords_cp2k(glob(f"{conf_path}/*-1.restart")[0],
+                                        spec_atms)
+            # Assign energy
+            for fil in os.listdir(conf_path):
+                if is_binary(conf_path+fil):
+                    continue
+                conf_energy = read_energy_cp2k(conf_path+fil)
+                if conf_energy is not None:
+                    ase_atms.info["energy"] = conf_energy
+                    break
+            ase_atms.info[run_type[:3]] = conf_dir
+            atoms_list.append(ase_atms)
+        elif code == 'vasp':
+            ase_atms = read_coords_vasp(f"{conf_path}/OUTCAR", spec_atms)
+            ase_atms.info["energy"] = ase_atms.get_total_energy() * 27.2113845
+            ase_atms.info[run_type[:3]] = conf_dir
+            atoms_list.append(ase_atms)
         else:
-            err_msg = f"Collect energies not implemented for '{code}'."
+            err_msg = f"Collect coordinates not implemented for '{code}'."
             logger.error(err_msg)
             raise NotImplementedError(err_msg)
+    return atoms_list
 
-    if len(energies) == 0:
-        err = f"No results found on {run_type}"
-        logger.error(err)
-        raise FileNotFoundError(err)
-
-    return np.array(energies)
diff --git a/modules/isolated.py b/modules/isolated.py
index 95bf1e5..db37d1a 100644
--- a/modules/isolated.py
+++ b/modules/isolated.py
@@ -138,7 +138,7 @@ def run_isolated(inp_vars):
     @return:
     """
     from modules.formats import adapt_format, confs_to_mol_list, \
-        rdkit_mol_to_ase_atoms, collect_confs, collect_energies
+        rdkit_mol_to_ase_atoms, collect_confs
     from modules.clustering import clustering, get_rmsd
     from modules.calculation import run_calc, check_finished_calcs
     from modules.refinement import select_stable_confs
@@ -163,11 +163,11 @@ def run_isolated(inp_vars):
         logger.error(err_msg)
         raise ValueError(err_msg)
     run_calc('isolated', inp_vars, ase_atms_list)
-    finished_calcs, unfinished_calcs = check_finished_calcs('isolated',
+    finished_calcs, failed_calcs = check_finished_calcs('isolated',
                                                             inp_vars['code'])
-    conf_list = collect_coords(finished_calcs, inp_vars['code'], 'isolated',
-                               inp_vars['special_atoms'])
-    most_stable_conf = select_confs(conf_list, finished_calcs, 0,
-                                    inp_vars['code'])[0]
-    logger.info('Finished the procedures for the isolated molecule section. '
-                f'Most stable conformers is {most_stable_conf}.')
+    conf_list = collect_confs(finished_calcs, inp_vars['code'], 'isolated',
+                              inp_vars['special_atoms'])
+    most_stable_conf = select_stable_confs(conf_list, 0)[0]
+    logger.info("Finished the procedures for the isolated molecule section. "
+                f"Most stable conformers is {most_stable_conf.info['id']}, "
+                f"with a total energy of {most_stable_conf.info['energy']} eV.")
diff --git a/modules/refinement.py b/modules/refinement.py
index cabf787..e218738 100644
--- a/modules/refinement.py
+++ b/modules/refinement.py
@@ -3,7 +3,7 @@ import logging
 logger = logging.getLogger('DockOnSurf')
 
 
-def select_confs(conf_list, calc_dirs, energy_cutoff, code):
+def select_stable_confs(conf_list, energy_cutoff):
     """From a list of atomic configurations selects the most stable ones.
 
     Given a list of ase.Atoms configurations and an energy cutoff, selects all
@@ -11,20 +11,10 @@ def select_confs(conf_list, calc_dirs, energy_cutoff, code):
     conformer plus the cutoff.
 
     @param conf_list: List of ase.Atoms objects of the conformers
-    @param calc_dirs: List of the directories of the finished calculations
-    @param energy_cutoff: The maximum energy above the most stable configuration
-    @param code: the code used to carry out the screening of structures.
-    @return: list of the the most stable configurations within the energy cutoff
+    @param energy_cutoff: The maximum energy above the most stable
+        configuration.
+    @return: list of the the most stable configurations within the energy cutoff.
     """
-    from copy import deepcopy
-    from modules.formats import collect_energies
-
-    conf_enrgs = collect_energies(calc_dirs, code, 'screening')
-
-    for i, conf in enumerate(conf_list):
-        conf.info['energy'] = conf_enrgs[i]
-        conf.info['id'] = calc_dirs[i]
-
     sorted_list = sorted(conf_list, key=lambda conf: conf.info['energy'])
     lowest_e = sorted_list[0].info['energy']
     return [conf for conf in sorted_list
@@ -37,7 +27,7 @@ def run_refinement(inp_vars):
     @param inp_vars: Calculation parameters from input file.
     """
     import os
-    from modules.formats import collect_coords
+    from modules.formats import collect_confs
     from modules.calculation import run_calc, check_finished_calcs
 
     logger.info('Carrying out procedures for the refinement of '
@@ -49,8 +39,8 @@ def run_refinement(inp_vars):
         logger.error(err)
         raise FileNotFoundError(err)
 
-    finished_calcs, unfinished_calcs = check_finished_calcs('screening',
-                                                            inp_vars['code'])
+    finished_calcs, failed_calcs = check_finished_calcs('screening',
+                                                        inp_vars['code'])
     if not finished_calcs:
         err_msg = "No calculations on 'screening' finished normally."
         logger.error(err_msg)
@@ -58,25 +48,24 @@ def run_refinement(inp_vars):
     logger.info(f"Found {len(finished_calcs)} structures of "
                 f"adsorbate-surface atomic configurations whose calculation"
                 f" finished normally.")
-    if len(unfinished_calcs) != 0:
-        logger.warning(f"Found {len(unfinished_calcs)} calculations more that "
-                       f"did not finish normally: {unfinished_calcs}. \n"
+    if len(failed_calcs) != 0:
+        logger.warning(f"Found {len(failed_calcs)} calculations more that "
+                       f"did not finish normally: {failed_calcs}. \n"
                        f"Using only the ones that finished normally: "
                        f"{finished_calcs}.")
 
-    conf_list = collect_coords(finished_calcs, inp_vars['code'], 'screening',
-                               inp_vars['special_atoms'])
-    selected_confs = select_confs(conf_list, finished_calcs,
-                                  inp_vars['energy_cutoff'], inp_vars['code'])
+    conf_list = collect_confs(finished_calcs, inp_vars['code'], 'screening',
+                              inp_vars['special_atoms'])
+    selected_confs = select_stable_confs(conf_list, inp_vars['energy_cutoff'])
     logger.info(f"Selected {len(selected_confs)} structures to carry out the"
                 f" refinement")
     run_calc('refinement', inp_vars, selected_confs)
-    finished_calcs, unfinished_calcs = check_finished_calcs('refinement',
-                                                            inp_vars['code'])
-    conf_list = collect_coords(finished_calcs, inp_vars['code'], 'refinement',
-                               inp_vars['special_atoms'])
-    most_stable_conf = select_confs(conf_list, finished_calcs, 0,
-                                    inp_vars['code'])[0]
-    logger.info('Finished the procedures for the refinement of '
-                'adsorbate-surface structures section. Most stable structure '
-                f"is {most_stable_conf.info['id']}")
+    finished_calcs, failed_calcs = check_finished_calcs('refinement',
+                                                        inp_vars['code'])
+    conf_list = collect_confs(finished_calcs, inp_vars['code'], 'refinement',
+                              inp_vars['special_atoms'])
+    most_stable_conf = select_stable_confs(conf_list, 0)[0]
+    logger.info("Finished the procedures for the refinement of "
+                "adsorbate-surface structures section. Most stable structure "
+                f"is {most_stable_conf.info['id']} with a Total energy of "
+                f"{most_stable_conf.info['energy']} eV.")
diff --git a/modules/screening.py b/modules/screening.py
index c268208..14ba2c6 100644
--- a/modules/screening.py
+++ b/modules/screening.py
@@ -5,70 +5,44 @@ import ase
 logger = logging.getLogger('DockOnSurf')
 
 
-def assign_prop(atoms: ase.Atoms, prop_name: str, prop_val):  # TODO Needed?
-    atoms.info[prop_name] = prop_val
-
-
-def select_confs(orig_conf_list: list, calc_dirs: list, magns: list,
-                 num_sel: int, code: str):
+def select_confs(conf_list: list, magns: list, num_sel: int):
     """Takes a list ase.Atoms and selects the most different magnitude-wise.
 
     Given a list of ase.Atoms objects and a list of magnitudes, it selects a
     number of the most different conformers according to every magnitude
     specified.
      
-    @param orig_conf_list: list of ase.Atoms objects to select among.
-    @param calc_dirs: List of directories where to read the energies from.
+    @param conf_list: list of ase.Atoms objects to select among.
     @param magns: list of str with the names of the magnitudes to use for the
         conformer selection. Supported magnitudes: 'energy', 'moi'.
     @param num_sel: number of conformers to select for every of the magnitudes.
-    @param code: The code that generated the magnitude information.
-         Supported codes: See formats.py
     @return: list of the selected ase.Atoms objects.
     """
-    from copy import deepcopy
-    from modules.formats import collect_energies
-
-    conf_list = deepcopy(orig_conf_list)
-    conf_enrgs, mois, selected_ids = [], [], []
+    selected_ids = []
     if num_sel >= len(conf_list):
         logger.warning('Number of conformers per magnitude is equal or larger '
                        'than the total number of conformers. Using all '
                        f'available conformers: {len(conf_list)}.')
         return conf_list
 
-    # Read properties
-    if 'energy' in magns:
-        if code == 'cp2k':
-            conf_enrgs = collect_energies(calc_dirs, code, 'isolated')
-        elif code == 'vasp':
-            conf_enrgs = np.array([conf.get_total_energy()
-                                   for conf in orig_conf_list])
+    # Assign mois
     if 'moi' in magns:
-        mois = np.array([conf.get_moments_of_inertia() for conf in conf_list])
-
-    # Assign values
-    for i, conf in enumerate(conf_list):
-        assign_prop(conf, 'idx', i)
-        assign_prop(conf, 'iso', calc_dirs[i])
-        if 'energy' in magns:
-            assign_prop(conf, 'energy', conf_enrgs[i])
-        if 'moi' in magns:
-            assign_prop(conf, 'moi', mois[i, 2])
+        for conf in conf_list:
+            conf.info["moi"] = conf.get_moments_of_inertia()[2]
 
     # pick ids
     for magn in magns:
         sorted_list = sorted(conf_list, key=lambda conf: abs(conf.info[magn]))
-        if sorted_list[-1].info['idx'] not in selected_ids:
-            selected_ids.append(sorted_list[-1].info['idx'])
+        if sorted_list[-1].info['iso'] not in selected_ids:
+            selected_ids.append(sorted_list[-1].info['iso'])
         if num_sel > 1:
             for i in range(0, len(sorted_list) - 1,
                            len(conf_list) // (num_sel - 1)):
-                if sorted_list[i].info['idx'] not in selected_ids:
-                    selected_ids.append(sorted_list[i].info['idx'])
+                if sorted_list[i].info['iso'] not in selected_ids:
+                    selected_ids.append(sorted_list[i].info['iso'])
 
     logger.info(f'Selected {len(selected_ids)} conformers for adsorption.')
-    return [conf_list[idx] for idx in selected_ids]
+    return [conf for conf in conf_list if conf.info["iso"] in selected_ids]
 
 
 def get_vect_angle(v1: list, v2: list, ref=None, degrees=True):
@@ -1026,7 +1000,7 @@ def run_screening(inp_vars):
     """
     import os
     import random
-    from modules.formats import collect_coords, adapt_format
+    from modules.formats import collect_confs, adapt_format
     from modules.calculation import run_calc, check_finished_calcs
 
     logger.info('Carrying out procedures for the screening of adsorbate-surface'
@@ -1042,29 +1016,26 @@ def run_screening(inp_vars):
             logger.error(err)
             raise FileNotFoundError(err)
 
-        correct_calcs, failed_calcs = check_finished_calcs('isolated',
-                                                           inp_vars['code'])
-        if not correct_calcs:
+        finished_calcs, failed_calcs = check_finished_calcs('isolated',
+                                                            inp_vars['code'])
+        if not finished_calcs:
             err_msg = "No calculations on 'isolated' finished normally."
             logger.error(err_msg)
             raise FileNotFoundError(err_msg)
 
-        logger.info(f"Found {len(correct_calcs)} structures of isolated "
+        logger.info(f"Found {len(finished_calcs)} structures of isolated "
                     f"conformers whose calculation finished normally.")
         if len(failed_calcs) != 0:
             logger.warning(
                 f"Found {len(failed_calcs)} calculations more that "
                 f"did not finish normally: {failed_calcs}. \n"
                 f"Using only the ones that finished normally: "
-                f"{correct_calcs}.")
-
-        conformer_atoms_list = collect_coords(correct_calcs, inp_vars['code'],
-                                              'isolated',
-                                              inp_vars['special_atoms'])
-        selected_confs = select_confs(conformer_atoms_list, correct_calcs,
-                                      inp_vars['select_magns'],
-                                      inp_vars['confs_per_magn'],
-                                      inp_vars['code'])
+                f"{finished_calcs}.")
+
+        conf_list = collect_confs(finished_calcs, inp_vars['code'], 'isolated',
+                                  inp_vars['special_atoms'])
+        selected_confs = select_confs(conf_list, inp_vars['select_magns'],
+                                      inp_vars['confs_per_magn'])
     surf = adapt_format('ase', inp_vars['surf_file'], inp_vars['special_atoms'])
     surf.info = {}
     surf_ads_list = adsorb_confs(selected_confs, surf, inp_vars)
