diff --git a/aiida_quantumespresso/calculations/namelists.py b/aiida_quantumespresso/calculations/namelists.py
index 75aa23ce2..be927e263 100644
--- a/aiida_quantumespresso/calculations/namelists.py
+++ b/aiida_quantumespresso/calculations/namelists.py
@@ -13,7 +13,6 @@
 from aiida.engine import CalcJob
 from aiida.orm import Dict
 from aiida.orm import RemoteData, FolderData, SinglefileData
-from aiida.plugins import ParserFactory
 
 from aiida_quantumespresso.calculations import _lowercase_dict, _uppercase_dict, _pop_parser_options
 from aiida_quantumespresso.utils.convert import convert_input_to_namelist_entry
@@ -39,7 +38,8 @@ class NamelistsCalculation(CalcJob):
     _default_namelists = ['INPUTPP']
     _blocked_keywords = []  # a list of tuples with key and value fixed
 
-    _retrieve_singlefile_list = []
+    _retrieve_temporary_list = ()
+    _retrieve_output_as_temp = False
 
     _DEFAULT_INPUT_FILE = 'aiida.in'
     _DEFAULT_OUTPUT_FILE = 'aiida.out'
@@ -174,12 +174,16 @@ def prepare_for_submission(self, folder):
 
         # Retrieve by default the output file and the xml file
         calcinfo.retrieve_list = []
-        calcinfo.retrieve_list.append(self.inputs.metadata.options.output_filename)
+        if not self._retrieve_output_as_temp:
+            calcinfo.retrieve_list.append(self.inputs.metadata.options.output_filename)
         settings_retrieve_list = settings.pop('ADDITIONAL_RETRIEVE_LIST', [])
         calcinfo.retrieve_list += settings_retrieve_list
         calcinfo.retrieve_list += self._internal_retrieve_list
 
-        calcinfo.retrieve_singlefile_list = self._retrieve_singlefile_list
+        calcinfo.retrieve_temporary_list = []
+        if self._retrieve_output_as_temp:
+            calcinfo.retrieve_temporary_list.append(self.inputs.metadata.options.output_filename)
+        calcinfo.retrieve_temporary_list += list(self._retrieve_temporary_list)
 
         # We might still have parser options in the settings dictionary: pop them.
         _pop_parser_options(self, settings)
diff --git a/aiida_quantumespresso/calculations/projwfc.py b/aiida_quantumespresso/calculations/projwfc.py
index 184cea752..3b3240b9f 100644
--- a/aiida_quantumespresso/calculations/projwfc.py
+++ b/aiida_quantumespresso/calculations/projwfc.py
@@ -3,6 +3,7 @@
 from aiida.orm import RemoteData, FolderData, Dict, XyData
 from aiida_quantumespresso.calculations.namelists import NamelistsCalculation
 
+
 class ProjwfcCalculation(NamelistsCalculation):
     """
     Projwfc.x code of the Quantum ESPRESSO distribution, handles the the
@@ -24,6 +25,7 @@ class ProjwfcCalculation(NamelistsCalculation):
     ]
     _default_parser = 'quantumespresso.projwfc'
     _internal_retrieve_list = [NamelistsCalculation._PREFIX + '.pdos*']
+    _retrieve_output_as_temp = True
 
     @classmethod
     def define(cls, spec):
@@ -31,23 +33,33 @@ def define(cls, spec):
         super(ProjwfcCalculation, cls).define(spec)
         spec.input('parent_folder', valid_type=(RemoteData, FolderData), help='The output folder of a pw.x calculation')
         spec.output('output_parameters', valid_type=Dict)
+        spec.output('lowdin', valid_type=Dict)
         spec.output('Dos', valid_type=XyData)
         # if spin
-        spec.output('projections_up',   valid_type=ProjectionData, required=False)
+        spec.output('projections_up', valid_type=ProjectionData, required=False)
         spec.output('projections_down', valid_type=ProjectionData, required=False)
-        spec.output('bands_up',   valid_type=BandsData, required=False)
+        spec.output('bands_up', valid_type=BandsData, required=False)
         spec.output('bands_down', valid_type=BandsData, required=False)
         # if non-spin
         spec.output('projections', valid_type=ProjectionData, required=False)
         spec.output('bands', valid_type=BandsData, required=False)
         spec.default_output_node = 'output_parameters'
         spec.exit_code(
-            100, 'ERROR_NO_RETRIEVED_FOLDER', message='The retrieved folder data node could not be accessed.')
+            100, 'ERROR_NO_RETRIEVED_FOLDER', message='The retrieved folder data node could not be accessed.'
+        )
+        spec.exit_code(
+            101, 'ERROR_NO_RETRIEVED_TEMPORARY_FOLDER', message='The retrieved temporary folder could not be accessed.'
+        )
         spec.exit_code(
-            110, 'ERROR_READING_OUTPUT_FILE', message='The output file could not be read from the retrieved folder.')
+            110, 'ERROR_READING_OUTPUT_FILE', message='The output file could not be read from the retrieved folder.'
+        )
         spec.exit_code(
-            111, 'ERROR_READING_PDOSTOT_FILE', message='The pdos_tot file could not be read from the retrieved folder.')
+            111, 'ERROR_READING_PDOSTOT_FILE', message='The pdos_tot file could not be read from the retrieved folder.'
+        )
         spec.exit_code(
-            112, 'ERROR_PARSING_PROJECTIONS', message='An exception was raised parsing bands and projections.')
+            112, 'ERROR_PARSING_PROJECTIONS', message='An exception was raised parsing bands and projections.'
+        )
+        spec.exit_code(113, 'ERROR_PARSING_LOWDIN', message='An exception was raised parsing Lowdin charges.')
         spec.exit_code(
-            130, 'ERROR_JOB_NOT_DONE', message='The computation did not finish properly (\'JOB DONE\' not found).')
+            130, 'ERROR_JOB_NOT_DONE', message='The computation did not finish properly (\'JOB DONE\' not found).'
+        )
diff --git a/aiida_quantumespresso/parsers/parse_raw/projwfc.py b/aiida_quantumespresso/parsers/parse_raw/projwfc.py
new file mode 100644
index 000000000..3041bfe7c
--- /dev/null
+++ b/aiida_quantumespresso/parsers/parse_raw/projwfc.py
@@ -0,0 +1,58 @@
+from __future__ import absolute_import
+import re
+
+
+def parse_lowdin_charges(out_file_lines, lowdin_lines=None):
+    """Parse the Lowdin Charge section of the ``projwfc.x`` output.
+
+    :param out_file_lines: The projwfc.x stdout file content
+    :type out_file_lines: list[str]
+    :param lowdin_lines: indices of lines where 'Lowdin Charges:' has been found
+    :return: A tuple of lowdin data dict and spill parameter float
+    """
+    # find start of Lowdin charge output
+    if lowdin_lines is None:
+        lowdin_lines = []
+        for i, line in enumerate(out_file_lines):
+            if line.strip() == 'Lowdin Charges:':
+                lowdin_lines.append(i)
+
+    if len(lowdin_lines) > 1:
+        raise IOError("'Lowdin Charges:' found on multiple lines: {}".format(lowdin_lines))
+
+    if not lowdin_lines:
+        return None, None
+
+    spill_parameter = None
+    started_section = False
+    atom_index = None
+    output = {}
+    for i, line in enumerate(out_file_lines[lowdin_lines[0] + 1:]):
+        lineno = lowdin_lines[0] + 2 + i
+        # break on empty line or spill parameter
+        if not line.strip():
+            if started_section:
+                break
+            started_section = True
+            continue
+        if line.strip().startswith('Spilling Parameter'):
+            spill_parameter = float(line.strip().split()[-1])
+            break
+        atom_line_match = re.match(r'^\s*Atom\s\#\s*(\d+)\:(.*)$', line)
+        if atom_line_match:
+            atom_index = int(atom_line_match.group(1))
+            line = atom_line_match.group(2)
+        if atom_index is None:
+            raise IOError('No atom index specified on or before line {}'.format(lineno))
+        charge_type_match = re.match(r'^\s*(total charge|spin up|spin down|polarization)\s*\=\s*(\-?[\d\.]+)', line)
+        if charge_type_match:
+            charge_type = charge_type_match.group(1).replace(' ', '_')
+            output.setdefault(atom_index, {}).setdefault('sum', {})[charge_type] = float(charge_type_match.group(2))
+        else:
+            raise IOError('No charge type specified on line {}'.format(lineno))
+        for orbital_type, value in re.findall(
+            r'(s|p|d|f|px|py|pz|dxy|dxz|dyz|dz2|dx2-y2|f[xyz23\-\+]+)\s*\=\s*(\-?[\d\.]+)', line
+        ):
+            output.setdefault(atom_index, {}).setdefault(charge_type, {})[orbital_type] = float(value)
+
+    return output, spill_parameter
diff --git a/aiida_quantumespresso/parsers/projwfc.py b/aiida_quantumespresso/parsers/projwfc.py
index e9f8ec841..64744c804 100644
--- a/aiida_quantumespresso/parsers/projwfc.py
+++ b/aiida_quantumespresso/parsers/projwfc.py
@@ -1,7 +1,8 @@
 # -*- coding: utf-8 -*-
 from __future__ import absolute_import, division
-import re
+import os
 import fnmatch
+import re
 import traceback
 
 import numpy as np
@@ -14,24 +15,25 @@
 
 from aiida_quantumespresso.parsers import QEOutputParsingError
 from aiida_quantumespresso.parsers import parse_raw_out_basic
+from aiida_quantumespresso.parsers.parse_raw.projwfc import parse_lowdin_charges
 
 
-def find_orbitals_from_statelines(out_info_dict):
-    """
-    This function reads in all the state_lines, that is, the lines describing
-    which atomic states, taken from the pseudopotential, are used for the
-    projection. Then it converts these state_lines into a set of orbitals.
+def find_orbitals_from_statelines(out_file_lines, statelines_idx, positions):
+    """Read all state-lines and convert into a set of orbitals.
+
+    State-lines are lines, taken from the pseudopotential,
+    describing which atomic states are used for the projection.
 
-    :param out_info_dict: contains various technical internals useful in parsing
+    :param out_file_lines: content of the projwfc stdout file
+    :param statelines_idx: out_file_lines indices containing state-lines
+    :param positions: list of structure site positions
     :return: orbitals, a list of orbitals suitable for setting ProjectionData
     """
