diff --git a/aiida_quantumespresso/calculations/pwimmigrant.py b/aiida_quantumespresso/calculations/pwimmigrant.py
index a8d3ca0e2..736a53020 100644
--- a/aiida_quantumespresso/calculations/pwimmigrant.py
+++ b/aiida_quantumespresso/calculations/pwimmigrant.py
@@ -199,7 +199,7 @@ class or using the ``set_`` methods (e.g. ``set_remote_workdir``).
             # Parse the input file.
             local_path = os.path.join(folder.abspath, self._INPUT_FILE_NAME)
             with open(local_path) as fin:
-                pwinputfile = pwinputparser.PwInputFile(fin)
+                pwinputfile = pwinputparser.PwInputFile(fin.read())
 
             # Determine PREFIX, if it hasn't already been set by the user.
             if self._PREFIX is None:
diff --git a/aiida_quantumespresso/cli/data/structure.py b/aiida_quantumespresso/cli/data/structure.py
index 90c74d7bb..a3190b866 100644
--- a/aiida_quantumespresso/cli/data/structure.py
+++ b/aiida_quantumespresso/cli/data/structure.py
@@ -21,7 +21,8 @@ def cmd_import(filename, dry_run):
     """Import a `StructureData` from a Quantum ESPRESSO input file."""
     from aiida_quantumespresso.tools.pwinputparser import PwInputFile
 
-    parser = PwInputFile(filename)
+    with open(filename, 'r') as input_file:
+        parser = PwInputFile(input_file.read())
     structure = parser.get_structuredata()
     formula = structure.get_formula()
 
diff --git a/aiida_quantumespresso/data/force_constants.py b/aiida_quantumespresso/data/force_constants.py
index 16f91a0bf..20d925af7 100644
--- a/aiida_quantumespresso/data/force_constants.py
+++ b/aiida_quantumespresso/data/force_constants.py
@@ -2,8 +2,8 @@
 """Sub class of `Data` to handle interatomic force constants produced by the Quantum ESPRESSO q2r.x code."""
 import numpy
 
-from qe_tools.constants import bohr_to_ang
 from aiida.orm import SinglefileData
+from qe_tools import CONSTANTS
 
 
 class ForceConstantsData(SinglefileData):
@@ -147,7 +147,9 @@ def parse_q2r_force_constants_file(lines, also_force_constants=False):
 
         # read cell data
         cell = tuple(
-            tuple(float(c) * celldm[0] * bohr_to_ang for c in l.split()) for l in lines[current_line:current_line + 3]
+            tuple(float(c) * celldm[0] * CONSTANTS.bohr_to_ang
+                  for c in l.split())
+            for l in lines[current_line:current_line + 3]
         )
         parsed_data['cell'] = cell
         current_line += 3
@@ -169,7 +171,7 @@ def parse_q2r_force_constants_file(lines, also_force_constants=False):
                 line[0] = atom_type_list[ityp - 1][0]  # string with element name
                 line[1] = atom_type_list[ityp - 1][1]  # element mass in amu_ry
                 # Convert atomic positions (in cartesian) from alat to Angstrom:
-                line[2:] = [pos * celldm[0] * bohr_to_ang for pos in line[2:]]
+                line[2:] = [pos * celldm[0] * CONSTANTS.bohr_to_ang for pos in line[2:]]
             atom_list.append(tuple(line))
             current_line += 1
 
diff --git a/aiida_quantumespresso/parsers/cp.py b/aiida_quantumespresso/parsers/cp.py
index 794a7e014..1057706a5 100644
--- a/aiida_quantumespresso/parsers/cp.py
+++ b/aiida_quantumespresso/parsers/cp.py
@@ -5,7 +5,7 @@
 from aiida.common import NotExistent
 from aiida.orm import Dict, TrajectoryData
 
-from qe_tools.constants import bohr_to_ang, hartree_to_ev, timeau_to_sec
+from qe_tools import CONSTANTS
 from aiida_quantumespresso.parsers.parse_raw.cp import parse_cp_raw_output, parse_cp_traj_stanzas
 from .base import Parser
 
@@ -66,9 +66,12 @@ def parse(self, **kwargs):
             return self.exit(self.exit_codes.ERROR_READING_POS_FILE)
 
         trajectories = [
-            ('positions', 'pos', bohr_to_ang, out_dict['number_of_atoms']),
-            ('cells', 'cel', bohr_to_ang, 3),
-            ('velocities', 'vel', bohr_to_ang / timeau_to_sec * 10**12, out_dict['number_of_atoms']),
+            ('positions', 'pos', CONSTANTS.bohr_to_ang, out_dict['number_of_atoms']),
+            ('cells', 'cel', CONSTANTS.bohr_to_ang, 3),
+            (
+                'velocities', 'vel', CONSTANTS.bohr_to_ang / CONSTANTS.timeau_to_sec * 10**12,
+                out_dict['number_of_atoms']
+            ),
         ]
 
         for name, extension, scale, elements in trajectories:
