diff --git a/python/solid_dmft/csc_flow.py b/python/solid_dmft/csc_flow.py
index 9d1a8dfa..7527921a 100644
--- a/python/solid_dmft/csc_flow.py
+++ b/python/solid_dmft/csc_flow.py
@@ -29,6 +29,7 @@
 import subprocess
 import shlex
 import os
+import gc
 import numpy as np
 
 # triqs
@@ -89,8 +90,8 @@ def _run_w90converter(seedname, tolerance):
 
     # Checks if creating of rot_mat succeeded
     if mpi.is_master_node():
-        with HDFArchive(seedname+'.h5', 'r') as archive:
-            assert archive['dft_input']['use_rotations'], 'Creation of rot_mat failed in W90 converter'
+        with HDFArchive(seedname+'.h5', 'r') as ar:
+            assert ar['dft_input']['use_rotations'], 'Creation of rot_mat failed in W90 converter'
     mpi.barrier()
 
 def _full_qe_run(seedname, dft_params, mode):
@@ -248,7 +249,8 @@ def csc_flow_control(general_params, solver_params, dft_params, gw_params, advan
     """
 
     # Removes legacy file vasp.suppress_projs if present
-    vasp.remove_legacy_projections_suppressed()
+    if dft_params['dft_code'] == 'vasp':
+        vasp.remove_legacy_projections_suppressed()
 
     # if GAMMA file already exists, load it by doing extra DFT iterations
     if dft_params['dft_code'] == 'vasp' and os.path.exists('GAMMA'):
@@ -259,9 +261,9 @@ def csc_flow_control(general_params, solver_params, dft_params, gw_params, advan
     # Reads in iteration offset if restarting
     iteration_offset = 0
     if mpi.is_master_node() and os.path.isfile(general_params['seedname']+'.h5'):
-        with HDFArchive(general_params['seedname']+'.h5', 'r') as archive:
-            if 'DMFT_results' in archive and 'iteration_count' in archive['DMFT_results']:
-                iteration_offset = archive['DMFT_results']['iteration_count']
+        with HDFArchive(general_params['seedname']+'.h5', 'r') as ar:
+            if 'DMFT_results' in ar and 'iteration_count' in ar['DMFT_results']:
+                iteration_offset = ar['DMFT_results']['iteration_count']
     iteration_offset = mpi.bcast(iteration_offset)
 
     iter_dmft = iteration_offset+1
@@ -359,6 +361,9 @@ def csc_flow_control(general_params, solver_params, dft_params, gw_params, advan
         end_time_dft = timer()
         mpi.report('  solid_dmft: DFT cycle took {:10.3f} seconds'.format(end_time_dft-start_time_dft))
 
+        del sum_k
+        gc.collect()
+
     # Kills background VASP process for clean end
     if mpi.is_master_node() and dft_params['dft_code'] == 'vasp':
         print('  solid_dmft: Stopping VASP\n', flush=True)
diff --git a/python/solid_dmft/dmft_cycle.py b/python/solid_dmft/dmft_cycle.py
index d3ddbc15..0c53a20f 100755
--- a/python/solid_dmft/dmft_cycle.py
+++ b/python/solid_dmft/dmft_cycle.py
@@ -27,6 +27,7 @@
 
 
 # system
+import gc
 from copy import deepcopy
 from timeit import default_timer as timer
 import numpy as np
@@ -143,7 +144,7 @@ def _determine_block_structure(sum_k, general_params, advanced_params, solver_ty
     # Applies manual selection of the solver struct
     if any(s is not None for s in advanced_params['pick_solver_struct']):
         mpi.report('selecting subset of orbital space for gf_struct_solver from input:')
-        mpi.report(advanced_params['pick_solver_struct'][icrsh])
+        mpi.report(advanced_params['pick_solver_struct'])
         sum_k.block_structure.pick_gf_struct_solver(advanced_params['pick_solver_struct'])
 
     # Applies the manual mapping to each inequivalent shell
@@ -218,7 +219,8 @@ def _calculate_rotation_matrix(general_params, sum_k):
     # in case sum_k.use_rotations == False before:
     sum_k.use_rotations = True
     # sum_k.eff_atomic_levels() needs to be recomputed if rot_mat were changed
-    if hasattr(sum_k, "Hsumk"): delattr(sum_k, "Hsumk")
+    if hasattr(sum_k, "Hsumk"):
+        delattr(sum_k, "Hsumk")
     mpi.report('Updating rotation matrices using dft {} eigenbasis to maximise sign'.format(general_params['set_rot']))
 
     # Prints matrices
@@ -313,31 +315,30 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
     iteration_offset = 0
 
     # determine chemical potential for bare DFT sum_k object
+    archive = general_params['jobname']+'/'+general_params['seedname']+'.h5'
     if mpi.is_master_node():
-        archive = HDFArchive(general_params['jobname']+'/'+general_params['seedname']+'.h5', 'a')
-        if 'DMFT_results' not in archive:
-            archive.create_group('DMFT_results')
-        if 'last_iter' not in archive['DMFT_results']:
-            archive['DMFT_results'].create_group('last_iter')
-        if 'DMFT_input' not in archive:
-            archive.create_group('DMFT_input')
-            archive['DMFT_input']['program'] = 'solid_dmft'
-            archive['DMFT_input'].create_group('solver')
-            archive['DMFT_input'].create_group('version')
-            archive['DMFT_input']['version']['triqs_hash'] = triqs_hash
-            archive['DMFT_input']['version']['triqs_version'] = triqs_version
-            archive['DMFT_input']['version']['solid_dmft_hash'] = solid_dmft_hash
-            archive['DMFT_input']['version']['solid_dmft_version'] = solid_dmft_version
-
-        if 'iteration_count' in archive['DMFT_results']:
-            iteration_offset = archive['DMFT_results/iteration_count']
-            sum_k.chemical_potential = archive['DMFT_results/last_iter/chemical_potential_post']
-            print(f'RESTARTING DMFT RUN at iteration {iteration_offset+1} using last self-energy')
-        else:
-            print('INITIAL DMFT RUN')
-        print('#'*80, '\n')
-    else:
-        archive = None
+        with HDFArchive(archive, 'a') as ar:
+            if 'DMFT_results' not in ar:
+                ar.create_group('DMFT_results')
+            if 'last_iter' not in ar['DMFT_results']:
+                ar['DMFT_results'].create_group('last_iter')
+            if 'DMFT_input' not in ar:
+                ar.create_group('DMFT_input')
+                ar['DMFT_input']['program'] = 'solid_dmft'
+                ar['DMFT_input'].create_group('solver')
+                ar['DMFT_input'].create_group('version')
+                ar['DMFT_input']['version']['triqs_hash'] = triqs_hash
+                ar['DMFT_input']['version']['triqs_version'] = triqs_version
+                ar['DMFT_input']['version']['solid_dmft_hash'] = solid_dmft_hash
+                ar['DMFT_input']['version']['solid_dmft_version'] = solid_dmft_version
+
+            if 'iteration_count' in ar['DMFT_results']:
+                iteration_offset = ar['DMFT_results/iteration_count']
+                sum_k.chemical_potential = ar['DMFT_results/last_iter/chemical_potential_post']
+                print(f'RESTARTING DMFT RUN at iteration {iteration_offset+1} using last self-energy')
+            else:
+                print('INITIAL DMFT RUN')
+            print('#'*80, '\n')
 
     iteration_offset = mpi.bcast(iteration_offset)
     sum_k.chemical_potential = mpi.bcast(sum_k.chemical_potential)
@@ -388,9 +389,10 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
 
     # check for previous broyden data oterhwise initialize it:
     if mpi.is_master_node() and  general_params['g0_mix_type'] == 'broyden':
-        if not 'broyler' in archive['DMFT_results']:
-            archive['DMFT_results']['broyler'] = [{'mu' : [],'V': [], 'dV': [], 'F': [], 'dF': []}
-                                                  for _ in range(sum_k.n_inequiv_shells)]
+        with HDFArchive(archive, 'a') as ar:
+            if 'broyler' not in ar['DMFT_results']:
+                ar['DMFT_results']['broyler'] = [{'mu' : [],'V': [], 'dV': [], 'F': [], 'dF': []}
+                                                      for _ in range(sum_k.n_inequiv_shells)]
 
     # Generates a rotation matrix to change the basis
     if general_params['set_rot'] is not None:
@@ -398,7 +400,8 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
         sum_k = _calculate_rotation_matrix(general_params, sum_k)
     # Saves rotation matrix to h5 archive:
     if mpi.is_master_node() and iteration_offset == 0:
-        archive['DMFT_input']['rot_mat'] = sum_k.rot_mat
+        with HDFArchive(archive, 'a') as ar:
+            ar['DMFT_input']['rot_mat'] = sum_k.rot_mat
     mpi.barrier()
 
     # Calculates density in default block structure
@@ -413,7 +416,8 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
     det_blocks = None
     # load previous block_structure if possible
     if mpi.is_master_node():
-        det_blocks = 'block_structure' not in archive['DMFT_input']
+        with HDFArchive(archive, 'a') as ar:
+            det_blocks = 'block_structure' not in ar['DMFT_input']
     det_blocks = mpi.bcast(det_blocks)
 
     # Previous rot_mat only not None if the rot_mat changed from load_sigma or previous run
@@ -449,16 +453,17 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
     else:
         # Master node checks if rot_mat stayed the same
         if mpi.is_master_node():
-            sum_k.block_structure = archive['DMFT_input']['block_structure']
-            sum_k.deg_shells = archive['DMFT_input/deg_shells']
-            previous_rot_mat = archive['DMFT_input']['rot_mat']
-            if not all(np.allclose(x, y) for x, y in zip(sum_k.rot_mat, previous_rot_mat)):
-                print('WARNING: rot_mat in current step is different from previous step.')
-                archive['DMFT_input']['rot_mat'] = sum_k.rot_mat
-            else:
-                previous_rot_mat = None
-            if 'solver_struct_ftps' in archive['DMFT_input']:
-                deg_orbs_ftps = archive['DMFT_input/solver_struct_ftps']
+            with HDFArchive(archive, 'a') as ar:
+                sum_k.block_structure = ar['DMFT_input']['block_structure']
+                sum_k.deg_shells = ar['DMFT_input/deg_shells']
+                previous_rot_mat = ar['DMFT_input']['rot_mat']
+                if not all(np.allclose(x, y) for x, y in zip(sum_k.rot_mat, previous_rot_mat)):
+                    print('WARNING: rot_mat in current step is different from previous step.')
+                    ar['DMFT_input']['rot_mat'] = sum_k.rot_mat
+                else:
+                    previous_rot_mat = None
+                if 'solver_struct_ftps' in ar['DMFT_input']:
+                    deg_orbs_ftps = ar['DMFT_input/solver_struct_ftps']
 
         sum_k.block_structure = mpi.bcast(sum_k.block_structure)
         sum_k.deg_shells = mpi.bcast(sum_k.deg_shells)
@@ -469,7 +474,8 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
     # Sumk doesn't hold corr_to_inequiv anymore, which is in block_structure now
     if sum_k.block_structure.corr_to_inequiv is None:
         if mpi.is_master_node():
-            sum_k.block_structure.corr_to_inequiv = archive['dft_input/corr_to_inequiv']
+            with HDFArchive(archive, 'r') as ar:
+                sum_k.block_structure.corr_to_inequiv = ar['dft_input/corr_to_inequiv']
         sum_k.block_structure = mpi.bcast(sum_k.block_structure)
 
     # Determination of shell_multiplicity
@@ -486,21 +492,23 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
     # Constructs interaction Hamiltonian and writes it to the h5 archive
     h_int, gw_params = interaction_hamiltonian.construct(sum_k, general_params, solver_type_per_imp, gw_params)
     if mpi.is_master_node():
-        archive['DMFT_input']['h_int'] = h_int
+        with HDFArchive(archive, 'a') as ar:
+            ar['DMFT_input']['h_int'] = h_int
 
     # If new calculation, writes input parameters and sum_k <-> solver mapping to archive
     if iteration_offset == 0:
         if mpi.is_master_node():
-            archive['DMFT_input']['general_params'] = dict_to_h5.prep_params_for_h5(general_params)
-            archive['DMFT_input']['solver_params'] = dict_to_h5.prep_params_for_h5(solver_params)
-            archive['DMFT_input']['dft_params'] = dict_to_h5.prep_params_for_h5(dft_params)
-            archive['DMFT_input']['advanced_params'] = dict_to_h5.prep_params_for_h5(advanced_params)
-
-            archive['DMFT_input']['block_structure'] = sum_k.block_structure
-            archive['DMFT_input']['deg_shells'] = sum_k.deg_shells
-            archive['DMFT_input']['shell_multiplicity'] = shell_multiplicity
-            if deg_orbs_ftps is not None:
-                archive['DMFT_input']['solver_struct_ftps'] = deg_orbs_ftps
+            with HDFArchive(archive, 'a') as ar:
+                ar['DMFT_input']['general_params'] = dict_to_h5.prep_params_for_h5(general_params)
+                ar['DMFT_input']['solver_params'] = dict_to_h5.prep_params_for_h5(solver_params)
+                ar['DMFT_input']['dft_params'] = dict_to_h5.prep_params_for_h5(dft_params)
+                ar['DMFT_input']['advanced_params'] = dict_to_h5.prep_params_for_h5(advanced_params)
+
+                ar['DMFT_input']['block_structure'] = sum_k.block_structure
+                ar['DMFT_input']['deg_shells'] = sum_k.deg_shells
+                ar['DMFT_input']['shell_multiplicity'] = shell_multiplicity
+                if deg_orbs_ftps is not None:
+                    ar['DMFT_input']['solver_struct_ftps'] = deg_orbs_ftps
     mpi.barrier()
 
     solvers = [None] * sum_k.n_inequiv_shells
@@ -520,14 +528,15 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
 
     # store solver hash to archive
     if mpi.is_master_node():
-        if 'version' not in archive['DMFT_input']:
-            archive['DMFT_input'].create_group('version')
-        for iineq, solver in enumerate(solvers):
-            if f'solver {iineq}' not in archive['DMFT_input']['version']:
-                archive['DMFT_input']['version'].create_group(f'solver {iineq}')
-            archive['DMFT_input']['version'][f'solver {iineq}']['name'] = solver.solver_params['type']
-            archive['DMFT_input']['version'][f'solver {iineq}']['hash'] = solver.git_hash
-            archive['DMFT_input']['version'][f'solver {iineq}']['version'] = solver.version
+        with HDFArchive(archive, 'a') as ar:
+            if 'version' not in ar['DMFT_input']:
+                ar['DMFT_input'].create_group('version')
+            for iineq, solver in enumerate(solvers):
+                if f'solver {iineq}' not in ar['DMFT_input']['version']:
+                    ar['DMFT_input']['version'].create_group(f'solver {iineq}')
+                ar['DMFT_input']['version'][f'solver {iineq}']['name'] = solver.solver_params['type']
+                ar['DMFT_input']['version'][f'solver {iineq}']['hash'] = solver.git_hash
+                ar['DMFT_input']['version'][f'solver {iineq}']['version'] = solver.version
 
     # Extracts local GF per *inequivalent* shell
     local_gf_dft = sum_k.extract_G_loc(broadening=broadening, with_Sigma=False, mu=dft_mu)
@@ -546,21 +555,19 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
     mpi.report('\n {} DMFT cycles requested. Starting with iteration  {}.\n'.format(n_iter, iteration_offset+1))
 
     # Prepares observable and conv dicts
-    observables = None
-    conv_obs = None
     if mpi.is_master_node():
         observables = prep_observables(archive, sum_k)
         conv_obs = convergence.prep_conv_obs(archive)
-    observables = mpi.bcast(observables)
-    conv_obs = mpi.bcast(conv_obs)
-
-    if mpi.is_master_node() and iteration_offset == 0:
-        write_header_to_file(general_params, sum_k)
-        observables = add_dft_values_as_zeroth_iteration(observables, general_params, solver_type_per_imp, dft_mu, dft_energy, sum_k,
-                                                         local_gf_dft, shell_multiplicity)
-        write_obs(observables, sum_k, general_params)
-        # write convergence file
-        convergence.prep_conv_file(general_params, sum_k)
+        if iteration_offset == 0:
+            write_header_to_file(general_params, sum_k)
+            observables = add_dft_values_as_zeroth_iteration(observables, general_params, solver_type_per_imp, dft_mu, dft_energy, sum_k,
+                                                             local_gf_dft, shell_multiplicity)
+            write_obs(observables, sum_k, general_params)
+            # write convergence file
+            convergence.prep_conv_file(general_params, sum_k)
+    else:
+        observables = None
+        conv_obs = None
 
     # The infamous DMFT self consistency cycle
     is_converged = False
@@ -581,6 +588,7 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
                                                  h_int, archive, shell_multiplicity, E_kin_dft,
                                                  observables, conv_obs, ops_chi_measure, dft_irred_kpt_indices, dft_energy, broadening,
                                                  is_converged, is_sampling=False)
+        gc.collect()
 
         if is_converged:
             break
@@ -611,10 +619,6 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
 
     mpi.barrier()
 
-    # close the h5 archive
-    if mpi.is_master_node():
-        del archive
-
     return is_converged, sum_k
 
 
@@ -635,8 +639,6 @@ def _dmft_step(sum_k, solvers, it, general_params, solver_params, gw_params,
     density_shell_pre = np.zeros(sum_k.n_inequiv_shells)
     density_mat_pre = [None] * sum_k.n_inequiv_shells
 
-    mpi.barrier()
-
     if sum_k.SO:
         printed = ((np.real, 'real'), (np.imag, 'imaginary'))
     else:
@@ -716,18 +718,18 @@ def _dmft_step(sum_k, solvers, it, general_params, solver_params, gw_params,
                                 Hloc_0 += spin_block[o1,o2].real/2 * (c_dag(spin,o1) * c(spin,o2) + c_dag(spin,o2) * c(spin,o1))
                 solvers[icrsh].Hloc_0 = Hloc_0
 
-         # store solver to h5 archive
-        if general_params['store_solver'] and mpi.is_master_node():
-            if 'solver' not in archive['DMFT_input']:
-                archive['DMFT_input'].create_group('solver')
-            archive['DMFT_input/solver'].create_group('it_'+str(it))
-            archive['DMFT_input/solver/it_'+str(it)]['S_'+str(icrsh)] = solvers[icrsh].triqs_solver
-
         # store DMFT input directly in last_iter
         if mpi.is_master_node():
-            archive['DMFT_results/last_iter']['G0_freq_{}'.format(icrsh)] = solvers[icrsh].G0_freq
-            if solver_type_per_imp[icrsh] == 'cthyb' and solvers[icrsh].solver_params['delta_interface']:
-                archive['DMFT_results/last_iter']['Delta_time_{}'.format(icrsh)] = solvers[icrsh].Delta_time
+            with HDFArchive(archive, 'a') as ar:
+                ar['DMFT_results/last_iter']['G0_freq_{}'.format(icrsh)] = solvers[icrsh].G0_freq
+                if solver_type_per_imp[icrsh] == 'cthyb' and solvers[icrsh].solver_params['delta_interface']:
+                    ar['DMFT_results/last_iter']['Delta_time_{}'.format(icrsh)] = solvers[icrsh].Delta_time
+                # store solver to h5 archive
+                if general_params['store_solver']:
+                    if 'solver' not in ar['DMFT_input']:
+                        ar['DMFT_input'].create_group('solver')
+                    ar['DMFT_input/solver'].create_group('it_'+str(it))
+                    ar['DMFT_input/solver/it_'+str(it)]['S_'+str(icrsh)] = solvers[icrsh].triqs_solver
 
         # setup of measurement of chi(SzSz(tau) if requested
         # TODO: move this into solver class?
@@ -768,7 +770,8 @@ def _dmft_step(sum_k, solvers, it, general_params, solver_params, gw_params,
 
         # update solver in h5 archive
         if general_params['store_solver'] and mpi.is_master_node():
-            archive['DMFT_input/solver/it_'+str(it)]['S_'+str(icrsh)] = solvers[icrsh].triqs_solver
+            with HDFArchive(archive, 'a') as ar:
+                ar['DMFT_input/solver/it_'+str(it)]['S_'+str(icrsh)] = solvers[icrsh].triqs_solver
 
     # Done with loop over impurities
 
@@ -841,35 +844,35 @@ def _dmft_step(sum_k, solvers, it, general_params, solver_params, gw_params,
 
         write_obs(observables, sum_k, general_params)
 
-        # write the new observable array to h5 archive
-        archive['DMFT_results']['observables'] = observables
-
-    # Computes convergence quantities and writes them to file
-    if mpi.is_master_node():
+        # Computes convergence quantities and writes them to file
         conv_obs = convergence.calc_convergence_quantities(sum_k, general_params, conv_obs, observables,
                                                            solvers, G0_freq_previous, G_loc_all, Sigma_freq_previous)
         convergence.write_conv(conv_obs, sum_k, general_params)
-        archive['DMFT_results']['convergence_obs'] = conv_obs
-    conv_obs = mpi.bcast(conv_obs)
 
-    mpi.report('*** iteration finished ***')
+        with HDFArchive(archive, 'a') as ar:
+            ar['DMFT_results']['observables'] = observables
+            ar['DMFT_results']['convergence_obs'] = conv_obs
 
-    # Checks for convergence
-    is_now_converged = convergence.check_convergence(sum_k.n_inequiv_shells, general_params, conv_obs)
-    if is_now_converged is None:
-        is_converged = False
-    else:
-        # if convergency criteria was already reached don't overwrite it!
-        is_converged = is_converged or is_now_converged
-
-    # Final prints
-    formatter.print_summary_observables(observables, sum_k.n_inequiv_shells,
-                                        sum_k.spin_block_names[sum_k.SO])
-    if general_params['calc_energies']:
-        formatter.print_summary_energetics(observables)
-    if general_params['magnetic'] and sum_k.SO == 0:
-        # if a magnetic calculation is done print out a summary of up/down occ
-        formatter.print_summary_magnetic_occ(observables, sum_k.n_inequiv_shells)
-    formatter.print_summary_convergence(conv_obs, general_params, sum_k.n_inequiv_shells)
+        mpi.report('*** iteration finished ***')
+
+        # Checks for convergence
+        is_now_converged = convergence.check_convergence(sum_k.n_inequiv_shells, general_params, conv_obs)
+        if is_now_converged is None:
+            is_converged = False
+        else:
+            # if convergency criteria was already reached don't overwrite it!
+            is_converged = is_converged or is_now_converged
+
+        # Final prints
+        formatter.print_summary_observables(observables, sum_k.n_inequiv_shells,
+                                            sum_k.spin_block_names[sum_k.SO])
+        if general_params['calc_energies']:
+            formatter.print_summary_energetics(observables)
+        if general_params['magnetic'] and sum_k.SO == 0:
+            # if a magnetic calculation is done print out a summary of up/down occ
+            formatter.print_summary_magnetic_occ(observables, sum_k.n_inequiv_shells)
+        formatter.print_summary_convergence(conv_obs, general_params, sum_k.n_inequiv_shells)
+
+    mpi.bcast(is_converged)
 
     return sum_k, solvers, observables, is_converged
diff --git a/python/solid_dmft/dmft_tools/afm_mapping.py b/python/solid_dmft/dmft_tools/afm_mapping.py
index 5aa93bff..150c1be5 100644
--- a/python/solid_dmft/dmft_tools/afm_mapping.py
+++ b/python/solid_dmft/dmft_tools/afm_mapping.py
@@ -24,6 +24,7 @@
 ################################################################################
 
 import numpy as np
+from h5 import HDFArchive
 import triqs.utility.mpi as mpi
 
 def determine(general_params, archive, n_inequiv_shells):
@@ -36,42 +37,43 @@ def determine(general_params, archive, n_inequiv_shells):
     afm_mapping = None
     if mpi.is_master_node():
         # Reads mapping from h5 archive if it exists already from a previous run
-        if 'afm_mapping' in archive['DMFT_input']:
-            afm_mapping = archive['DMFT_input']['afm_mapping']
-        elif len(general_params['magmom']) == n_inequiv_shells:
-            # find equal or opposite spin imps, where we use the magmom array to
-            # identity those with equal numbers or opposite
-            # [copy Yes/False, from where, switch up/down channel]
-            afm_mapping = [None] * n_inequiv_shells
-            abs_moms = np.abs(general_params['magmom'])
-
-            for icrsh in range(n_inequiv_shells):
-                # if the moment was seen before ...
-                previous_occurences = np.nonzero(np.isclose(abs_moms[:icrsh], abs_moms[icrsh]))[0]
-                if previous_occurences.size > 0:
-                    # find the source imp to copy from
-                    source = np.min(previous_occurences)
-                    # determine if we need to switch up and down channel
-                    switch = np.isclose(general_params['magmom'][icrsh], -general_params['magmom'][source])
-
-                    afm_mapping[icrsh] = [True, source, switch]
-                else:
-                    afm_mapping[icrsh] = [False, icrsh, False]
-
-
-            print('AFM calculation selected, mapping self energies as follows:')
-            print('imp  [copy sigma, source imp, switch up/down]')
-            print('---------------------------------------------')
-            for i, elem in enumerate(afm_mapping):
-                print('{}: {}'.format(i, elem))
-            print('')
-
-            archive['DMFT_input']['afm_mapping'] = afm_mapping
-
-        # if anything did not work set afm_order false
-        else:
-            print('WARNING: couldn\'t determine afm mapping. No mapping used.')
-            general_params['afm_order'] = False
+        with HDFArchive(archive, 'a') as ar:
+            if 'afm_mapping' in ar['DMFT_input']:
+                afm_mapping = ar['DMFT_input']['afm_mapping']
+            elif len(general_params['magmom']) == n_inequiv_shells:
+                # find equal or opposite spin imps, where we use the magmom array to
+                # identity those with equal numbers or opposite
+                # [copy Yes/False, from where, switch up/down channel]
+                afm_mapping = [None] * n_inequiv_shells
+                abs_moms = np.abs(general_params['magmom'])
+
+                for icrsh in range(n_inequiv_shells):
+                    # if the moment was seen before ...
+                    previous_occurences = np.nonzero(np.isclose(abs_moms[:icrsh], abs_moms[icrsh]))[0]
+                    if previous_occurences.size > 0:
+                        # find the source imp to copy from
+                        source = np.min(previous_occurences)
+                        # determine if we need to switch up and down channel
+                        switch = np.isclose(general_params['magmom'][icrsh], -general_params['magmom'][source])
+
+                        afm_mapping[icrsh] = [True, source, switch]
+                    else:
+                        afm_mapping[icrsh] = [False, icrsh, False]
+
+
+                print('AFM calculation selected, mapping self energies as follows:')
+                print('imp  [copy sigma, source imp, switch up/down]')
+                print('---------------------------------------------')
+                for i, elem in enumerate(afm_mapping):
+                    print('{}: {}'.format(i, elem))
+                print('')
+
+                ar['DMFT_input']['afm_mapping'] = afm_mapping
+
+            # if anything did not work set afm_order false
+            else:
+                print('WARNING: couldn\'t determine afm mapping. No mapping used.')
+                general_params['afm_order'] = False
 
     general_params['afm_order'] = mpi.bcast(general_params['afm_order'])
     if general_params['afm_order']:
diff --git a/python/solid_dmft/dmft_tools/convergence.py b/python/solid_dmft/dmft_tools/convergence.py
index 52116e62..d0f6951a 100644
--- a/python/solid_dmft/dmft_tools/convergence.py
+++ b/python/solid_dmft/dmft_tools/convergence.py
@@ -29,6 +29,7 @@
 import numpy as np
 
 # triqs
+from h5 import HDFArchive
 from triqs.gf import MeshImFreq, MeshImTime, MeshReFreq, BlockGf
 from solid_dmft.dmft_tools import common
 
@@ -183,13 +184,14 @@ def prep_conv_obs(h5_archive):
         conv array for calculation
     """
 
-    # determine number of impurities
-    n_inequiv_shells = h5_archive['dft_input']['n_inequiv_shells']
+    with HDFArchive(h5_archive, 'r') as ar:
+        # determine number of impurities
+        n_inequiv_shells = ar['dft_input']['n_inequiv_shells']
 
-    # check for previous iterations
-    conv_prev = []
-    if 'convergence_obs' in h5_archive['DMFT_results']:
-        conv_prev = h5_archive['DMFT_results']['convergence_obs']
+        # check for previous iterations
+        conv_prev = []
+        if 'convergence_obs' in ar['DMFT_results']:
+            conv_prev = ar['DMFT_results']['convergence_obs']
 
     # prepare observable dicts
     if len(conv_prev) > 0:
diff --git a/python/solid_dmft/dmft_tools/greens_functions_mixer.py b/python/solid_dmft/dmft_tools/greens_functions_mixer.py
index 2e5cf47d..dbf01144 100644
--- a/python/solid_dmft/dmft_tools/greens_functions_mixer.py
+++ b/python/solid_dmft/dmft_tools/greens_functions_mixer.py
@@ -28,6 +28,7 @@
 
 # triqs
 import triqs.utility.mpi as mpi
+from h5 import HDFArchive
 from triqs.gf import MeshImFreq, MeshImTime, Gf, make_gf_from_fourier
 from triqs.gf.descriptors import Fourier
 from triqs.gf.tools import inverse
@@ -194,14 +195,15 @@ def mix_g0(solver, general_params, icrsh, archive, G0_freq_previous, it, deg_she
         mpi.report('\n############\n!!!! WARNING !!!! broyden mixing is still in early testing stage ! Use with caution.\n############\n')
         # TODO implement broyden mixing for Sigma
         if mpi.is_master_node():
-            broyler = archive['DMFT_results']['broyler']
-            # calculate the next G0 via broyden scheme
-            G0_broyden_update, broyler[icrsh] = _broyden_update(it, broyler[icrsh], general_params,
-                                                                deg_shell, solver.G0_freq,
-                                                                G0_freq_previous)
-            # store broyden update to h5 archive
-            archive['DMFT_results']['broyler'] = broyler
-            solver.G0_freq << G0_broyden_update
+            with HDFArchive(archive, 'a') as ar:
+                broyler = ar['DMFT_results']['broyler']
+                # calculate the next G0 via broyden scheme
+                G0_broyden_update, broyler[icrsh] = _broyden_update(it, broyler[icrsh], general_params,
+                                                                    deg_shell, solver.G0_freq,
+                                                                    G0_freq_previous)
+                # store broyden update to h5 archive
+                ar['DMFT_results']['broyler'] = broyler
+                solver.G0_freq << G0_broyden_update
 
         solver.G0_freq << mpi.bcast(solver.G0_freq)
 
diff --git a/python/solid_dmft/dmft_tools/initial_self_energies.py b/python/solid_dmft/dmft_tools/initial_self_energies.py
index 320a6fcd..ff7b4718 100755
--- a/python/solid_dmft/dmft_tools/initial_self_energies.py
+++ b/python/solid_dmft/dmft_tools/initial_self_energies.py
@@ -299,18 +299,19 @@ def _load_sigma_from_h5(h5_archive, iteration):
     internal_path = 'DMFT_results/'
     internal_path += 'last_iter' if iteration == -1 else 'it_{}'.format(iteration)
 
-    n_inequiv_shells = h5_archive['dft_input']['n_inequiv_shells']
+    with HDFArchive(h5_archive, 'r') as ar:
+        n_inequiv_shells = ar['dft_input']['n_inequiv_shells']
 
-    # Loads previous self-energies and DC
-    self_energies = [h5_archive[internal_path]['Sigma_freq_{}'.format(iineq)]
-                     for iineq in range(n_inequiv_shells)]
-    last_g0 = [h5_archive[internal_path]['G0_freq_{}'.format(iineq)]
-                     for iineq in range(n_inequiv_shells)]
-    dc_imp = h5_archive[internal_path]['DC_pot']
-    dc_energy = h5_archive[internal_path]['DC_energ']
+        # Loads previous self-energies and DC
+        self_energies = [ar[internal_path]['Sigma_freq_{}'.format(iineq)]
+                         for iineq in range(n_inequiv_shells)]
+        last_g0 = [ar[internal_path]['G0_freq_{}'.format(iineq)]
+                         for iineq in range(n_inequiv_shells)]
+        dc_imp = ar[internal_path]['DC_pot']
+        dc_energy = ar[internal_path]['DC_energ']
 
-    # Loads density_matrix to recalculate DC if dc_dmft
-    density_matrix = h5_archive[internal_path]['dens_mat_post']
+        # Loads density_matrix to recalculate DC if dc_dmft
+        density_matrix = ar[internal_path]['dens_mat_post']
 
     print('Loaded Sigma_imp0...imp{} '.format(n_inequiv_shells-1)
           + ('at last it ' if iteration == -1 else 'at it {} '.format(iteration)))
@@ -495,9 +496,8 @@ def determine_dc_and_initial_sigma(general_params, gw_params, advanced_params, s
         # Loads Sigma from different calculation
         elif general_params['load_sigma']:
             print('\nFrom {}:'.format(general_params['path_to_sigma']), end=' ')
-            with HDFArchive(general_params['path_to_sigma'], 'r') as sigma_archive:
-                (loaded_sigma, loaded_dc_imp, _,
-                 _, loaded_density_matrix) = _load_sigma_from_h5(sigma_archive, general_params['load_sigma_iter'])
+            (loaded_sigma, loaded_dc_imp, _,
+             _, loaded_density_matrix) = _load_sigma_from_h5(general_params['path_to_sigma'], general_params['load_sigma_iter'])
 
             # Recalculate double counting in case U, J or DC formula changed
             if general_params['dc']:
@@ -538,6 +538,7 @@ def determine_dc_and_initial_sigma(general_params, gw_params, advanced_params, s
                             start_sigma[icrsh][spin_channel] <<  dc_pot[spin_channel] - fac*np.eye(n_orb)
                         else:
                             start_sigma[icrsh][spin_channel] << dc_pot[spin_channel] + fac*np.eye(n_orb)
+
                 else:
                     for spin_channel in sum_k.gf_struct_solver[icrsh].keys():
                         start_sigma[icrsh][spin_channel] << dc_pot[spin_channel]
diff --git a/python/solid_dmft/dmft_tools/interaction_hamiltonian.py b/python/solid_dmft/dmft_tools/interaction_hamiltonian.py
index 590bbc3c..0a7fe487 100644
--- a/python/solid_dmft/dmft_tools/interaction_hamiltonian.py
+++ b/python/solid_dmft/dmft_tools/interaction_hamiltonian.py
@@ -327,8 +327,8 @@ def _construct_dynamic(sum_k, general_params, icrsh):
     mpi.report('###### Dynamic U calculation ######, load parameters from input archive.')
     U_onsite = None
     if mpi.is_master_node():
-        with HDFArchive(general_params['jobname']+'/'+general_params['seedname']+'.h5', 'r') as archive:
-            U_onsite = archive['dynamic_U']['U_scr']
+        with HDFArchive(general_params['jobname']+'/'+general_params['seedname']+'.h5', 'r') as ar:
+            U_onsite = ar['dynamic_U']['U_scr']
     U_onsite = mpi.bcast(U_onsite)
 
     n_orb = common.get_n_orbitals(sum_k)[icrsh]['up']
diff --git a/python/solid_dmft/dmft_tools/manipulate_chemical_potential.py b/python/solid_dmft/dmft_tools/manipulate_chemical_potential.py
index 236aa552..1a897ffa 100755
--- a/python/solid_dmft/dmft_tools/manipulate_chemical_potential.py
+++ b/python/solid_dmft/dmft_tools/manipulate_chemical_potential.py
@@ -31,6 +31,7 @@
 import numpy as np
 
 import triqs.utility.mpi as mpi
+from h5 import HDFArchive
 from triqs.gf import BlockGf, GfImFreq, GfImTime, Fourier, MeshImFreq
 try:
     if mpi.is_master_node():
@@ -242,7 +243,8 @@ def _set_mu_to_gap_middle_with_maxent(general_params, sum_k, gf_lattice_iw, arch
     # Writes spectral function to archive
     if archive is not None:
         unpacked_results = maxent_gf_latt._unpack_maxent_results(maxent_results, mesh)
-        archive['DMFT_results/last_iter']['Alatt_w'] = unpacked_results
+        with HDFArchive(archive, 'a') as ar:
+            ar['DMFT_results/last_iter']['Alatt_w'] = unpacked_results
 
     # Checks if spectral function at Fermi energy below threshold
     spectral_func_threshold = general_params['beta']/np.pi * general_params['mu_gap_gb2_threshold']
@@ -304,7 +306,8 @@ def set_initial_mu(general_params, sum_k, iteration_offset, archive, broadening)
     # If continuing calculation and not updating mu, loads sold value
     if iteration_offset % general_params['mu_update_freq'] != 0:
         if mpi.is_master_node():
-            sum_k.chemical_potential = archive['DMFT_results/last_iter/chemical_potential_pre']
+            with HDFArchive(archive, 'r') as ar:
+                sum_k.chemical_potential = ar['DMFT_results/last_iter/chemical_potential_pre']
         sum_k.chemical_potential = mpi.bcast(sum_k.chemical_potential)
         mpi.report('Chemical potential not updated this step, '
                    + 'reusing loaded one of {:.3f} eV'.format(sum_k.chemical_potential))
@@ -314,7 +317,8 @@ def set_initial_mu(general_params, sum_k, iteration_offset, archive, broadening)
     # chemical_potential_pre from the last run
     previous_mu = None
     if mpi.is_master_node():
-        previous_mu = archive['DMFT_results/last_iter/chemical_potential_pre']
+        with HDFArchive(archive, 'r') as ar:
+            previous_mu = ar['DMFT_results/last_iter/chemical_potential_pre']
     previous_mu = mpi.bcast(previous_mu)
 
     # Runs maxent if spectral weight too low and occupation is close to desired one
diff --git a/python/solid_dmft/dmft_tools/observables.py b/python/solid_dmft/dmft_tools/observables.py
index dc0d25b9..8364ff4c 100644
--- a/python/solid_dmft/dmft_tools/observables.py
+++ b/python/solid_dmft/dmft_tools/observables.py
@@ -31,6 +31,7 @@
 
 # triqs
 import triqs.utility.mpi as mpi
+from h5 import HDFArchive
 from triqs.gf import Gf, MeshImTime
 from triqs.atom_diag import trace_rho_op
 from triqs.gf.descriptors import Fourier
@@ -52,13 +53,14 @@ def prep_observables(h5_archive, sum_k):
         observable array for calculation
     """
 
-    # determine number of impurities
-    n_inequiv_shells = h5_archive['dft_input']['n_inequiv_shells']
+    with HDFArchive(h5_archive, 'r') as ar:
+        # determine number of impurities
+        n_inequiv_shells = ar['dft_input']['n_inequiv_shells']
 
-    # check for previous iterations
-    obs_prev = []
-    if 'observables' in h5_archive['DMFT_results']:
-        obs_prev = h5_archive['DMFT_results']['observables']
+        # check for previous iterations
+        obs_prev = []
+        if 'observables' in ar['DMFT_results']:
+            obs_prev = ar['DMFT_results']['observables']
 
     # prepare observable dicts
     if len(obs_prev) > 0:
diff --git a/python/solid_dmft/dmft_tools/results_to_archive.py b/python/solid_dmft/dmft_tools/results_to_archive.py
index f53f35ae..628f590d 100644
--- a/python/solid_dmft/dmft_tools/results_to_archive.py
+++ b/python/solid_dmft/dmft_tools/results_to_archive.py
@@ -23,6 +23,7 @@
 ################################################################################
 
 import os
+from h5 import HDFArchive
 import triqs.utility.mpi as mpi
 
 
@@ -142,34 +143,35 @@ def write(archive, sum_k, general_params, solver_params, solvers, map_imp_solver
                                        previous_mu, density_mat_pre, density_mat, deltaN, dens)
 
     # Saves the results to last_iter
-    archive['DMFT_results']['iteration_count'] = it
-    for key, value in write_to_h5.items():
-        archive['DMFT_results/last_iter'][key] = value
-
-    # Permanently saves to h5 archive every h5_save_freq iterations
-    if ((not is_sampling and it % general_params['h5_save_freq'] == 0)
-            or (is_sampling and it % general_params['sampling_h5_save_freq'] == 0)):
-
-        archive['DMFT_results'].create_group('it_{}'.format(it))
+    with HDFArchive(archive, 'a') as ar:
+        ar['DMFT_results']['iteration_count'] = it
         for key, value in write_to_h5.items():
-            # Full density matrix only written to last_iter - it is large
-            if 'full_dens_mat_' not in key and 'h_loc_diag_' not in key:
-                archive['DMFT_results/it_{}'.format(it)][key] = value
-
-        # Saves CSC input
-        if general_params['csc']:
-            for dft_var in ['dft_update', 'dft_input', 'dft_misc_input']:
-                if dft_var in archive:
-                    archive['DMFT_results/it_{}'.format(it)].create_group(dft_var)
-                    for key, value in archive[dft_var].items():
-                        # do only store changing elements
-                        if key not in ['symm_kpath', 'kpts_cart']:
-                            archive['DMFT_results/it_{}'.format(it)][dft_var][key] = value
-            for band_elem in ['_bands.dat', '_bands.dat.gnu', '_bands.projwfc_up', '_band.dat']:
-                if os.path.isfile('./{}{}'.format(general_params['seedname'], band_elem)):
-                    os.rename('./{}{}'.format(general_params['seedname'], band_elem),
-                              './{}{}_it{}'.format(general_params['seedname'], band_elem, it))
-            for w90_elem in ['_hr.dat', '.wout']:
-                if os.path.isfile('./{}{}'.format(general_params['seedname'], w90_elem)):
-                    os.rename('./{}{}'.format(general_params['seedname'], w90_elem),
-                              './{}_it{}{}'.format(general_params['seedname'], it, w90_elem))
+            ar['DMFT_results/last_iter'][key] = value
+
+        # Permanently saves to h5 archive every h5_save_freq iterations
+        if ((not is_sampling and it % general_params['h5_save_freq'] == 0)
+                or (is_sampling and it % general_params['sampling_h5_save_freq'] == 0)):
+
+            ar['DMFT_results'].create_group('it_{}'.format(it))
+            for key, value in write_to_h5.items():
+                # Full density matrix only written to last_iter - it is large
+                if 'full_dens_mat_' not in key and 'h_loc_diag_' not in key:
+                    ar['DMFT_results/it_{}'.format(it)][key] = value
+
+            # Saves CSC input
+            if general_params['csc']:
+                for dft_var in ['dft_update', 'dft_input', 'dft_misc_input']:
+                    if dft_var in ar:
+                        ar['DMFT_results/it_{}'.format(it)].create_group(dft_var)
+                        for key, value in ar[dft_var].items():
+                            # do only store changing elements
+                            if key not in ['symm_kpath', 'kpts_cart']:
+                                ar['DMFT_results/it_{}'.format(it)][dft_var][key] = value
+                for band_elem in ['_bands.dat', '_bands.dat.gnu', '_bands.projwfc_up', '_band.dat']:
+                    if os.path.isfile('./{}{}'.format(general_params['seedname'], band_elem)):
+                        os.rename('./{}{}'.format(general_params['seedname'], band_elem),
+                                  './{}{}_it{}'.format(general_params['seedname'], band_elem, it))
+                for w90_elem in ['_hr.dat', '.wout']:
+                    if os.path.isfile('./{}{}'.format(general_params['seedname'], w90_elem)):
+                        os.rename('./{}{}'.format(general_params['seedname'], w90_elem),
+                                  './{}_it{}{}'.format(general_params['seedname'], it, w90_elem))
diff --git a/python/solid_dmft/dmft_tools/solver.py b/python/solid_dmft/dmft_tools/solver.py
index c5e22f7d..21b90254 100755
--- a/python/solid_dmft/dmft_tools/solver.py
+++ b/python/solid_dmft/dmft_tools/solver.py
@@ -62,11 +62,7 @@
     from solid_dmft.dmft_tools.solvers.ftps_interface import FTPSInterface
     interfaces_dict['ftps'] = FTPSInterface
 
-
-
-
 # generic dispatch function for the solver classes
-
 def create_solver(general_params, solver_params, sum_k, icrsh, h_int, iteration_offset,
                     deg_orbs_ftps, gw_params=None, advanced_params=None):
         '''
diff --git a/test/python/lno_hubbardI_mag/ref.h5 b/test/python/lno_hubbardI_mag/ref.h5
index 6a6cab04..62e24045 100644
Binary files a/test/python/lno_hubbardI_mag/ref.h5 and b/test/python/lno_hubbardI_mag/ref.h5 differ
diff --git a/test/python/test_afm_mapping.py b/test/python/test_afm_mapping.py
index 499888f6..3cd8ec7e 100755
--- a/test/python/test_afm_mapping.py
+++ b/test/python/test_afm_mapping.py
@@ -7,6 +7,7 @@
 import sys
 sys.path.append('./src')
 
+from h5 import HDFArchive
 from solid_dmft.dmft_tools import afm_mapping
 from helper import are_iterables_equal
 
@@ -15,7 +16,10 @@
 class test_afm_mapping(unittest.TestCase):
     def test_determine_afm_mapping(self):
         general_params = {'magmom': [+1, -1, +1, -1], 'afm_order': True}
-        archive = {'DMFT_input': {}}
+        archive = 'test.h5'
+        with HDFArchive(archive, 'w') as ar:
+            ar.create_group('DMFT_input')
+            ar['DMFT_input']= {}
         n_inequiv_shells = 4
 
         expected_general_params = general_params.copy()
@@ -28,7 +32,9 @@ def test_determine_afm_mapping(self):
         assert are_iterables_equal(general_params, expected_general_params)
 
         general_params = {'magmom': [+1, -1, +2, +2], 'afm_order': True}
-        archive = {'DMFT_input': {}}
+        with HDFArchive(archive, 'w') as ar:
+            ar.create_group('DMFT_input')
+            ar['DMFT_input']= {}
         n_inequiv_shells = 4
 
         expected_general_params = general_params.copy()
@@ -42,8 +48,10 @@ def test_determine_afm_mapping(self):
 
         # Reading in the afm_mapping from the archive
         general_params = {'magmom': [+1, -1, +2], 'afm_order': True}
-        archive = {'DMFT_input': {'afm_mapping': [[False, 0, False], [False, 1, False],
-                                                  [False, 2, False]]}}
+        with HDFArchive(archive, 'w') as ar:
+            ar.create_group('DMFT_input')
+            ar['DMFT_input']= {'afm_mapping': [[False, 0, False], [False, 1, False],
+                                                  [False, 2, False]]}
         n_inequiv_shells = 3
 
         expected_general_params = general_params.copy()
@@ -56,7 +64,9 @@ def test_determine_afm_mapping(self):
         assert are_iterables_equal(general_params, expected_general_params)
 
         general_params = {'magmom': [+1, -1, +2, +2], 'afm_order': True}
-        archive = {'DMFT_input': {}}
+        with HDFArchive(archive, 'w') as ar:
+            ar.create_group('DMFT_input')
+            ar['DMFT_input']= {}
         n_inequiv_shells = 3
 
         expected_general_params = general_params.copy()