-    out_file = out_info_dict['out_file']
     atomnum_re = re.compile(r'atom (.*?)\(')
     element_re = re.compile(r'\((.*?)\)')
     lnum_re = re.compile(r'l=(.*?)m=')
     mnum_re = re.compile(r'm=(.*?)\)')
-    wfc_lines = out_info_dict['wfc_lines']
-    state_lines = [out_file[wfc_line] for wfc_line in wfc_lines]
+    state_lines = [out_file_lines[wfc_line] for wfc_line in statelines_idx]
     state_dicts = []
     for state_line in state_lines:
         try:
@@ -48,7 +50,7 @@ def find_orbitals_from_statelines(out_info_dict):
 
     # here is some logic to figure out the value of radial_nodes to use
     new_state_dicts = []
-    for i in range(len(state_dicts)):
+    for i in range(len(state_dicts)):  # pylint: disable=consider-using-enumerate
         radial_nodes = 0
         state_dict = state_dicts[i].copy()
         for j in range(i - 1, -1, -1):
@@ -59,92 +61,84 @@ def find_orbitals_from_statelines(out_info_dict):
     state_dicts = new_state_dicts
 
     # here is some logic to assign positions based on the atom_index
-    structure = out_info_dict['structure']
     for state_dict in state_dicts:
         site_index = state_dict.pop('atomnum')
-        state_dict['position'] = structure.sites[site_index].position
+        state_dict['position'] = positions[site_index]
 
     # here we set the resulting state_dicts to a new set of orbitals
     orbitals = []
-    RealhydrogenOrbital = OrbitalFactory('realhydrogen')
+    real_hydrogen_orbital = OrbitalFactory('realhydrogen')
     for state_dict in state_dicts:
-        orbitals.append(RealhydrogenOrbital(**state_dict))
+        orbitals.append(real_hydrogen_orbital(**state_dict))
 
     return orbitals
 
 
-def spin_dependent_subparser(out_info_dict):
-    """
-    This find the projection and bands arrays from the out_file and
-    out_info_dict. Used to handle the different possible spin-cases in
-    a convenient manner.
+def spin_dependent_subparser(info_dict, parent_kpoints):
+    """Create projection and bands arrays from the parsed data and parent calculation.
+
+    Used to handle the different possible spin-cases in a convenient manner.
 
-    :param out_info_dict: contains various technical internals useful in parsing
-    :return: ProjectionData, BandsData parsed from out_file
+    :param info_dict: dictionary of input and parsed data
+    :param parent_kpoints: kpoints from the parent calculation
+    :return: ProjectionData, BandsData
     """
 
-    out_file = out_info_dict['out_file']
-    spin_down = out_info_dict['spin_down']
-    od = out_info_dict  #using a shorter name for convenience
+    out_file_lines = info_dict['out_file_lines']
+    spin_down, k_states, num_bands = info_dict['spin_down'], info_dict['k_states'], info_dict['num_bands']
+
     #   regular expressions needed for later parsing
-    WaveFraction1_re = re.compile(r'\=(.*?)\*')  # state composition 1
-    WaveFractionremain_re = re.compile(r'\+(.*?)\*')  # state comp 2
-    FunctionId_re = re.compile(r'\#(.*?)\]')  # state identity
+    regex_wavefrac1 = re.compile(r'\=(.*?)\*')  # state composition 1
+    regex_wavefrac2 = re.compile(r'\+(.*?)\*')  # state comp 2
+    regex_func_id = re.compile(r'\#(.*?)\]')  # state identity
     # primes arrays for the later parsing
-    num_wfc = len(od['wfc_lines'])
-    bands = np.zeros([od['k_states'], od['num_bands']])
-    projection_arrays = np.zeros([od['k_states'], od['num_bands'], num_wfc])
+    bands = np.zeros([k_states, num_bands])
+    projection_arrays = np.zeros([k_states, num_bands, len(info_dict['wfc_lines'])])
 
     try:
-        for i in range(od['k_states']):
+        for i in range(k_states):
             if spin_down:
-                i += od['k_states']
+                i += k_states
             # grabs band energy
-            for j in range(i * od['num_bands'], (i + 1) * od['num_bands'], 1):
-                out_ind = od['e_lines'][j]
+            for j in range(i * num_bands, (i + 1) * num_bands, 1):
+                out_ind = info_dict['e_lines'][j]
                 try:
                     # post ~6.3 output format "e ="
-                    val = out_file[out_ind].split()[2]
+                    val = out_file_lines[out_ind].split()[2]
                     float(val)
                 except ValueError:
                     # pre ~6.3 output format? "==== e("
-                    val = out_file[out_ind].split()[4]
-                bands[i % od['k_states']][j % od['num_bands']] = val
+                    val = out_file_lines[out_ind].split()[4]
+                bands[i % k_states][j % num_bands] = val
                 #subloop grabs pdos
                 wave_fraction = []
                 wave_id = []
-                for k in range(od['e_lines'][j] + 1, od['psi_lines'][j], 1):
-                    out_line = out_file[k]
-                    wave_fraction += WaveFraction1_re.findall(out_line)
-                    wave_fraction += WaveFractionremain_re.findall(out_line)
-                    wave_id += FunctionId_re.findall(out_line)
+                for k in range(info_dict['e_lines'][j] + 1, info_dict['psi_lines'][j], 1):
+                    out_line = out_file_lines[k]
+                    wave_fraction += regex_wavefrac1.findall(out_line)
+                    wave_fraction += regex_wavefrac2.findall(out_line)
+                    wave_id += regex_func_id.findall(out_line)
                 if len(wave_id) != len(wave_fraction):
                     raise IndexError
-                for l in range(len(wave_id)):
+                for l in range(len(wave_id)):  # pylint: disable=consider-using-enumerate,invalid-name
                     wave_id[l] = int(wave_id[l])
                     wave_fraction[l] = float(wave_fraction[l])
                     #sets relevant values in pdos_array
-                    projection_arrays[i % od['k_states']][j % od['num_bands']][wave_id[l] - 1] = wave_fraction[l]
+                    projection_arrays[i % k_states][j % num_bands][wave_id[l] - 1] = wave_fraction[l]
     except IndexError:
-        raise QEOutputParsingError('the standard out file does not ' 'comply with the official ' 'documentation.')
+        raise QEOutputParsingError('the standard out file does not comply with the official documentation.')
 
     bands_data = BandsData()
-    # Attempts to retrieve the kpoints from the parent calc
-    parent_calc = out_info_dict['parent_calc']
-    try:
-        parent_kpoints = parent_calc.get_incoming(link_label_filter='kpoints').one().node
-    except ValueError:
-        raise QEOutputParsingError('The parent had no input kpoints! Cannot parse from this!')
     try:
-        if len(od['k_vect']) != len(parent_kpoints.get_kpoints()):
+        if len(info_dict['k_vect']) != len(parent_kpoints.get_kpoints()):
             raise AttributeError
         bands_data.set_kpointsdata(parent_kpoints)
     except AttributeError:
-        bands_data.set_kpoints(od['k_vect'].astype(float))
+        bands_data.set_kpoints(info_dict['k_vect'].astype(float))
 
     bands_data.set_bands(bands, units='eV')
 
-    orbitals = out_info_dict['orbitals']
+    orbitals = info_dict['orbitals']
     if len(orbitals) != np.shape(projection_arrays[0, 0, :])[0]:
         raise QEOutputParsingError(
             'There was an internal parsing error, '
@@ -160,9 +154,9 @@ def spin_dependent_subparser(out_info_dict):
         if np.shape(projection) != np.shape(bands):
             raise AttributeError('Projections not the same shape as the bands')
 
-    #insert here some logic to assign pdos to the orbitals
-    pdos_arrays = spin_dependent_pdos_subparser(out_info_dict)
-    energy_arrays = [out_info_dict['energy']] * len(orbitals)
+    # insert here some logic to assign pdos to the orbitals
+    pdos_arrays = spin_dependent_pdos_subparser(info_dict)
+    energy_arrays = [info_dict['energy']] * len(orbitals)
     projection_data.set_projectiondata(
         orbitals,
         list_of_projections=projections,
@@ -175,12 +169,13 @@ def spin_dependent_subparser(out_info_dict):
 
 
 def natural_sort_key(sort_key, _nsre=re.compile('([0-9]+)')):
-    """
-    Pass to ``key`` for ``str.sort`` to achieve natural sorting.
+    """Convert a sorting key to one suitable for natural sorting.
+
     For example, ``["2", "11", "1"]`` will be
     sorted to ``["1", "2", "11"]`` instead of ``["1", "11", "2"]``
+
     :param sort_key: Original key to be processed
-    :return: A list of string and integers.
+    :return: list[str or int]
     """
     keys = []
     for text in _nsre.split(sort_key):