@@ -109,26 +112,26 @@ def parse(self, **kwargs):
                 #print "New version"
                 raw_trajectory['steps'] = numpy.array(matrix[:, 0], dtype=int)
                 raw_trajectory['evp_times'] = matrix[:, 1]  # TPS, ps
-                raw_trajectory['electronic_kinetic_energy'] = matrix[:, 2] * hartree_to_ev  # EKINC, eV
+                raw_trajectory['electronic_kinetic_energy'] = matrix[:, 2] * CONSTANTS.hartree_to_ev  # EKINC, eV
                 raw_trajectory['cell_temperature'] = matrix[:, 3]  # TEMPH, K
                 raw_trajectory['ionic_temperature'] = matrix[:, 4]  # TEMPP, K
-                raw_trajectory['scf_total_energy'] = matrix[:, 5] * hartree_to_ev  # ETOT, eV
-                raw_trajectory['enthalpy'] = matrix[:, 6] * hartree_to_ev  # ENTHAL, eV
-                raw_trajectory['enthalpy_plus_kinetic'] = matrix[:, 7] * hartree_to_ev  # ECONS, eV
-                raw_trajectory['energy_constant_motion'] = matrix[:, 8] * hartree_to_ev  # ECONT, eV
-                raw_trajectory['volume'] = matrix[:, 9] * (bohr_to_ang**3)  # volume, angstrom^3
+                raw_trajectory['scf_total_energy'] = matrix[:, 5] * CONSTANTS.hartree_to_ev  # ETOT, eV
+                raw_trajectory['enthalpy'] = matrix[:, 6] * CONSTANTS.hartree_to_ev  # ENTHAL, eV
+                raw_trajectory['enthalpy_plus_kinetic'] = matrix[:, 7] * CONSTANTS.hartree_to_ev  # ECONS, eV
+                raw_trajectory['energy_constant_motion'] = matrix[:, 8] * CONSTANTS.hartree_to_ev  # ECONT, eV
+                raw_trajectory['volume'] = matrix[:, 9] * (CONSTANTS.bohr_to_ang**3)  # volume, angstrom^3
                 raw_trajectory['pressure'] = matrix[:, 10]  # out_press, GPa
             else:
                 #print "Old version"
                 raw_trajectory['steps'] = numpy.array(matrix[:, 0], dtype=int)
-                raw_trajectory['electronic_kinetic_energy'] = matrix[:, 1] * hartree_to_ev  # EKINC, eV
+                raw_trajectory['electronic_kinetic_energy'] = matrix[:, 1] * CONSTANTS.hartree_to_ev  # EKINC, eV
                 raw_trajectory['cell_temperature'] = matrix[:, 2]  # TEMPH, K
                 raw_trajectory['ionic_temperature'] = matrix[:, 3]  # TEMPP, K
-                raw_trajectory['scf_total_energy'] = matrix[:, 4] * hartree_to_ev  # ETOT, eV
-                raw_trajectory['enthalpy'] = matrix[:, 5] * hartree_to_ev  # ENTHAL, eV
-                raw_trajectory['enthalpy_plus_kinetic'] = matrix[:, 6] * hartree_to_ev  # ECONS, eV
-                raw_trajectory['energy_constant_motion'] = matrix[:, 7] * hartree_to_ev  # ECONT, eV
-                raw_trajectory['volume'] = matrix[:, 8] * (bohr_to_ang**3)  # volume, angstrom^3
+                raw_trajectory['scf_total_energy'] = matrix[:, 4] * CONSTANTS.hartree_to_ev  # ETOT, eV
+                raw_trajectory['enthalpy'] = matrix[:, 5] * CONSTANTS.hartree_to_ev  # ENTHAL, eV
+                raw_trajectory['enthalpy_plus_kinetic'] = matrix[:, 6] * CONSTANTS.hartree_to_ev  # ECONS, eV
+                raw_trajectory['energy_constant_motion'] = matrix[:, 7] * CONSTANTS.hartree_to_ev  # ECONT, eV
+                raw_trajectory['volume'] = matrix[:, 8] * (CONSTANTS.bohr_to_ang**3)  # volume, angstrom^3
                 raw_trajectory['pressure'] = matrix[:, 9]  # out_press, GPa
                 raw_trajectory['evp_times'] = matrix[:, 10]  # TPS, ps
 
diff --git a/aiida_quantumespresso/parsers/matdyn.py b/aiida_quantumespresso/parsers/matdyn.py
index 2d53549dc..7a9132a8e 100644
--- a/aiida_quantumespresso/parsers/matdyn.py
+++ b/aiida_quantumespresso/parsers/matdyn.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 from aiida import orm
 from aiida.common import exceptions