@@ -192,10 +187,9 @@ def natural_sort_key(sort_key, _nsre=re.compile('([0-9]+)')):
 
 
 def spin_dependent_pdos_subparser(out_info_dict):
-    """
-    Finds and labels the pdos arrays associated with the out_info_dict
+    """Find and label the pdos arrays associated with the out_info_dict.
 
-    :param out_info_dict: contains various technical internals useful in parsing
+    :param out_info_dict: dictionary of input and parsed data
     :return: (pdos_name, pdos_array) tuples for all the specific pdos
     """
     spin = out_info_dict['spin']
@@ -210,8 +204,6 @@ def spin_dependent_pdos_subparser(out_info_dict):
     else:
         mult_factor = 1
         first_array = 2
-    mf = mult_factor
-    fa = first_array
     pdos_file_names = [k for k in pdos_atm_array_dict]
     pdos_file_names.sort(key=natural_sort_key)
     out_arrays = []
@@ -219,12 +211,32 @@ def spin_dependent_pdos_subparser(out_info_dict):
     # both are produced in the same order (thus the sorted file_names)
     for name in pdos_file_names:
         this_array = pdos_atm_array_dict[name]
-        for i in range(fa, np.shape(this_array)[1], mf):
+        for i in range(first_array, np.shape(this_array)[1], mult_factor):
             out_arrays.append(this_array[:, i])
 
     return out_arrays
 
 
+def scan_stdout(lines):
+    """Iterate through the projwf.x stdout file and record indices of required lines."""
+    data = {'k_lines': [], 'e_lines': [], 'psi_lines': [], 'wfc_lines': [], 'lowdin_lines': []}
+    for i, line in enumerate(lines):
+        if 'k =' in line:
+            data['k_lines'].append(i)
+        if '==== e(' in line or line.strip().startswith('e ='):
+            # The energy format in output was changed in QE6.3
+            # this check supports old and new format
+            data['e_lines'].append(i)
+        if '|psi|^2' in line:
+            data['psi_lines'].append(i)
+        if 'state #' in line:
+            data['wfc_lines'].append(i)
+        if line.strip() == 'Lowdin Charges:':
+            data['lowdin_lines'].append(i)
+
+    return data
+
+
 class ProjwfcParser(Parser):
     """
     This class is the implementation of the Parser class for projwfc.x in
@@ -245,17 +257,23 @@ def parse(self, **kwargs):
         except NotExistent:
             return self.exit_codes.ERROR_NO_RETRIEVED_FOLDER
 
+        # Check that the retrieved temporary folder is there
+        try:
+            temp_folder = kwargs['retrieved_temporary_folder']
+        except KeyError:
+            return self.exit_codes.ERROR_NO_RETRIEVED_TEMPORARY_FOLDER
+
         # Read standard out
         try:
-            filename_stdout = self.node.get_option('output_filename')  # or get_attribute(), but this is clearer
-            with out_folder.open(filename_stdout, 'r') as fil:
-                out_file = fil.readlines()
+            filename_stdout = self.node.get_option('output_filename')
+            with open(os.path.join(temp_folder, filename_stdout), 'r') as fil:
+                out_file_lines = fil.readlines()
         except OSError:
             return self.exit_codes.ERROR_READING_OUTPUT_FILE
 
+        # check that the computation completed successfully
         job_done = False
-        for i in range(len(out_file)):
-            line = out_file[-i]
+        for line in reversed(out_file_lines):
             if 'JOB DONE' in line:
                 job_done = True
                 break
@@ -264,7 +282,7 @@ def parse(self, **kwargs):
             return self.exit_codes.ERROR_JOB_NOT_DONE
 
         # Parse basic info and warnings, and output them as output_parmeters
-        parsed_data = parse_raw_out_basic(out_file, 'PROJWFC')
+        parsed_data = parse_raw_out_basic(out_file_lines, 'PROJWFC')
         for message in parsed_data['warnings']:
             self.logger.error(message)
         output_params = Dict(dict=parsed_data)
@@ -292,52 +310,82 @@ def parse(self, **kwargs):
 
         # finding the bands and projections
         # we create a dictionary the progressively accumulates more info
-        out_info_dict = {}
-        out_info_dict['out_file'] = out_file
-        out_info_dict['energy'] = energy
-        out_info_dict['pdos_atm_array_dict'] = pdos_atm_array_dict
+        out_info_dict = {'out_file_lines': out_file_lines, 'energy': energy, 'pdos_atm_array_dict': pdos_atm_array_dict}
+        out_info_dict.update(scan_stdout(out_file_lines))
         try:
             new_nodes_list = self._parse_bands_and_projections(out_info_dict)
         except QEOutputParsingError as err:
             self.logger.error('Error parsing bands and projections: {}'.format(err))
             traceback.print_exc()
             return self.exit_codes.ERROR_PARSING_PROJECTIONS
+
         for linkname, node in new_nodes_list:
             self.out(linkname, node)
 
-        Dos_out = XyData()
-        Dos_out.set_x(energy, 'Energy', 'eV')
-        Dos_out.set_y(dos, 'Dos', 'states/eV')
-        self.out('Dos', Dos_out)
+        dos_out = XyData()
+        dos_out.set_x(energy, 'Energy', 'eV')
+        dos_out.set_y(dos, 'Dos', 'states/eV')
+        self.out('Dos', dos_out)
+
+        try:
+            node = self._parse_lodwin_charges(out_info_dict)
+        except IOError as err:
+            self.logger.error('Error parsing Lowdin charges: {}'.format(err))
+            traceback.print_exc()
+            return self.exit_codes.ERROR_PARSING_LOWDIN
+
+        self.out('lowdin', node)
+
+    def retrieve_parent_node(self, linkname, outgoing=True):
+        """Retrieve a node from the calculation that created the parent_folder."""
+        try:
+            parent_calc = (
+                self.node.inputs.parent_folder.get_incoming(node_class=CalcJobNode,
+                                                            link_type=LinkType.CREATE).one().node
+            )
+        except ValueError as err:
+            raise QEOutputParsingError('Could not retrieve incoming calculation of parent_folder: {}'.format(err))
+        try:
+            if outgoing:
+                node = parent_calc.get_outgoing(link_label_filter=linkname).one().node
+            else:
+                node = parent_calc.get_incoming(link_label_filter=linkname).one().node
+        except ValueError:
+            raise QEOutputParsingError(
+                'The parent calculation had no {} {} node'.format('outgoing' if outgoing else 'incoming', linkname)
+            )
+        return node
+
+    def _parse_lodwin_charges(self, out_info_dict):
+        """Parse the Lodwin charge data from the output file."""
+        data, spill_parameter = parse_lowdin_charges(out_info_dict['out_file_lines'], out_info_dict['lowdin_lines'])
+        # we store the uuid of the structure that these charges relate to, which will be
+        # an input if the PwCalculation was an scf/nscf or output if relax/vc-relax
+        try:
+            structure = self.retrieve_parent_node('output_structure', outgoing=True)
+        except QEOutputParsingError:
+            structure = self.retrieve_parent_node('structure', outgoing=False)
+        try:
+            site_data = [data[i + 1] for i in range(len(structure.sites))]
+        except KeyError:
+            raise IOError('The lowdin atom numbers do not match those in the input structure.')
+        return Dict(
+            dict={
+                'site_data': site_data,
+                'spill_parameter': spill_parameter,
+                'structure_uuid': structure.uuid
+            }
+        )
 
     def _parse_bands_and_projections(self, out_info_dict):
-        """
-        Function that parses the standard output into bands and projection data.
+        """Parse the standard output into bands and projection data.
 