-from qe_tools.constants import invcm_to_THz
+from qe_tools import CONSTANTS
 
 from aiida_quantumespresso.calculations.matdyn import MatdynCalculation
 from .base import Parser
@@ -114,7 +114,7 @@ def parse_raw_matdyn_phonon_file(phonon_frequencies):
     for i in range(num_kpoints):
         for j in range(num_bands):
             try:
-                freq_matrix[i, j] = corrected_data[counter] * invcm_to_THz  # from cm-1 to THz
+                freq_matrix[i, j] = corrected_data[counter] * CONSTANTS.invcm_to_THz  # from cm-1 to THz
             except ValueError:
                 parsed_data['warnings'].append('Error while parsing the frequencies')
             except IndexError:
diff --git a/aiida_quantumespresso/parsers/parse_raw/neb.py b/aiida_quantumespresso/parsers/parse_raw/neb.py
index 5b0d31293..fbfb13304 100644
--- a/aiida_quantumespresso/parsers/parse_raw/neb.py
+++ b/aiida_quantumespresso/parsers/parse_raw/neb.py
@@ -5,7 +5,7 @@
 specific functionalities. The parsing will try to convert whatever it can in some dictionary, which by operative
 decision doesn't have much structure encoded, [the values are simple ]
 """
-from qe_tools.constants import bohr_to_ang
+from qe_tools import CONSTANTS
 
 from aiida_quantumespresso.parsers import QEOutputParsingError, get_parser_info
 from aiida_quantumespresso.parsers.parse_raw import convert_qe_time_to_sec
@@ -148,10 +148,10 @@ def parse_neb_text_output(data, input_dict={}):
     for count, line in enumerate(lines):
         if 'initial path length' in line:
             initial_path_length = float(line.split('=')[1].split('bohr')[0])
-            parsed_data['initial_path_length'] = initial_path_length * bohr_to_ang
+            parsed_data['initial_path_length'] = initial_path_length * CONSTANTS.bohr_to_ang
         elif 'initial inter-image distance' in line:
             initial_image_dist = float(line.split('=')[1].split('bohr')[0])
-            parsed_data['initial_image_dist'] = initial_image_dist * bohr_to_ang
+            parsed_data['initial_image_dist'] = initial_image_dist * CONSTANTS.bohr_to_ang
         elif 'string_method' in line:
             parsed_data['string_method'] = line.split('=')[1].strip()
         elif 'restart_mode' in line:
@@ -235,9 +235,9 @@ def parse_neb_text_output(data, input_dict={}):
                 iteration_data['climbing_image_auto'].append([int(_) for _ in line.split('=')[1].split(',')])
             elif 'path length' in line:
                 path_length = float(line.split('=')[1].split('bohr')[0])
-                iteration_data['path_length'].append(path_length * bohr_to_ang)
+                iteration_data['path_length'].append(path_length * CONSTANTS.bohr_to_ang)
             elif 'inter-image distance' in line:
                 image_dist = float(line.split('=')[1].split('bohr')[0])
-                iteration_data['image_dist'].append(image_dist * bohr_to_ang)
+                iteration_data['image_dist'].append(image_dist * CONSTANTS.bohr_to_ang)
 
     return parsed_data, dict(iteration_data), list(critical_warnings.values())
diff --git a/aiida_quantumespresso/parsers/parse_raw/ph.py b/aiida_quantumespresso/parsers/parse_raw/ph.py
index 2117e18e6..4e5f61124 100644
--- a/aiida_quantumespresso/parsers/parse_raw/ph.py
+++ b/aiida_quantumespresso/parsers/parse_raw/ph.py
@@ -6,7 +6,7 @@
 """
 import numpy
 
-from qe_tools.constants import *
+from qe_tools import CONSTANTS
 
 from aiida_quantumespresso.parsers import QEOutputParsingError, get_parser_info
 from aiida_quantumespresso.parsers.parse_raw.base import convert_qe_time_to_sec