-        :param out_info_dict: used to pass technical internal variables
-                              to helper functions in compact form
-        :return: append_nodes_list a list containing BandsData and
-                 ProjectionData parsed from standard_out
+        :param out_info_dict: used to pass technical internal variables to helper functions in compact form
+        :return: a list of (linkname, node) containing BandsData and ProjectionData parsed from standard_out
         """
-        out_file = out_info_dict['out_file']  # Note: we expect a list of lines
-        out_info_dict['k_lines'] = []
-        out_info_dict['e_lines'] = []
-        out_info_dict['psi_lines'] = []
-        out_info_dict['wfc_lines'] = []
+        out_file_lines = out_info_dict['out_file_lines']
         append_nodes_list = []
 
-        for i, line in enumerate(out_file):
-            if 'k =' in line:
-                out_info_dict['k_lines'].append(i)
-            if '==== e(' in line or line.strip().startswith('e ='):
-                # The energy format in output was changed in QE6.3
-                # this check supports old and new format
-                out_info_dict['e_lines'].append(i)
-            if '|psi|^2' in line:
-                out_info_dict['psi_lines'].append(i)
-            if 'state #' in line:
-                out_info_dict['wfc_lines'].append(i)
-
         # Basic check
         if len(out_info_dict['e_lines']) != len(out_info_dict['psi_lines']):
             raise QEOutputParsingError('e-lines and psi-lines are in different number')
@@ -346,32 +394,9 @@ def _parse_bands_and_projections(self, out_info_dict):
         # calculates the number of bands
         out_info_dict['num_bands'] = len(out_info_dict['psi_lines']) // len(out_info_dict['k_lines'])
 
-        # Uses the parent input parameters, and checks if the parent used
-        # spin calculations. Try to replace with a query, if possible.
-        try:
-            parent_calc = (
-                self.node.inputs.parent_folder.get_incoming(node_class=CalcJobNode,
-                                                            link_type=LinkType.CREATE).one().node
-            )
-        except ValueError as e:
-            raise QEOutputParsingError('Could not get parent calculation of input folder: {}'.format(e))
-        out_info_dict['parent_calc'] = parent_calc
-        try:
-            parent_param = parent_calc.get_outgoing(link_label_filter='output_parameters').one().node
-        except ValueError:
-            raise QEOutputParsingError('The parent had no output_parameters! Cannot parse from this!')
-        try:
-            structure = parent_calc.get_incoming(link_label_filter='structure').one().node
-        except ValueError:
-            raise QEOutputParsingError('The parent had no input structure! Cannot parse from this!')
-        try:
-            nspin = parent_param.get_dict()['number_of_spin_components']
-            if nspin != 1:
-                spin = True
-            else:
-                spin = False
-        except KeyError:
-            spin = False
+        # Use the parent input parameters, to checks if the calculation has spin.
+        parent_param = self.retrieve_parent_node('output_parameters')
+        spin = parent_param.get_dict().get('number_of_spin_components', 1) != 1
         out_info_dict['spin'] = spin
 
         # changes k-numbers to match spin
@@ -382,12 +407,15 @@ def _parse_bands_and_projections(self, out_info_dict):
                 raise QEOutputParsingError('Internal formatting error regarding spin')
             out_info_dict['k_states'] = out_info_dict['k_states'] // 2
 
-        #   adds in the k-vector for each kpoint
-        k_vect = [out_file[out_info_dict['k_lines'][i]].split()[2:] for i in range(out_info_dict['k_states'])]
+        # adds in the k-vector for each kpoint
+        k_vect = [out_file_lines[out_info_dict['k_lines'][i]].split()[2:] for i in range(out_info_dict['k_states'])]
         out_info_dict['k_vect'] = np.array(k_vect)
-        out_info_dict['structure'] = structure
-        out_info_dict['orbitals'] = find_orbitals_from_statelines(out_info_dict)
+        structure = self.retrieve_parent_node('structure', outgoing=False)
+        out_info_dict['orbitals'] = find_orbitals_from_statelines(
+            out_file_lines, out_info_dict['wfc_lines'], [s.position for s in structure.sites]
+        )
 
+        kpoints = self.retrieve_parent_node('kpoints', outgoing=False)
         if spin:
             # I had to guess what the ordering of the spin is, because
             # the projwfc.x documentation doesn't say, but looking at the
@@ -402,14 +430,14 @@ def _parse_bands_and_projections(self, out_info_dict):
             # spin up states are written first, then spin down
             #
             out_info_dict['spin_down'] = False
-            bands_data1, projection_data1 = spin_dependent_subparser(out_info_dict)
+            bands_data1, projection_data1 = spin_dependent_subparser(out_info_dict, kpoints)
             append_nodes_list += [('projections_up', projection_data1), ('bands_up', bands_data1)]
             out_info_dict['spin_down'] = True
-            bands_data2, projection_data2 = spin_dependent_subparser(out_info_dict)
+            bands_data2, projection_data2 = spin_dependent_subparser(out_info_dict, kpoints)
             append_nodes_list += [('projections_down', projection_data2), ('bands_down', bands_data2)]
         else:
             out_info_dict['spin_down'] = False
-            bands_data, projection_data = spin_dependent_subparser(out_info_dict)
+            bands_data, projection_data = spin_dependent_subparser(out_info_dict, kpoints)
             append_nodes_list += [('projections', projection_data), ('bands', bands_data)]
 
         return append_nodes_list
diff --git a/tests/calculations/test_projwfc.py b/tests/calculations/test_projwfc.py
new file mode 100644
index 000000000..73a74e2f6
--- /dev/null
+++ b/tests/calculations/test_projwfc.py
@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+# pylint: disable=unused-argument
+"""Tests for the `ProjwfcCalculation` class."""
+from __future__ import absolute_import
+
+from aiida import orm
+from aiida.common import datastructures
+
+from aiida_quantumespresso.utils.resources import get_default_options
+
+
+def test_projwfc_default(
+    fixture_database, fixture_computer_localhost, fixture_sandbox_folder, generate_calc_job, generate_code_localhost,
+    file_regression
+):
+    """Test a default `ProjwfcCalculation`."""
+    entry_point_name = 'quantumespresso.projwfc'
+
+    parameters = {'PROJWFC': {'emin': -1, 'emax': 1, 'DeltaE': 0.01, 'ngauss': 0, 'degauss': 0.01}}
+    parent_folder = orm.RemoteData(computer=fixture_computer_localhost, remote_path='path/on/remote')
+    parent_folder.store()
+
+    inputs = {
+        'code': generate_code_localhost(entry_point_name, fixture_computer_localhost),
+        'parameters': orm.Dict(dict=parameters),
+        'parent_folder': parent_folder,
+        'metadata': {
+            'options': get_default_options()
+        }
+    }
+
+    calc_info = generate_calc_job(fixture_sandbox_folder, entry_point_name, inputs)
+    code_info = calc_info.codes_info[0]
+
+    # Check the attributes of the returned `CodeInfo`
+    assert isinstance(code_info, datastructures.CodeInfo)
+    assert code_info.stdin_name == 'aiida.in'
+    assert code_info.stdout_name == 'aiida.out'
+    assert code_info.cmdline_params == []
+    # Check the attributes of the returned `CalcInfo`
+    assert isinstance(calc_info, datastructures.CalcInfo)
+    assert sorted(calc_info.local_copy_list) == sorted([])
+    assert sorted(calc_info.remote_copy_list
+                  ) == sorted([(parent_folder.computer.uuid, 'path/on/remote/./out/', './out/')])
+    assert sorted(calc_info.retrieve_list) == sorted(['aiida.pdos*'])
+    assert sorted(calc_info.retrieve_temporary_list) == sorted(['aiida.out'])
+
+    with fixture_sandbox_folder.open('aiida.in') as handle:
+        input_written = handle.read()
+
+    # Checks on the files written to the sandbox folder as raw input
+    assert sorted(fixture_sandbox_folder.get_content_list()) == sorted(['aiida.in'])
+    file_regression.check(input_written, encoding='utf-8', extension='.in')
diff --git a/tests/calculations/test_projwfc/test_projwfc_default.in b/tests/calculations/test_projwfc/test_projwfc_default.in
new file mode 100644
index 000000000..592c7bfd8
--- /dev/null
+++ b/tests/calculations/test_projwfc/test_projwfc_default.in
@@ -0,0 +1,15 @@
+&PROJWFC
+  degauss =   1.0000000000d-02
+  deltae =   1.0000000000d-02
+  emax = 1
+  emin = -1
+  kresolveddos = .false.
+  lbinary_data = .false.
+  lsym = .true.
+  lwrite_overlaps = .false.
+  ngauss = 0
+  outdir = './out/'
+  plotboxes = .false.
+  prefix = 'aiida'
+  tdosinboxes = .false.
+/
diff --git a/tests/conftest.py b/tests/conftest.py
index 8f719dea2..1a1ece237 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -58,7 +58,8 @@ def fixture_computer_localhost(fixture_work_directory):
         hostname='localhost',
         transport_type='local',
         scheduler_type='direct',
-        workdir=fixture_work_directory).store()
+        workdir=fixture_work_directory
+    ).store()
     computer.set_default_mpiprocs_per_machine(1)
     yield computer
 
@@ -126,7 +127,7 @@ def _generate_calc_job_node(entry_point_name, computer, test_name=None, inputs=N
         node.set_option('max_wallclock_seconds', 1800)
 
         if attributes:
-            node.set_attributes(attributes)
+            node.set_attribute_many(attributes)
 
         if inputs:
             for link_label, input_node in flatten_inputs(inputs):
@@ -137,7 +138,9 @@ def _generate_calc_job_node(entry_point_name, computer, test_name=None, inputs=N
 
         if test_name is not None:
             basepath = os.path.dirname(os.path.abspath(__file__))
-            filepath = os.path.join(basepath, 'parsers', 'fixtures', entry_point_name[len('quantumespresso.'):], test_name)
+            filepath = os.path.join(
+                basepath, 'parsers', 'fixtures', entry_point_name[len('quantumespresso.'):], test_name
+            )
 
             retrieved = orm.FolderData()
             retrieved.put_object_from_tree(filepath)
@@ -234,3 +237,59 @@ def _generate_parser(entry_point_name):
         return ParserFactory(entry_point_name)
 
     return _generate_parser
+
+
+@pytest.fixture
+def parse_from_node():
+    """Fixture to parse calcjob outputs directly from a `CalcJobNode`."""
+
+    # TODO this can be removed once aiidateam/aiida-core#3061 is implemented  # pylint: disable=fixme
+    def _parse_from_node(cls, node, store_provenance=True, retrieved_temp=None):
+        """Parse the outputs directly from the `CalcJobNode`.
+
+        :param cls: a `Parser` instance
+        :param node: a `CalcJobNode` instance
+        :param store_provenance: bool, if True will store the parsing as a `CalcFunctionNode` in the provenance
+        :param retrieved_temp: None or str, abspath to the temporary folder.
+        :return: a tuple of the parsed results and the `CalcFunctionNode` representing the process of parsing
+        """
+        from aiida.engine import calcfunction, Process
+        from aiida.orm import Str
+
+        parser = cls(node=node)
+
+        @calcfunction
+        def parse_calcfunction(**kwargs):
+            """A wrapper function that will turn calling the `Parser.parse` method into a `CalcFunctionNode`.
+
+            :param kwargs: keyword arguments that are passed to `Parser.parse` after it has been constructed
+            """
+            if 'retrieved_temporary_folder' in kwargs:
+                string = kwargs.pop('retrieved_temporary_folder').value
+                kwargs['retrieved_temporary_folder'] = string
+
+            exit_code = parser.parse(**kwargs)
+            outputs = parser.outputs
+
+            if exit_code and exit_code.status:
+                process = Process.current()
+                process.out_many(outputs)
+                return exit_code
+
+            return dict(outputs)
+
+        inputs = {'metadata': {'store_provenance': store_provenance}}
+        inputs.update(parser.get_outputs_for_parsing())
+        if retrieved_temp is not None:
+            inputs['retrieved_temporary_folder'] = Str(retrieved_temp)
+
+        return parse_calcfunction.run_get_node(**inputs)
+
+    return _parse_from_node
+
+
+@pytest.fixture
+def parser_fixture_path():
+    """Fixture to obtain the absolute path of the ``test/parsers/fixtures`` folder."""
+    basepath = os.path.dirname(os.path.abspath(__file__))
+    return os.path.join(basepath, 'parsers', 'fixtures')
diff --git a/tests/parsers/fixtures/projwfc/spin/aiida.out b/tests/parsers/fixtures/projwfc/spin/aiida.out
new file mode 100644
index 000000000..7d19eeb0a
--- /dev/null
+++ b/tests/parsers/fixtures/projwfc/spin/aiida.out
@@ -0,0 +1,175 @@
+
+
+# mpiexec: full path to program is /apps/espresso/5.3.0/bin/projwfc.x
+# mpiexec: program arguments are: -npools 2 -i generic.pdos.in
+# mpiexec: 32 ranks allocated via PBS select
+# mpiexec: 1 OpenMP threads / rank allocated by PBS select
+# mpiexec: Job has exclusive use of the allocated nodes. Enabling process-pinning
+# mpiexec: launching program...
+
+# mpiexec: machinefile configured as:
+cx1-136-2-1.cx1.hpc.ic.ac.uk:16
+cx1-136-3-1.cx1.hpc.ic.ac.uk:16
+#
+
+     Program PROJWFC v.5.3.0 (svn rev. 11974) starts on 26Mar2018 at 20:12:51 
+
+     This program is part of the open-source Quantum ESPRESSO suite
+     for quantum simulation of materials; please cite
+         "P. Giannozzi et al., J. Phys.:Condens. Matter 21 395502 (2009);
+          URL http://www.quantum-espresso.org", 
+     in publications or presentations arising from this work. More details at
+     http://www.quantum-espresso.org/quote
+
+     Parallel version (MPI & OpenMP), running on      32 processor cores
+     Number of MPI processes:                32
+     Threads/MPI process:                     1
+     K-points division:     npool     =       2
+     R & G space division:  proc/nbgrp/npool/nimage =      16
+
+     Reading data from directory:
+     ./iron_sulfide.save
+
+   Info: using nr1, nr2, nr3 values from input
+
+   Info: using nr1, nr2, nr3 values from input
+
+     IMPORTANT: XC functional enforced from input :
+     Exchange-correlation      =  SLA  PW   PBX  PBC ( 1  4  3  4 0 0)
+     Any further DFT definition will be discarded
+     Please, verify this is what you really want
+
+     Message from routine readpp:
+     file ./iron_sulfide.save/Fe.pbe-spn-kjpaw_psl.0.2.1.UPF not found
+               file Fe.pbe-spn-kjpaw_psl.0.2.1.UPF: wavefunction(s)  3D renormalized
+     Message from routine readpp:
+     file ./iron_sulfide.save/Fe.pbe-spn-kjpaw_psl.0.2.1.UPF not found
+               file Fe.pbe-spn-kjpaw_psl.0.2.1.UPF: wavefunction(s)  3D renormalized
+     Message from routine readpp:
+     file ./iron_sulfide.save/S.pbe-n-kjpaw_psl.0.1.UPF not found
+               file S.pbe-n-kjpaw_psl.0.1.UPF: wavefunction(s)  3P renormalized
+ 
+     Parallelization info
+     --------------------
+     sticks:   dense  smooth     PW     G-vecs:    dense   smooth      PW
+     Min         119     119     34                 3933     3933     608
+     Max         122     122     35                 3935     3935     611
+     Sum        1917    1917    553                62949    62949    9755
+ 
+     Generating pointlists ...
+     new r_m :   0.2143 (alat units)  1.9656 (a.u.) for type    1
+     new r_m :   0.2143 (alat units)  1.9656 (a.u.) for type    2
+     new r_m :   0.2143 (alat units)  1.9656 (a.u.) for type    3
+
+     Check: negative/imaginary core charge=   -0.000232    0.000000
+
+     negative rho (up, down):  0.000E+00 3.075E-03
+
+     Gaussian broadening (read from input): ngauss,degauss=   0    0.010000
+
+
+     Calling projwave .... 
+
+     Atomic states used for projection
+     (read from pseudopotential files):
+
+     state #   1: atom   1 (Fe1), wfc  1 (l=0 m= 1)
+     state #   2: atom   1 (Fe1), wfc  2 (l=0 m= 1)
+     state #   3: atom   1 (Fe1), wfc  3 (l=1 m= 1)
+     state #   4: atom   1 (Fe1), wfc  3 (l=1 m= 2)
+     state #   5: atom   1 (Fe1), wfc  3 (l=1 m= 3)
+     state #   6: atom   1 (Fe1), wfc  4 (l=1 m= 1)
+     state #   7: atom   1 (Fe1), wfc  4 (l=1 m= 2)
+     state #   8: atom   1 (Fe1), wfc  4 (l=1 m= 3)
+     state #   9: atom   1 (Fe1), wfc  5 (l=2 m= 1)
+     state #  10: atom   1 (Fe1), wfc  5 (l=2 m= 2)
+
+<NOTE The central section of this file has been removed, for testing purposes.>
+
+==== e(  51) =    14.21482 eV ==== 
+     psi = 0.328*[#  35]+0.328*[#  48]+0.106*[#  58]+0.106*[#  62]+0.027*[#  55]+
+          +0.027*[#  56]+0.027*[#  67]+0.027*[#  68]+0.007*[#  39]+0.007*[#  52]+
+          +0.002*[#  10]+0.002*[#  11]+0.002*[#  23]+0.002*[#  24]+
+    |psi|^2 = 0.998
+==== e(  52) =    14.21483 eV ==== 
+     psi = 0.657*[#   9]+0.212*[#  66]+0.054*[#  59]+0.054*[#  64]+0.015*[#  26]+
+          +0.003*[#  37]+0.003*[#  49]+
+    |psi|^2 = 0.998
+
+Lowdin Charges: 
+
+     Atom #   1: total charge =  15.8849, s =  2.4347, p =  7.0900, d =  6.3602, 
+                 spin up      =   9.7893, s =  1.2506, 
+                 spin up      =   9.7893, p =  3.5691, pz=  1.1900, px=  1.1895, py=  1.1895, 
+                 spin up      =   9.7893, d =  4.9696, dz2=  0.9893, dxz=  0.9968, dyz=  0.9968, dx2-y2=  0.9899, dxy=  0.9968, 
+                 spin down    =   6.0955, s =  1.1840, 
+                 spin down    =   6.0955, p =  3.5208, pz=  1.1740, px=  1.1734, py=  1.1734, 
+                 spin down    =   6.0955, d =  1.3907, dz2=  0.1575, dxz=  0.3592, dyz=  0.3592, dx2-y2=  0.1557, dxy=  0.3590, 
+                 polarization =   3.6938, s =  0.0666, p =  0.0483, d =  3.5789, 
+     Atom #   2: total charge =  15.8864, s =  2.4347, p =  7.0900, d =  6.3618, 
+                 spin up      =   9.7896, s =  1.2506, 
+                 spin up      =   9.7896, p =  3.5691, pz=  1.1900, px=  1.1895, py=  1.1895, 
+                 spin up      =   9.7896, d =  4.9699, dz2=  0.9895, dxz=  0.9968, dyz=  0.9968, dx2-y2=  0.9901, dxy=  0.9968, 
+                 spin down    =   6.0968, s =  1.1840, 
+                 spin down    =   6.0968, p =  3.5208, pz=  1.1740, px=  1.1734, py=  1.1734, 
+                 spin down    =   6.0968, d =  1.3919, dz2=  0.1575, dxz=  0.3596, dyz=  0.3596, dx2-y2=  0.1557, dxy=  0.3596, 
+                 polarization =   3.6928, s =  0.0666, p =  0.0483, d =  3.5779, 
+     Atom #   3: total charge =  15.8863, s =  2.4347, p =  7.0900, d =  6.3617, 
+                 spin up      =   9.7894, s =  1.2506, 
+                 spin up      =   9.7894, p =  3.5691, pz=  1.1900, px=  1.1895, py=  1.1895, 
+                 spin up      =   9.7894, d =  4.9697, dz2=  0.9893, dxz=  0.9968, dyz=  0.9968, dx2-y2=  0.9900, dxy=  0.9968, 
+                 spin down    =   6.0968, s =  1.1840, 
+                 spin down    =   6.0968, p =  3.5208, pz=  1.1740, px=  1.1734, py=  1.1734, 
+                 spin down    =   6.0968, d =  1.3920, dz2=  0.1575, dxz=  0.3597, dyz=  0.3601, dx2-y2=  0.1557, dxy=  0.3589, 
+                 polarization =   3.6926, s =  0.0666, p =  0.0483, d =  3.5777, 
+     Atom #   4: total charge =  15.8863, s =  2.4347, p =  7.0900, d =  6.3617, 
+                 spin up      =   9.7894, s =  1.2506, 
+                 spin up      =   9.7894, p =  3.5691, pz=  1.1900, px=  1.1895, py=  1.1895, 
+                 spin up      =   9.7894, d =  4.9697, dz2=  0.9893, dxz=  0.9968, dyz=  0.9968, dx2-y2=  0.9900, dxy=  0.9968, 
+                 spin down    =   6.0968, s =  1.1840, 
+                 spin down    =   6.0968, p =  3.5208, pz=  1.1740, px=  1.1734, py=  1.1734, 
+                 spin down    =   6.0968, d =  1.3920, dz2=  0.1575, dxz=  0.3601, dyz=  0.3597, dx2-y2=  0.1557, dxy=  0.3589, 
+                 polarization =   3.6926, s =  0.0666, p =  0.0483, d =  3.5777, 
+     Atom #   5: total charge =   6.0441, s =  1.5290, p =  4.5151, d =  0.0000, 
+                 spin up      =   3.1726, s =  0.7848, 
+                 spin up      =   3.1726, p =  2.3878, pz=  0.7950, px=  0.7964, py=  0.7964, 
+                 spin up      =   3.1726, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 spin down    =   2.8715, s =  0.7441, 
+                 spin down    =   2.8715, p =  2.1274, pz=  0.7093, px=  0.7090, py=  0.7090, 
+                 spin down    =   2.8715, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 polarization =   0.3011, s =  0.0407, p =  0.2604, d =  0.0000, 
+     Atom #   6: total charge =   6.0440, s =  1.5290, p =  4.5150, d =  0.0000, 
+                 spin up      =   3.1725, s =  0.7848, 
+                 spin up      =   3.1725, p =  2.3877, pz=  0.7947, px=  0.7963, py=  0.7967, 
+                 spin up      =   3.1725, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 spin down    =   2.8715, s =  0.7441, 
+                 spin down    =   2.8715, p =  2.1274, pz=  0.7093, px=  0.7090, py=  0.7090, 
+                 spin down    =   2.8715, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 polarization =   0.3010, s =  0.0407, p =  0.2603, d =  0.0000, 
+     Atom #   7: total charge =   6.0440, s =  1.5290, p =  4.5150, d =  0.0000, 
+                 spin up      =   3.1725, s =  0.7848, 
+                 spin up      =   3.1725, p =  2.3877, pz=  0.7947, px=  0.7967, py=  0.7963, 
+                 spin up      =   3.1725, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 spin down    =   2.8715, s =  0.7441, 
+                 spin down    =   2.8715, p =  2.1274, pz=  0.7093, px=  0.7090, py=  0.7090, 
+                 spin down    =   2.8715, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 polarization =   0.3010, s =  0.0407, p =  0.2603, d =  0.0000, 
+     Atom #   8: total charge =   6.0440, s =  1.5290, p =  4.5150, d =  0.0000, 
+                 spin up      =   3.1725, s =  0.7848, 
+                 spin up      =   3.1725, p =  2.3877, pz=  0.7946, px=  0.7965, py=  0.7965, 
+                 spin up      =   3.1725, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 spin down    =   2.8715, s =  0.7441, 
+                 spin down    =   2.8715, p =  2.1274, pz=  0.7093, px=  0.7090, py=  0.7090, 
+                 spin down    =   2.8715, d =  0.0000, dz2=  0.0000, dxz=  0.0000, dyz=  0.0000, dx2-y2=  0.0000, dxy=  0.0000, 
+                 polarization =   0.3010, s =  0.0407, p =  0.2603, d =  0.0000, 
+     Spilling Parameter:   0.0032
+ 
+     PROJWFC      :  4m27.51s CPU     5m45.70s WALL
+
+ 
+   This run was terminated on:  20:18:36  26Mar2018            
+
+=------------------------------------------------------------------------------=
+   JOB DONE.
+=------------------------------------------------------------------------------=
+Job terminated normally
diff --git a/tests/parsers/test_projwfc.py b/tests/parsers/test_projwfc.py
index 6f30aa003..9f6afa18a 100644
--- a/tests/parsers/test_projwfc.py
+++ b/tests/parsers/test_projwfc.py
@@ -3,10 +3,24 @@
 """Tests for the `ProjwfcParser`."""
 from __future__ import absolute_import
 
+import os
 import pytest
 
 from aiida import orm
 from aiida.common import AttributeDict, LinkType
+from aiida.plugins import CalculationFactory
+
+from aiida_quantumespresso.parsers.parse_raw.projwfc import parse_lowdin_charges
+
+
+@pytest.mark.parametrize('test_file', ('default', 'spin'))
+def test_parse_lowdin_charges(test_file, parser_fixture_path, data_regression):
+    """Test parsing of lowdin charges from projwfc.x, for non-spin/spin cases."""
+    path = os.path.join(parser_fixture_path, 'projwfc', test_file, 'aiida.out')
+    with open(path) as handle:
+        data, spill_param = parse_lowdin_charges(handle.read().splitlines())
+    data['spill'] = spill_param
+    data_regression.check(data)
 
 
 @pytest.fixture
@@ -24,34 +38,53 @@ def projwfc_inputs(generate_calc_job_node, fixture_computer_localhost, generate_
     params = orm.Dict(dict={'number_of_spin_components': 1})
     params.add_incoming(parent_calcjob, link_type=LinkType.CREATE, link_label='output_parameters')
     params.store()
+    parameters = {'PROJWFC': {'emin': -1, 'emax': 1, 'DeltaE': 0.01, 'ngauss': 0, 'degauss': 0.01}}
     inputs = {
+        'parameters': orm.Dict(dict=parameters),
         'parent_folder': parent_calcjob.outputs.remote_folder,
     }
 
     return AttributeDict(inputs)
 
 
+def test_pw_link_spec():
+    """Test the ``PwCalculation`` input/output link names are as required for ``ProjwfcParser``.
+
+    ``ProjwfcParser`` relies on extracting data from the parent ``PwCalculation``.
+    This test safeguards against changes in the link names not being propagated to this parser.
+    """
+    pw_calc = CalculationFactory('quantumespresso.pw')
+    pw_spec = pw_calc.spec()
+    assert 'structure' in pw_spec.inputs, list(pw_spec.inputs.keys())
+    assert 'kpoints' in pw_spec.inputs, list(pw_spec.inputs.keys())
+    assert 'output_parameters' in pw_spec.outputs, list(pw_spec.outputs.keys())
+    assert 'output_structure' in pw_spec.outputs, list(pw_spec.outputs.keys())
+
+
 def test_projwfc_default(
-    fixture_database, fixture_computer_localhost, generate_calc_job_node, generate_parser, projwfc_inputs,
-    data_regression
+    fixture_database, fixture_computer_localhost, generate_calc_job_node, generate_parser, parse_from_node,
+    parser_fixture_path, projwfc_inputs, data_regression
 ):
     """Test ``ProjwfcParser`` on the results of a simple ``projwfc.x`` calculation."""
     entry_point_calc_job = 'quantumespresso.projwfc'
     entry_point_parser = 'quantumespresso.projwfc'
 
-    node = generate_calc_job_node(entry_point_calc_job, fixture_computer_localhost, 'default', projwfc_inputs)
+    node = generate_calc_job_node(entry_point_calc_job, fixture_computer_localhost, 'default', inputs=projwfc_inputs)
     parser = generate_parser(entry_point_parser)
-    results, calcfunction = parser.parse_from_node(node, store_provenance=False)
+    results, calcfunction = parse_from_node(
+        parser, node, store_provenance=False, retrieved_temp=os.path.join(parser_fixture_path, 'projwfc', 'default')
+    )
 
     assert calcfunction.is_finished, calcfunction.exception
     assert calcfunction.is_finished_ok, calcfunction.exit_message
 
-    for link_name in ['output_parameters', 'Dos', 'bands', 'projections']:
+    for link_name in ['output_parameters', 'Dos', 'bands', 'projections', 'lowdin']:
         assert link_name in results, list(results.keys())
 
     data_regression.check({
         'Dos': results['Dos'].attributes,
         'bands': results['bands'].attributes,
         'projections':
-        {k: v for k, v in results['projections'].attributes.items() if k not in ['reference_bandsdata_uuid']}
+        {k: v for k, v in results['projections'].attributes.items() if k not in ['reference_bandsdata_uuid']},
+        'lowdin': {k: v for k, v in results['lowdin'].attributes.items() if k not in ['structure_uuid']}
     })
diff --git a/tests/parsers/test_projwfc/test_parse_lowdin_charges_default_.yml b/tests/parsers/test_projwfc/test_parse_lowdin_charges_default_.yml
new file mode 100644
index 000000000..a9c05169a
--- /dev/null
+++ b/tests/parsers/test_projwfc/test_parse_lowdin_charges_default_.yml
@@ -0,0 +1,19 @@
+1:
+  sum:
+    total_charge: 3.9663
+  total_charge:
+    p: 2.8328
+    px: 0.9443
+    py: 0.9443
+    pz: 0.9443
+    s: 1.1336
+2:
+  sum:
+    total_charge: 3.9663
+  total_charge:
+    p: 2.8328
+    px: 0.9443
+    py: 0.9443
+    pz: 0.9443
+    s: 1.1336
+spill: 0.0084
diff --git a/tests/parsers/test_projwfc/test_parse_lowdin_charges_spin_.yml b/tests/parsers/test_projwfc/test_parse_lowdin_charges_spin_.yml
new file mode 100644
index 000000000..2168582fc
--- /dev/null
+++ b/tests/parsers/test_projwfc/test_parse_lowdin_charges_spin_.yml
@@ -0,0 +1,305 @@
+1:
+  polarization:
+    d: 3.5789
+    p: 0.0483
+    s: 0.0666
+  spin_down:
+    d: 1.3907
+    dx2-y2: 0.1557
+    dxy: 0.359
+    dxz: 0.3592
+    dyz: 0.3592
+    dz2: 0.1575
+    p: 3.5208
+    px: 1.1734
+    py: 1.1734
+    pz: 1.174
+    s: 1.184
+  spin_up:
+    d: 4.9696
+    dx2-y2: 0.9899
+    dxy: 0.9968
+    dxz: 0.9968
+    dyz: 0.9968
+    dz2: 0.9893
+    p: 9.7893
+    px: 1.1895
+    py: 1.1895
+    pz: 1.19
+    s: 1.2506
+  sum:
+    polarization: 3.6938
+    spin_down: 6.0955
+    spin_up: 9.7893
+    total_charge: 15.8849
+  total_charge:
+    d: 6.3602
+    p: 7.09
+    s: 2.4347
+2:
+  polarization:
+    d: 3.5779
+    p: 0.0483
+    s: 0.0666
+  spin_down:
+    d: 1.3919
+    dx2-y2: 0.1557
+    dxy: 0.3596
+    dxz: 0.3596
+    dyz: 0.3596
+    dz2: 0.1575
+    p: 3.5208
+    px: 1.1734
+    py: 1.1734
+    pz: 1.174
+    s: 1.184
+  spin_up:
+    d: 4.9699
+    dx2-y2: 0.9901
+    dxy: 0.9968
+    dxz: 0.9968
+    dyz: 0.9968
+    dz2: 0.9895
+    p: 9.7896
+    px: 1.1895
+    py: 1.1895
+    pz: 1.19
+    s: 1.2506
+  sum:
+    polarization: 3.6928
+    spin_down: 6.0968
+    spin_up: 9.7896
+    total_charge: 15.8864
+  total_charge:
+    d: 6.3618
+    p: 7.09
+    s: 2.4347
+3:
+  polarization:
+    d: 3.5777
+    p: 0.0483
+    s: 0.0666
+  spin_down:
+    d: 1.392
+    dx2-y2: 0.1557
+    dxy: 0.3589
+    dxz: 0.3597
+    dyz: 0.3601
+    dz2: 0.1575
+    p: 3.5208
+    px: 1.1734
+    py: 1.1734
+    pz: 1.174
+    s: 1.184
+  spin_up:
+    d: 4.9697
+    dx2-y2: 0.99
+    dxy: 0.9968
+    dxz: 0.9968
+    dyz: 0.9968
+    dz2: 0.9893
+    p: 9.7894
+    px: 1.1895
+    py: 1.1895
+    pz: 1.19
+    s: 1.2506
+  sum:
+    polarization: 3.6926
+    spin_down: 6.0968
+    spin_up: 9.7894
+    total_charge: 15.8863
+  total_charge:
+    d: 6.3617
+    p: 7.09
+    s: 2.4347
+4:
+  polarization:
+    d: 3.5777
+    p: 0.0483
+    s: 0.0666
+  spin_down:
+    d: 1.392
+    dx2-y2: 0.1557
+    dxy: 0.3589
+    dxz: 0.3601
+    dyz: 0.3597
+    dz2: 0.1575
+    p: 3.5208
+    px: 1.1734
+    py: 1.1734
+    pz: 1.174
+    s: 1.184
+  spin_up:
+    d: 4.9697
+    dx2-y2: 0.99
+    dxy: 0.9968
+    dxz: 0.9968
+    dyz: 0.9968
+    dz2: 0.9893
+    p: 9.7894
+    px: 1.1895
+    py: 1.1895
+    pz: 1.19
+    s: 1.2506
+  sum:
+    polarization: 3.6926
+    spin_down: 6.0968
+    spin_up: 9.7894
+    total_charge: 15.8863
+  total_charge:
+    d: 6.3617
+    p: 7.09
+    s: 2.4347
+5:
+  polarization:
+    d: 0.0
+    p: 0.2604
+    s: 0.0407
+  spin_down:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 2.1274
+    px: 0.709
+    py: 0.709
+    pz: 0.7093
+    s: 0.7441
+  spin_up:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 3.1726
+    px: 0.7964
+    py: 0.7964
+    pz: 0.795
+    s: 0.7848
+  sum:
+    polarization: 0.3011
+    spin_down: 2.8715
+    spin_up: 3.1726
+    total_charge: 6.0441
+  total_charge:
+    d: 0.0
+    p: 4.5151
+    s: 1.529
+6:
+  polarization:
+    d: 0.0
+    p: 0.2603
+    s: 0.0407
+  spin_down:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 2.1274
+    px: 0.709
+    py: 0.709
+    pz: 0.7093
+    s: 0.7441
+  spin_up:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 3.1725
+    px: 0.7963
+    py: 0.7967
+    pz: 0.7947
+    s: 0.7848
+  sum:
+    polarization: 0.301
+    spin_down: 2.8715
+    spin_up: 3.1725
+    total_charge: 6.044
+  total_charge:
+    d: 0.0
+    p: 4.515
+    s: 1.529
+7:
+  polarization:
+    d: 0.0
+    p: 0.2603
+    s: 0.0407
+  spin_down:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 2.1274
+    px: 0.709
+    py: 0.709
+    pz: 0.7093
+    s: 0.7441
+  spin_up:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 3.1725
+    px: 0.7967
+    py: 0.7963
+    pz: 0.7947
+    s: 0.7848
+  sum:
+    polarization: 0.301
+    spin_down: 2.8715
+    spin_up: 3.1725
+    total_charge: 6.044
+  total_charge:
+    d: 0.0
+    p: 4.515
+    s: 1.529
+8:
+  polarization:
+    d: 0.0
+    p: 0.2603
+    s: 0.0407
+  spin_down:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 2.1274
+    px: 0.709
+    py: 0.709
+    pz: 0.7093
+    s: 0.7441
+  spin_up:
+    d: 0.0
+    dx2-y2: 0.0
+    dxy: 0.0
+    dxz: 0.0
+    dyz: 0.0
+    dz2: 0.0
+    p: 3.1725
+    px: 0.7965
+    py: 0.7965
+    pz: 0.7946
+    s: 0.7848
+  sum:
+    polarization: 0.301
+    spin_down: 2.8715
+    spin_up: 3.1725
+    total_charge: 6.044
+  total_charge:
+    d: 0.0
+    p: 4.515
+    s: 1.529
+spill: 0.0032
diff --git a/tests/parsers/test_projwfc/test_projwfc_default.yml b/tests/parsers/test_projwfc/test_projwfc_default.yml
index 3dd573efd..ebf3e4794 100644
--- a/tests/parsers/test_projwfc/test_projwfc_default.yml
+++ b/tests/parsers/test_projwfc/test_projwfc_default.yml
@@ -17,6 +17,25 @@ bands:
   - 36
   - 3
   units: eV
+lowdin:
+  site_data:
+  - sum:
+      total_charge: 3.9663
+    total_charge:
+      p: 2.8328
+      px: 0.9443
+      py: 0.9443
+      pz: 0.9443
+      s: 1.1336
+  - sum:
+      total_charge: 3.9663
+    total_charge:
+      p: 2.8328
+      px: 0.9443
+      py: 0.9443
+      pz: 0.9443
+      s: 1.1336
+  spill_parameter: 0.0084
 projections:
   array|energy_array_0:
   - 22