@@ -282,7 +282,7 @@ def parse_ph_dynmat(data, logs, lattice_parameter=None, also_eigenvectors=False,
                 header_dict['ibrav'] = int(pieces[2])
                 header_dict['celldm'] = [float(i) for i in pieces[3:]]
                 # In angstrom
-                alat = header_dict['celldm'][0] * bohr_to_ang
+                alat = header_dict['celldm'][0] * CONSTANTS.bohr_to_ang
                 if abs(alat) < 1.e-5:
                     raise QEOutputParsingError(
                         'Lattice constant=0! Probably you are using an '
@@ -318,7 +318,7 @@ def parse_ph_dynmat(data, logs, lattice_parameter=None, also_eigenvectors=False,
                 try:
                     if int(pieces[0]) != idx:
                         raise QEOutputParsingError('Error with the indices of the species')
-                    species.append([pieces[1].strip(), float(pieces[2]) / amu_Ry])
+                    species.append([pieces[1].strip(), float(pieces[2]) / CONSTANTS.amu_Ry])
                 except ValueError:
                     raise QEOutputParsingError('Error parsing the species')
 
diff --git a/aiida_quantumespresso/parsers/parse_raw/pw.py b/aiida_quantumespresso/parsers/parse_raw/pw.py
index 75605a440..7c3e72811 100644
--- a/aiida_quantumespresso/parsers/parse_raw/pw.py
+++ b/aiida_quantumespresso/parsers/parse_raw/pw.py
@@ -7,7 +7,7 @@
 """
 import re
 import numpy
-from qe_tools.constants import ry_to_ev, bohr_to_ang, ry_si, bohr_si
+from qe_tools import CONSTANTS
 
 from aiida_quantumespresso.parsers import QEOutputParsingError
 from aiida_quantumespresso.parsers.parse_raw import convert_qe_time_to_sec
@@ -334,8 +334,8 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
                 elif 'Smooth grid' in line:
                     smooth_FFT_grid = [int(g) for g in line.split('(')[1].split(')')[0].split(',')]
                     break
-            alat *= bohr_to_ang
-            volume *= bohr_to_ang**3
+            alat *= CONSTANTS.bohr_to_ang
+            volume *= CONSTANTS.bohr_to_ang**3
             parsed_data['lattice_parameter_initial'] = alat
             parsed_data['number_of_bands'] = nbnd
             try:
@@ -514,19 +514,19 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
                         )
 
                     if 'alat' in lattice[0].lower():
-                        a1 = [alat * bohr_to_ang * float(s) for s in a1]
-                        a2 = [alat * bohr_to_ang * float(s) for s in a2]
-                        a3 = [alat * bohr_to_ang * float(s) for s in a3]
+                        a1 = [alat * CONSTANTS.bohr_to_ang * float(s) for s in a1]
+                        a2 = [alat * CONSTANTS.bohr_to_ang * float(s) for s in a2]
+                        a3 = [alat * CONSTANTS.bohr_to_ang * float(s) for s in a3]
                         lattice_parameter_b = float(lattice[1])
                         if abs(lattice_parameter_b - alat) > lattice_tolerance:
                             raise QEOutputParsingError(
                                 'Lattice parameters mismatch! ' + '{} vs {}'.format(lattice_parameter_b, alat)
                             )
                     elif 'bohr' in lattice[0].lower():
-                        lattice_parameter_b *= bohr_to_ang
-                        a1 = [bohr_to_ang * float(s) for s in a1]
-                        a2 = [bohr_to_ang * float(s) for s in a2]
-                        a3 = [bohr_to_ang * float(s) for s in a3]
+                        lattice_parameter_b *= CONSTANTS.bohr_to_ang
+                        a1 = [CONSTANTS.bohr_to_ang * float(s) for s in a1]
+                        a2 = [CONSTANTS.bohr_to_ang * float(s) for s in a2]
+                        a3 = [CONSTANTS.bohr_to_ang * float(s) for s in a3]
                     trajectory_data.setdefault('lattice_vectors_relax', []).append([a1, a2, a3])
 
                 except Exception:
@@ -550,7 +550,7 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
                         if metric == 'alat':
                             tau = [alat * float(s) for s in tau]
                         elif metric == 'bohr':
-                            tau = [bohr_to_ang * float(s) for s in tau]
+                            tau = [CONSTANTS.bohr_to_ang * float(s) for s in tau]
                         positions.append(tau)
                     trajectory_data.setdefault(this_key, []).append(positions)
                 except Exception:
@@ -590,7 +590,7 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
             # If for some step this line is not printed, the later check with the scf_accuracy array length should catch it
             elif 'estimated scf accuracy' in line:
                 try:
-                    value = float(line.split()[-2]) * ry_to_ev
+                    value = float(line.split()[-2]) * CONSTANTS.ry_to_ev
                     trajectory_data.setdefault('scf_accuracy', []).append(value)
                 except Exception:
                     logs.warning.append('Error while parsing scf accuracy.')
@@ -654,8 +654,8 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
             elif '!' in line:
                 try:
 
-                    En = float(line.split('=')[1].split('Ry')[0]) * ry_to_ev
-                    E_acc = float(data_step[count + 2].split('<')[1].split('Ry')[0]) * ry_to_ev
+                    En = float(line.split('=')[1].split('Ry')[0]) * CONSTANTS.ry_to_ev
+                    E_acc = float(data_step[count + 2].split('<')[1].split('Ry')[0]) * CONSTANTS.ry_to_ev
 
                     for key, value in [['energy', En], ['energy_accuracy', E_acc]]:
                         trajectory_data.setdefault(key, []).append(value)
@@ -737,7 +737,7 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
                         if 'atom ' in line2:
                             line2 = line2.split('=')[1].split()
                             # CONVERT FORCES IN eV/Ang
-                            vec = [float(s) * ry_to_ev / bohr_to_ang for s in line2]
+                            vec = [float(s) * CONSTANTS.ry_to_ev / CONSTANTS.bohr_to_ang for s in line2]
                             forces.append(vec)
                         if len(forces) == nat:
                             break
@@ -750,7 +750,7 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
 
             elif 'Total force =' in line:
                 try:  # note that I can't check the units: not written in output!
-                    value = float(line.split('=')[1].split('Total')[0]) * ry_to_ev / bohr_to_ang
+                    value = float(line.split('=')[1].split('Total')[0]) * CONSTANTS.ry_to_ev / CONSTANTS.bohr_to_ang
                     trajectory_data.setdefault('total_force', []).append(value)
                     parsed_data['total_force' + units_suffix] = default_force_units
                 except Exception:
@@ -767,7 +767,7 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
                         raise QEOutputParsingError('Error while parsing stress: unexpected units.')
                     for k in range(3):
                         line2 = data_step[count2 + k + 1].split()
-                        vec = [float(s) * 10**(-9) * ry_si / (bohr_si)**3 for s in line2[0:3]]
+                        vec = [float(s) * 10**(-9) * CONSTANTS.ry_si / (CONSTANTS.bohr_si)**3 for s in line2[0:3]]
                         stress.append(vec)
                     trajectory_data.setdefault('stress', []).append(stress)
                     parsed_data['stress' + units_suffix] = default_stress_units
@@ -885,6 +885,6 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
 
 def grep_energy_from_line(line):
     try:
-        return float(line.split('=')[1].split('Ry')[0]) * ry_to_ev
+        return float(line.split('=')[1].split('Ry')[0]) * CONSTANTS.ry_to_ev
     except Exception:
         raise QEOutputParsingError('Error while parsing energy')
diff --git a/aiida_quantumespresso/parsers/parse_xml/pw/legacy.py b/aiida_quantumespresso/parsers/parse_xml/pw/legacy.py
index d6146ca7e..e4e03c263 100644
--- a/aiida_quantumespresso/parsers/parse_xml/pw/legacy.py
+++ b/aiida_quantumespresso/parsers/parse_xml/pw/legacy.py
@@ -6,7 +6,7 @@
 
 from aiida_quantumespresso.parsers import QEOutputParsingError
 from aiida_quantumespresso.utils.mapping import get_logging_container
-from qe_tools.constants import ry_to_ev, hartree_to_ev, bohr_to_ang
+from qe_tools import CONSTANTS
 
 units_suffix = '_units'
 default_energy_units = 'eV'
@@ -205,16 +205,16 @@ def parse_pw_xml_pre_6_2(xml_file, dir_with_bands):
     if parsed_data.get('two_fermi_energies', False):
         tagname = 'FERMI_ENERGY_UP'
         parsed_data[tagname.replace('-','_').lower()] = \
-            parse_xml_child_float(tagname,target_tags) * hartree_to_ev
+            parse_xml_child_float(tagname,target_tags) * CONSTANTS.hartree_to_ev
         parsed_data[tagname.lower() + units_suffix] = default_energy_units
         tagname = 'FERMI_ENERGY_DOWN'
         parsed_data[tagname.replace('-','_').lower()] = \
-            parse_xml_child_float(tagname,target_tags) * hartree_to_ev
+            parse_xml_child_float(tagname,target_tags) * CONSTANTS.hartree_to_ev
         parsed_data[tagname.lower() + units_suffix] = default_energy_units
     else:
         tagname = 'FERMI_ENERGY'
         parsed_data[tagname.replace('-','_').lower()] = \
-            parse_xml_child_float(tagname,target_tags) * hartree_to_ev
+            parse_xml_child_float(tagname,target_tags) * CONSTANTS.hartree_to_ev
         parsed_data[tagname.lower() + units_suffix] = default_energy_units
 
     #CARD MAGNETIZATION_INIT
@@ -303,7 +303,7 @@ def read_bands_and_occupations(eigenval_n):
                     tagname = 'EIGENVALUES'
                     a = eig_dom.getElementsByTagName(tagname)[0]
                     b = a.childNodes[0]
-                    value_e = [float(s) * hartree_to_ev for s in b.data.split()]
+                    value_e = [float(s) * CONSTANTS.hartree_to_ev for s in b.data.split()]
 
                     tagname = 'OCCUPATIONS'
                     a = eig_dom.getElementsByTagName(tagname)[0]
@@ -520,7 +520,7 @@ def xml_card_cell(parsed_data, dom):
             )
         )
     if metric == 'bohr':
-        value *= bohr_to_ang
+        value *= CONSTANTS.bohr_to_ang
     parsed_data[tagname.replace('-', '_').lower()] = value
 
     tagname = 'CELL_DIMENSIONS'
@@ -558,7 +558,7 @@ def xml_card_cell(parsed_data, dom):
             d = c.data.replace('\n', '').split()
             value = [float(i) for i in d]
             if metric == 'bohr':
-                value = [bohr_to_ang * float(s) for s in value]
+                value = [CONSTANTS.bohr_to_ang * float(s) for s in value]
             lattice_vectors.append(value)
 
         volume = cell_volume(lattice_vectors[0], lattice_vectors[1], lattice_vectors[2])
@@ -682,7 +682,7 @@ def xml_card_ions(parsed_data, dom, lattice_vectors, volume):
             if metric == 'alat':
                 tau = [parsed_data['lattice_parameter_xml'] * float(s) for s in tau]
             elif metric == 'bohr':
-                tau = [bohr_to_ang * float(s) for s in tau]
+                tau = [CONSTANTS.bohr_to_ang * float(s) for s in tau]
             atomlist.append([chem_symbol, tau])
             tagname2 = 'if_pos'
             b = a.getAttribute(tagname2)
@@ -745,9 +745,9 @@ def xml_card_planewaves(parsed_data, dom, calctype):
             raise QEOutputParsingError('Units {} are not supported by parser'.format(units))
     else:
         if 'hartree' in units:
-            conv_fac = hartree_to_ev
+            conv_fac = CONSTANTS.hartree_to_ev
         else:
-            conv_fac = ry_to_ev
+            conv_fac = CONSTANTS.ry_to_ev
 
         tagname = 'WFC_CUTOFF'
         parsed_data[tagname.lower()] = parse_xml_child_float(tagname, target_tags) * conv_fac
@@ -882,7 +882,7 @@ def xml_card_exchangecorrelation(parsed_data, dom):
                 a = [_ for _ in target_tags.childNodes if _.nodeName == tagname][0]
                 b = a.childNodes[0]
                 c = b.data.replace('\n', ' ').split()  # note the need of a white space!
-                value = [float(i) * ry_to_ev for i in c]
+                value = [float(i) * CONSTANTS.ry_to_ev for i in c]
                 parsed_data[tagname.lower()] = value
             except Exception:
                 raise QEOutputParsingError('Error parsing tag '+\
diff --git a/aiida_quantumespresso/parsers/parse_xml/pw/parse.py b/aiida_quantumespresso/parsers/parse_xml/pw/parse.py
index bbfa93201..c5bfdd5c3 100644
--- a/aiida_quantumespresso/parsers/parse_xml/pw/parse.py
+++ b/aiida_quantumespresso/parsers/parse_xml/pw/parse.py
@@ -5,7 +5,7 @@
 
 from xmlschema import XMLSchema
 from xmlschema.etree import ElementTree
-from qe_tools.constants import hartree_to_ev, bohr_to_ang, ry_to_ev
+from qe_tools import CONSTANTS
 
 from aiida_quantumespresso.utils.mapping import get_logging_container
 
@@ -106,9 +106,9 @@ def parse_pw_xml_post_6_2(xml):
     outputs = xml_dictionary['output']
 
     lattice_vectors = [
-        [x * bohr_to_ang for x in outputs['atomic_structure']['cell']['a1']],
-        [x * bohr_to_ang for x in outputs['atomic_structure']['cell']['a2']],
-        [x * bohr_to_ang for x in outputs['atomic_structure']['cell']['a3']],
+        [x * CONSTANTS.bohr_to_ang for x in outputs['atomic_structure']['cell']['a1']],
+        [x * CONSTANTS.bohr_to_ang for x in outputs['atomic_structure']['cell']['a2']],
+        [x * CONSTANTS.bohr_to_ang for x in outputs['atomic_structure']['cell']['a3']],
     ]
 
     has_electric_field = inputs.get('electric_field', {}).get('electric_potential', None) == 'sawtooth_potential'
@@ -262,8 +262,8 @@ def parse_pw_xml_post_6_2(xml):
         inversion_symmetry,  # the old tag was INVERSION_SYMMETRY and was set to (from the code): "invsym    if true the system has inversion symmetry"
         'number_of_bravais_symmetries': nrot,  # lattice symmetries
         'number_of_symmetries': nsym,  # crystal symmetries
-        'wfc_cutoff': inputs['basis']['ecutwfc'] * hartree_to_ev,
-        'rho_cutoff': outputs['basis_set']['ecutrho'] * hartree_to_ev,  # not always printed in input->basis
+        'wfc_cutoff': inputs['basis']['ecutwfc'] * CONSTANTS.hartree_to_ev,
+        'rho_cutoff': outputs['basis_set']['ecutrho'] * CONSTANTS.hartree_to_ev,  # not always printed in input->basis
         'smooth_fft_grid': [value for _, value in sorted(outputs['basis_set']['fft_smooth'].items())],
         'dft_exchange_correlation': inputs['dft']['functional'],  # TODO: also parse optional elements of 'dft' tag
         # WARNING: this is different between old XML and new XML
@@ -274,7 +274,7 @@ def parse_pw_xml_post_6_2(xml):
     # alat is technically an optional attribute according to the schema,
     # but I don't know what to do if it's missing. atomic_structure is mandatory.
     output_alat_bohr = outputs['atomic_structure']['@alat']
-    output_alat_angstrom = output_alat_bohr * bohr_to_ang
+    output_alat_angstrom = output_alat_bohr * CONSTANTS.bohr_to_ang
 
     # Band structure
     if 'band_structure' in outputs:
@@ -292,9 +292,9 @@ def parse_pw_xml_post_6_2(xml):
 
             # Versions below 19.03.04 (Quantum ESPRESSO<=6.4.1) incorrectly print degauss in Ry instead of Hartree
             if xml_version < StrictVersion('19.03.04'):
-                degauss *= ry_to_ev
+                degauss *= CONSTANTS.ry_to_ev
             else:
-                degauss *= hartree_to_ev
+                degauss *= CONSTANTS.hartree_to_ev
 
             xml_data['degauss'] = degauss
             xml_data['smearing_type'] = smearing_xml['$']
@@ -357,7 +357,7 @@ def parse_pw_xml_post_6_2(xml):
                 band_occupations[0].append(ks_state['occupations']['$'][0:num_bands_up])
                 band_occupations[1].append(ks_state['occupations']['$'][num_bands_up:num_bands])
 
-        band_eigenvalues = np.array(band_eigenvalues) * hartree_to_ev
+        band_eigenvalues = np.array(band_eigenvalues) * CONSTANTS.hartree_to_ev
         band_occupations = np.array(band_occupations)
 
         if not spins:
@@ -387,7 +387,7 @@ def parse_pw_xml_post_6_2(xml):
                 xml_data[key] = value
 
         if 'fermi_energy' in band_structure:
-            xml_data['fermi_energy'] = band_structure['fermi_energy'] * hartree_to_ev
+            xml_data['fermi_energy'] = band_structure['fermi_energy'] * CONSTANTS.hartree_to_ev
 
         bands_dict = {
             'occupations': band_occupations,
@@ -497,7 +497,7 @@ def parse_pw_xml_post_6_2(xml):
         # - individual electronic phases and weights
 
     # TODO: We should put the `non_periodic_cell_correction` string in (?)
-    atoms = [[atom['@name'], [coord * bohr_to_ang
+    atoms = [[atom['@name'], [coord * CONSTANTS.bohr_to_ang
                               for coord in atom['$']]]
              for atom in outputs['atomic_structure']['atomic_positions']['atom']]
     species = outputs['atomic_species']['species']
diff --git a/aiida_quantumespresso/tools/base.py b/aiida_quantumespresso/tools/base.py
index bdd36c20d..9a7900598 100644
--- a/aiida_quantumespresso/tools/base.py
+++ b/aiida_quantumespresso/tools/base.py
@@ -42,7 +42,7 @@ def get_structuredata(self):
         """, re.X | re.I
         )
 
-        data = self.get_structure_from_qeinput()
+        data = self.structure
         species = self.atomic_species
 
         structure = StructureData()
diff --git a/aiida_quantumespresso/tools/cpinputparser.py b/aiida_quantumespresso/tools/cpinputparser.py
index 4d08672c6..f52cad87f 100644
--- a/aiida_quantumespresso/tools/cpinputparser.py
+++ b/aiida_quantumespresso/tools/cpinputparser.py
@@ -4,7 +4,7 @@
 from .base import StructureParseMixin
 
 
-class CpInputFile(StructureParseMixin, BaseCpInputFile):
+class CpInputFile(StructureParseMixin, BaseCpInputFile):  # pylint: disable=too-few-public-methods
     """Parser of Quantum ESPRESSO cp.x input file into AiiDA nodes.
 
     .. note:: This mixes in :class:`~aiida_quantumespresso.tools.base.StructureParseMixin` which adds the functionality
diff --git a/aiida_quantumespresso/tools/dbexporters/tcod_plugins/pw.py b/aiida_quantumespresso/tools/dbexporters/tcod_plugins/pw.py
index 55da93159..7826ff814 100644
--- a/aiida_quantumespresso/tools/dbexporters/tcod_plugins/pw.py
+++ b/aiida_quantumespresso/tools/dbexporters/tcod_plugins/pw.py
@@ -260,7 +260,7 @@ def get_integration_Methfessel_Paxton_order(cls, calc, **kwargs):  # pylint: dis
     @classmethod
     def get_kinetic_energy_cutoff_wavefunctions(cls, calc, **kwargs):  # pylint: disable=invalid-name,unused-argument
         """Return kinetic energy cutoff for wavefunctions in eV."""
-        from qe_tools.constants import ry_to_ev
+        from qe_tools import CONSTANTS
         parameters = calc.inputs.parameters
         ecutwfc = None
         try:
@@ -269,7 +269,7 @@ def get_kinetic_energy_cutoff_wavefunctions(cls, calc, **kwargs):  # pylint: dis
             pass
         if ecutwfc is None:
             return None
-        return ecutwfc * ry_to_ev
+        return ecutwfc * CONSTANTS.ry_to_ev
 
     @classmethod
     def get_kinetic_energy_cutoff_charge_density(cls, calc, **kwargs):  # pylint: disable=invalid-name,unused-argument
@@ -278,10 +278,10 @@ def get_kinetic_energy_cutoff_charge_density(cls, calc, **kwargs):  # pylint: di
         .. note :: by default returns 4 * ecutwfc, as indicated in
             http://www.quantum-espresso.org/wp-content/uploads/Doc/INPUT_PW.html
         """
-        from qe_tools.constants import ry_to_ev
+        from qe_tools import CONSTANTS
         parameters = calc.inputs.parameters
         try:
-            return parameters.get_dict()['SYSTEM']['ecutrho'] * ry_to_ev
+            return parameters.get_dict()['SYSTEM']['ecutrho'] * CONSTANTS.ry_to_ev
         except KeyError:
             pass
         ecutwfc = cls.get_kinetic_energy_cutoff_wavefunctions(calc)
@@ -296,10 +296,10 @@ def get_kinetic_energy_cutoff_EEX(cls, calc, **kwargs):  # pylint: disable=inval
         .. note :: by default returns ecutrho, as indicated in
             http://www.quantum-espresso.org/wp-content/uploads/Doc/INPUT_PW.html
         """
-        from qe_tools.constants import ry_to_ev
+        from qe_tools import CONSTANTS
         parameters = calc.inputs.parameters
         try:
-            return parameters.get_dict()['SYSTEM']['ecutfock'] * ry_to_ev
+            return parameters.get_dict()['SYSTEM']['ecutfock'] * CONSTANTS.ry_to_ev
         except KeyError:
             pass
         return cls.get_kinetic_energy_cutoff_charge_density(calc)
diff --git a/aiida_quantumespresso/tools/pwinputparser.py b/aiida_quantumespresso/tools/pwinputparser.py
index 26d4a0b3f..30cb23fd3 100644
--- a/aiida_quantumespresso/tools/pwinputparser.py
+++ b/aiida_quantumespresso/tools/pwinputparser.py
@@ -91,7 +91,7 @@ def create_builder_from_file(input_folder, input_file_name, code, metadata, pseu
         input_folder = Folder(input_folder)
 
     with input_folder.open(input_file_name) as input_file:
-        parsed_file = PwInputFile(input_file)
+        parsed_file = PwInputFile(input_file.read())
 
     builder.structure = parsed_file.get_structuredata()
     builder.kpoints = parsed_file.get_kpointsdata()
diff --git a/aiida_quantumespresso/utils/validation/trajectory.py b/aiida_quantumespresso/utils/validation/trajectory.py
index 95b6d6212..607b3c06a 100644
--- a/aiida_quantumespresso/utils/validation/trajectory.py
+++ b/aiida_quantumespresso/utils/validation/trajectory.py
@@ -46,12 +46,12 @@ def verify_convergence_forces(trajectory, index=-1, threshold=None):
     :return: `True` if threshold is valid, `False` otherwise
     :raises ValueError: if the `forces` array or given index does not exist
     """
-    from qe_tools.constants import ry_to_ev, bohr_to_ang
+    from qe_tools import CONSTANTS
 
     if threshold is None:
         return None
 
-    threshold *= ry_to_ev / bohr_to_ang  # Convert to eV / Å
+    threshold *= CONSTANTS.ry_to_ev / CONSTANTS.bohr_to_ang  # Convert to eV / Å
 
     try:
         forces = trajectory.get_array('forces')[index]
diff --git a/setup.json b/setup.json
index 8b7edf048..25fa8eec7 100644
--- a/setup.json
+++ b/setup.json
@@ -87,7 +87,7 @@
     "install_requires": [
         "aiida_core[atomic_tools]~=1.2",
         "packaging",
-        "qe-tools~=1.1",
+        "qe-tools~=2.0rc1",
         "xmlschema~=1.2"
     ],
     "license": "MIT License",
diff --git a/tests/conftest.py b/tests/conftest.py
index fa3f3376f..36cc3a02f 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -139,11 +139,12 @@ def _generate_calc_job_node(
             node.set_attribute_many(attributes)
 
         if filepath_folder:
-            from qe_tools.utils.exceptions import ParsingError
+            from qe_tools.exceptions import ParsingError
             from aiida_quantumespresso.tools.pwinputparser import PwInputFile
             try:
-                parsed_input = PwInputFile(filepath_input)
-            except ParsingError:
+                with open(filepath_input, 'r') as input_file:
+                    parsed_input = PwInputFile(input_file.read())
+            except (ParsingError, FileNotFoundError):
                 pass
             else:
                 inputs['structure'] = parsed_input.get_structuredata()