diff --git a/CMakeLists.txt b/CMakeLists.txt
index b6400873f4..a27699b7dc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1023,6 +1023,18 @@ if(NOT QMC_NO_SLOW_CUSTOM_TESTING_COMMANDS)
   endif()
 endif()
 
+#-------------------------------------------------------------------
+# Check if GPAW is available
+#-------------------------------------------------------------------
+if(NOT QMC_NO_SLOW_CUSTOM_TESTING_COMMANDS)
+  check_python_reqs(gpaw "" HAVE_GPAW)
+  if(NOT ${HAVE_GPAW})
+    message(STATUS "Unable to import GPAW python module. GPAW converter tests will not be run.")
+  else()
+    message(STATUS "Successfully imported GPAW python modules.")
+  endif()
+endif()
+
 #-------------------------------------------------------------------
 # Check if QE executables available at QE_BIN or on the PATH
 #-------------------------------------------------------------------
diff --git a/src/QMCTools/CMakeLists.txt b/src/QMCTools/CMakeLists.txt
index f4915d8ef2..1d791b5ffb 100644
--- a/src/QMCTools/CMakeLists.txt
+++ b/src/QMCTools/CMakeLists.txt
@@ -53,6 +53,9 @@ if(USE_OBJECT_TARGET)
   target_link_libraries(qmcfinitesize qmcparticle qmcutil)
 endif()
 
+file(CREATE_LINK ${CMAKE_CURRENT_SOURCE_DIR}/gpaw4qmcpack.py ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/gpaw4qmcpack.py SYMBOLIC)
+install(FILES gpaw4qmcpack.py DESTINATION bin)
+
 foreach(EXE_TARGET convert4qmc qmc-extract-eshdf-kvectors qmc-get-supercell qmc-check-affinity convertpw4qmc qmcfinitesize)
   add_test_target_in_output_location(${EXE_TARGET} bin)
   install(TARGETS ${EXE_TARGET} RUNTIME DESTINATION bin)
diff --git a/src/QMCTools/Eshdf.py b/src/QMCTools/Eshdf.py
new file mode 100644
index 0000000000..ffffaa6195
--- /dev/null
+++ b/src/QMCTools/Eshdf.py
@@ -0,0 +1,245 @@
+######################################################################################
+## This file is distributed under the University of Illinois/NCSA Open Source License.
+## See LICENSE file in top directory for details.
+##
+## Copyright (c) 2021 QMCPACK developers.
+##
+## File developed by: Juha Tiihonen, tiihonen@iki.fi, University of Jyvaskyla
+##
+## File created by: Juha Tiihonen, tiihonen@iki.fi, University of Jyvaskyla 
+#######################################################################################
+
+from numpy import array
+
+# This file defines python classes to load/write data needed by the ESHDF file format.
+#
+# The main logic of the class structure is to inherit from Base -> Representation -> Application
+#
+# currently implemented:
+#
+# Base:
+#   EshdfFileBase
+# Representation:
+#   EshdfFilePw
+# Application:
+#   EshdfFilePwGpaw  <- found in Eshdf_gpaw.py
+
+
+# define dtypes for the HDF5
+int_dtype     = '<i4'   # default integer dtype
+str_dtype     = 'S'     # default string dtype
+float_dtype   = 'f'     # default float dtype
+eshdf_version = [2,1,0] # current version
+
+
+# Class for atomic species
+class AtomicSpecies():
+    id      = None
+    mass    = None
+    Z       = None
+    name    = None
+    valence = None
+    def __init__(
+            self,
+            id,
+            name,
+            mass,
+            Z,
+            valence
+            ):
+        id           = id
+        self.mass    = mass
+        self.Z       = Z
+        self.name    = namei
+        self.valence = valence
+    #end def
+#end class
+
+
+# Base class for an ESHDF file
+class EshdfFileBase():
+    
+    # class flags and defaults
+    loaded      = False
+    outfile     = 'eshdf.h5'
+
+    # init null or default ESHDF data
+
+    # /application must be overridden
+    code    = ''      # string
+    version = [0,0,1] # 3-int list
+    # /atoms must be overridden
+    num_atoms   = 0
+    num_species = 0
+    positions   = [[]] # num_atoms x 3 array
+    species     = []   # list of AtomicSpecies (num_species)
+    species_ids = []   # list of integers (num_atoms)
+    # /electrons must be overridden
+    num_elec = 0
+    num_spin = 0
+    # /format
+    fmt     = 'ES-HDF'
+    # /supercell must be overridden
+    supercell = [[1,0,0],[0,1,0],[0,0,1]] # 3x3 matrix
+    # /version
+    version = eshdf_version
+
+    # standard init: infile, outfile, **kwargs
+    def __init__(
+        self,
+        infile    = None,
+        outfile   = None,
+        **kwargs
+        ):
+        self.infile  = infile
+        self.outfile = outfile
+        if infile is not None:
+            self.load(infile=infile,**kwargs)
+        #end if
+    #end def
+
+    # standard load: infile, **kwargs
+    #   override this in implementations to read data
+    def load(
+        self,
+        infile  = None,
+        **kwargs,
+        ):
+        print('Loading of data not implemented.')
+        self.loaded = True
+    #end def
+
+    # check integrity of the data
+    def check_integrity(self):
+        #TODO: various tests
+        return True
+    #end def
+
+    # standard write: outfile
+    #   default: write to HDF5
+    def write(self, outfile = None):
+        if not self.loaded:
+            raise(DataNotLoaded('Must load data first'))
+        #end if
+        if outfile is None:
+            outfile = self.outfile
+        #end if
+        if self.check_integrity():
+            self._write_h5(outfile)
+        #end if
+    #end def
+
+    # write to HDF5
+    def _write_h5(self, outfile):
+        import h5py as h5
+        with h5.File(outfile,mode='w') as f:
+            self._write_h5_application(f)
+            self._write_h5_atoms(f)
+            self._write_h5_electrons(f)
+            self._write_h5_format(f)
+            self._write_h5_supercell(f)
+            self._write_h5_version(f)
+        #end with
+    #end def
+
+    def _write_h5_application(self,f):
+        g = f.create_group('application')
+        g.create_dataset('code',shape=1,data=array(self.application,dtype=str_dtype))
+        g.create_dataset('version',shape=3,data=array(self.app_version,dtype=int_dtype))
+    #end def
+
+    def _write_h5_atoms(self,f):
+        g  = f.create_group('atoms')
+        g.create_dataset('number_of_atoms',data=array([self.num_atoms],dtype=int_dtype))
+        g.create_dataset('number_of_species',data=array([self.num_species],dtype=int_dtype))
+        g.create_dataset('positions',data=self.positions)
+        for s,sid in enumerate(self.species):
+            gs = g.create_group('species_'+str(s))
+            gs.create_dataset('atomic_number',data=array([self.species_Z[s]]),dtype=int_dtype)
+            gs.create_dataset('mass',data=array([self.species_mass[s]],dtype=float_dtype))
+            gs.create_dataset('name',data=array([self.species_name[s]],dtype=str_dtype))
+            gs.create_dataset('valence_charge',data=array([self.species_valence[s]],dtype=float_dtype))
+        #end for
+        g.create_dataset('species_ids',data=array(self.species_ids),dtype=int_dtype)
+    #end def
+
+    # writing of the orbital information must be overridden based on the representation
+    #   Pw, real-space, etc
+    def _write_h5_electrons(self,f):
+        g = f.create_group('electrons')
+        g.create_dataset('number_of_electrons',data=self.num_elec)
+        g.create_dataset('number_of_spins',data=array([self.num_spin],dtype=int_dtype))
+        print('To write kpoint/orbital information, must not use the base class EshdfFileBase\nbut an appropriate representation, such as EshdfFilePw!')
+    #end def
+
+    def _write_h5_format(self,f):
+        f.create_dataset('format',data=array(self.fmt,dtype=str_dtype))
+    #end def
+
+    def _write_h5_supercell(self,f):
+        g = f.create_group('supercell')
+        g.create_dataset('primitive_vectors',data=self.supercell)
+    #end def
+    
+    def _write_h5_version(self,f):
+        f.create_dataset('version',data=array(self.version,dtype=int_dtype))
+    #end def
+
+#end class
+
+
+# Class for an ESHDF file in the PW representation
+class EshdfFilePw(EshdfFileBase):
+
+    # /electrons
+    #  WIP
+    density   = False
+    rho_gv    = []
+    rho_g     = []
+    kpts      = []
+    gvectors  = []
+    numsym    = []
+    eps       = []
+    num_bands = [[]]
+    psi_g     = [[[]]]
+    symgroup  = []
+    W_kpt     = []
+
+    def _write_h5_electrons(self,f):
+        g = f.create_group('electrons')
+        if self.density:
+            gd = g.create_group('density')
+            gd.create_dataset('gvectors',data=self.rho_gv)
+            gd.create_dataset('number_of_gvectors',data=array([self.rho_gv.shape[0]],dtype=int_dtype))
+            for s in range(self.num_spin):
+                gd0 = gd.create_group('spin_{}'.format(s))
+                gd0.create_dataset('density_g',data=self.rho_g[s])
+            #end for
+        #end if
+        for k,kpt in enumerate(self.kpts):
+            kp0 = g.create_group('kpoint_{}'.format(k))
+            if k==0:
+                kp0.create_dataset('gvectors',data=self.gvectors)
+                kp0.create_dataset('number_of_gvectors',data=array([self.gvectors.shape[0]],dtype=int_dtype))
+            #end if
+            kp0.create_dataset('numsym',data=array([self.numsym[k]],dtype=int_dtype))
+            kp0.create_dataset('reduced_k',data=kpt)
+            for s in range(self.num_spin):
+                kp00 = kp0.create_group('spin_{}'.format(s))
+                kp00.create_dataset('eigenvalues',data=self.eps[k][s])
+                kp00.create_dataset('number_of_states',data=array([self.num_bands[k][s]],dtype=int_dtype))
+                for n in range(self.num_bands[k][s]):
+                    kp00s = kp00.create_group('state_{}'.format(n))
+                    kp00s.create_dataset('psi_g',data=self.psi_g[k][s][n])
+                #end for
+            #end for
+            kp0.create_dataset('symgroup',data=array([self.symgroup[k]],dtype=int_dtype))
+            kp0.create_dataset('weight',data=array([self.W_kpt[k]],dtype=float_dtype))
+        #end for
+        g.create_dataset('number_of_electrons',data=self.num_elec)
+        g.create_dataset('number_of_kpoints',data=array([len(self.kpts)],dtype=int_dtype))
+        g.create_dataset('number_of_spins',data=array([self.num_spin],dtype=int_dtype))
+    #end def
+
+#end class
+
diff --git a/src/QMCTools/Eshdf_gpaw.py b/src/QMCTools/Eshdf_gpaw.py
new file mode 100644
index 0000000000..2223577d49
--- /dev/null
+++ b/src/QMCTools/Eshdf_gpaw.py
@@ -0,0 +1,283 @@
+######################################################################################
+## This file is distributed under the University of Illinois/NCSA Open Source License.
+## See LICENSE file in top directory for details.
+##
+## Copyright (c) 2021 QMCPACK developers.
+##
+## File developed by: Juha Tiihonen, tiihonen@iki.fi, University of Jyvaskyla
+##
+## File created by: Juha Tiihonen, tiihonen@iki.fi, University of Jyvaskyla 
+#######################################################################################
+
+# This file implements Eshdf bindings (Eshdf.py) for reading GPAW orbitals
+
+from numpy import array,pi,dot,real,imag,concatenate,unique,zeros
+from Eshdf import EshdfFilePw
+
+# Calculate vector norm of 2-float representation of a complex vector
+def calculate_norm(v):
+    return ((v**2).sum())**0.5
+#end def
+
+# Helper function (obsolete) to match gvector g in an array of gvectors gset
+# and return the value of the associated density coefficients from dset
+def match_gv(g,gset,dset):
+    g0 = gset[:,0]==g[0]
+    g1 = gset[:,1]==g[1]
+    g2 = gset[:,2]==g[2]
+    return dset[logical_and(logical_and(g0,g1),g2),:][0]
+#end def
+
+# Add the negative half of gvectors removed by symmetry (full array expected by QMCPACK)
+def desymmetrize_z(gv,psi):
+    gv2  = -gv[1:].copy()
+    psi2 = psi[1:].copy()
+    psi2[:,1] *= -1
+    gv2  = concatenate((gv, gv2), axis=0)
+    psi2 = concatenate((psi,psi2),axis=0)
+    return gv2,psi2
+#end def
+
+# Eshdf File class from GPAW
+class EshdfFilePwGpaw(EshdfFilePw):
+
+    # if infile and outfile provided, load automatically
+    def __init__(
+        self,
+        infile    = None,
+        outfile   = None,
+        **kwargs
+        ):
+        self.infile  = infile
+        self.outfile = outfile
+        if infile is not None:
+            self.load(infile=infile,**kwargs)
+        #end if
+    #end def
+
+    # initiate loading of data from the provided objects
+    #   Compute and add charge density, if density=True
+    def load(
+        self,
+        infile  = None,
+        calc    = None,
+        atoms   = None,
+        density = False,
+        **kwargs,
+        ):
+        self.density = density
+        if infile is not None:
+            if '.gpw' in infile:
+                self.load_gpaw_restart(infile,**kwargs)
+            else:
+                print('Unsupported file {} not loaded.'.format(infile))
+            #end if
+        elif calc is not None and atoms is not None:
+            self.load_gpaw_calc(calc,atoms)
+        else:
+            print('Could not load data for ESHDF.')
+        #end if
+    #end def
+
+    # load data from GPAW restart file
+    def load_gpaw_restart(
+        self,
+        infile,
+        **kwargs,
+        ):
+        from gpaw import restart
+        atoms,calc = restart(infile)
+        try: # check if the wavefunction exists
+            wfs = calc.wfs
+        except:
+            print('ERROR: No wavefunctions found in the GPAW restart file {}'.format(filename))
+            print("Make sure to write the GPAW file with mode='all'!")
+            raise NoWavefunctionFound
+        #end try
+        if self.density: 
+            # run SCF to re-populate density
+            print('Running SCF to re-populate density')
+            calc.scf.converged = False
+            calc.results.pop('energy')
+            E = calc.get_potential_energy()
+        #end def
+        self.load_gpaw_calc(calc,atoms)
+    #end def
+
+    # GPAW: load data from calc and atoms objects
+    #   Currently only supports PW mode
+    def load_gpaw_calc(
+        self,
+        calc,
+        atoms,
+        ):
+        mode = calc.wfs.mode
+        if mode=='pw':
+            self._load_application()
+            self._load_atoms(calc,atoms)
+            self._load_electrons(calc,atoms)
+        else:
+            print('Could not load data from GPAW Calculator in mode={}. Only pw mode is currently supported.'.format(mode))
+        #end if
+        self.loaded = True
+    #end def
+
+    # GPAW: load "application" data
+    def _load_application(self):
+        from gpaw import __version__ as gpaw_version
+        from re import sub
+        self.application = 'gpaw'
+        self.app_version = sub('\D',' ',gpaw_version).split()[:3]
+    #end def
+
+    # GPAW: load "atoms" data
+    def _load_atoms(self,calc,atoms):
+        from ase.units import Bohr
+        num_atoms = atoms.get_global_number_of_atoms()
+        supercell = calc.wfs.gd.cell_cv
+        syms      = atoms.get_chemical_symbols()
+        masses    = atoms.get_masses()
+        positions = atoms.get_positions()/Bohr # convert from A to Bohr
+        Zs        = atoms.get_atomic_numbers()
+        # figure out species data
+        sym_set         = set(syms)
+        num_species     = len(sym_set)
+        species_ids     = []
+        species_mass    = []
+        species_Z       = []
+        species_name    = []
+        species_valence = []
+        for s,name in enumerate(sym_set):
+            i = syms.index(name)
+            species_mass.append(masses[i])
+            species_Z.append(Zs[i])
+            species_valence.append(Zs[i]) # fixme for PP calculations
+            species_name.append(name)
+            for sym in syms:
+                if sym==name:
+                    species_ids.append(s)
+                #end if
+            #end for
+        #end for
+        self.species         = sym_set
+        self.num_species     = num_species
+        self.num_atoms       = num_atoms
+        self.positions       = positions
+        self.species_ids     = species_ids
+        self.species_mass    = species_mass
+        self.species_valence = species_valence
+        self.species_Z       = species_Z
+        self.species_name    = species_name
+        self.supercell       = supercell
+    #end def
+
+    # GPAW: transform reciprocal vectors to Miller indices
+    def _get_pw_indices(self,pd,k):
+        return array(dot(pd.get_reciprocal_vectors(k,add_q=False),pd.gd.cell_cv/2/pi).round(0),dtype=int)
+    #end def
+
+    # GPAW: load "electrons" data
+    def _load_electrons(self,calc,atoms):
+        mesh     = calc.wfs.gd.N_c
+        num_spin = calc.wfs.nspins
+        N_tot    = calc.wfs.nvalence
+        if num_spin==1:
+            N_up = N_tot // 2
+            N_dn = N_tot - N_up
+        else: # spin-polarized # FIXME
+            N_up = N_tot // 2
+            N_dn = N_tot - N_up
+        #end if
+        # kpoints
+        num_kpt  = calc.wfs.kd.get_count()
+        W_kpt    = calc.wfs.kd.weight_k
+        try:
+            kpts    = calc.wfs.kd.get_bz_q_points()
+        except:
+            kpts    = calc.wfs.kd.bzk_kc
+        #end try
+
+        # Combine the GPAW plane wave data for each separate kpoint into one big array, gvall
+        gvall = self._get_pw_indices(calc.wfs.pd,0)
+        imin  = 0
+        imax  = len(gvall)
+        iis   = [[imin,imax]] # indices to mark the beginning and end of each kpoint
+        for k in range(1,num_kpt):
+            gv    = self._get_pw_indices(calc.wfs.pd,k)
+            gvall = concatenate((gvall,gv),axis=0)
+            imin  = imax
+            imax += len(gv)
+            iis.append([imin,imax])
+        #end for
+        # Gv is a new gvector basis with only unique vectors
+        # indices help to associate the pw correct coefficients later on
+        Gv,indices = unique(gvall,axis=0,return_inverse=True)
+
+        data_numsym     = []
+        data_symgroup   = []
+        data_reduced_k  = []
+        data_occupation = []
+        data_num_bands  = []
+        data_eps        = []
+        data_psi_g      = []
+        for k,kpt in enumerate(kpts):
+            numsym = 1   # FIXME
+            symgroup = 1 # FIXME
+            data_numsym.append(numsym)
+            data_symgroup.append(symgroup)
+            ii = iis[k]
+            # TODO: checks for spin polarized calculations?
+            datas_eps       = []
+            datas_psi_g     = []
+            datas_num_bands = []
+            for s in range(num_spin):
+                occ       = calc.wfs.kpt_u[k].f_n
+                eps       = calc.wfs.kpt_u[k].eps_n
+                num_bands = len(eps)
+                u         = k*num_spin+s
+                datasn_psi_g = []
+                for n in range(num_bands):
+                    # get wavefunction in PW
+                    wf_c  = calc.wfs.kpt_u[u].psit_nG[n]
+                    wf    = array([real(wf_c),imag(wf_c)])
+                    psi   = zeros((len(Gv),2))
+                    psi[indices[ii[0]:ii[1]]] = wf.T/calculate_norm(wf) # normalize the crude way
+                    datasn_psi_g.append(psi)
+                #end for
+                datas_num_bands.append(num_bands)
+                datas_eps.append(eps)
+                datas_psi_g.append(datasn_psi_g)
+            #end for
+            data_eps.append(datas_eps)
+            data_num_bands.append(datas_num_bands)
+            data_psi_g.append(datas_psi_g)
+        #end for
+        self.gvectors   = Gv
+        self.num_elec   = [N_up,N_dn]
+        self.kpts       = kpts
+        self.numsym     = data_numsym
+        self.symgroup   = data_symgroup
+        self.reduced_k  = data_reduced_k
+        self.W_kpt      = W_kpt
+        self.num_bands  = data_num_bands
+        self.num_spin   = num_spin
+        self.eps        = data_eps
+        self.psi_g      = data_psi_g
+        if self.density:
+            # FIXME: spin polarized
+            # read data
+            rho_gv   = self._get_pw_indices(calc.density.pd2,0)
+            rho_psic = calc.density.nt_Q
+            # process
+            rho_psi  = array([real(rho_psic),imag(rho_psic)]).T
+            gv,rho   = desymmetrize_z(rho_gv,rho_psi)
+            n        = N_tot/rho[0,0]/calc.density.gd.volume # normalize at the gamma point
+            # save
+            self.rho_gv   = gv
+            self.rho_g    = rho*n
+        #end if
+    #end def
+
+    # file writing functions are defined in the base classes
+
+#end class
diff --git a/src/QMCTools/README.md b/src/QMCTools/README.md
new file mode 100644
index 0000000000..9b7c2150a5
--- /dev/null
+++ b/src/QMCTools/README.md
@@ -0,0 +1,44 @@
+# ESHDF Python class
+
+ESHDF (Electronic Structure Hierarchical Data Format) is the file
+specification used to provide orbital information for QMCPACK. It is based on
+the Hierarchical Data Format version 5 (HDF5; 
+https://www.hdfgroup.org/solutions/hdf5/).
+
+A Python class for ESHDF is defined in ´eshdf.py´, using ´h5py´ Python
+library. The class hierarchy is intended for reading, managing and writing of
+ESHDF data.
+
+The classes are inherited as follows: `Base -> Representation -> Application.`
+The Base class ´EshdfFileBase´ contains virtual implementation of the core
+methods. The Representation classes add methods and variables for managing
+specific choices of the representation. The Application class implements
+bindings to a specific application.
+
+
+## Currently implemented classes
+
+Base:
+- EshdfFileBase
+Representation:
+- EshdfFilePw (plane-wave)
+Application:
+- EshdfFilePwGpaw (´Eshdf_gpaw.py´)
+
+This is work-in-progress, and other features are coming up.
+
+
+# GPAW4QMCPACK
+
+Converter of electronic orbitals from GPAW (https://wiki.fysik.dtu.dk/gpaw/)
+is implemented in ´gpaw4qmcpack.py´, using the ESHDF Python class.
+
+The converter is run on the command line. It requires a GPAW restart file
+(.gpw) with orbitals included with ´mode=all´, and produces an ESHDF file (.h5)
+that can be used in QMCPACK.
+
+If the MPC correction will be needed, the charge density information can be 
+added with `-d`, but it must be computed on the runtime. Therefore, the 
+density is not included by default.
+
+This is work-in-progress, and other features and documentation are coming up.
diff --git a/src/QMCTools/gpaw4qmcpack.py b/src/QMCTools/gpaw4qmcpack.py
new file mode 100755
index 0000000000..936d03d6a7
--- /dev/null
+++ b/src/QMCTools/gpaw4qmcpack.py
@@ -0,0 +1,29 @@
+#!/usr/bin/env python3
+######################################################################################
+## This file is distributed under the University of Illinois/NCSA Open Source License.
+## See LICENSE file in top directory for details.
+##
+## Copyright (c) 2021 QMCPACK developers.
+##
+## File developed by: Juha Tiihonen, tiihonen@iki.fi, University of Jyvaskyla
+##
+## File created by: Juha Tiihonen, tiihonen@iki.fi, University of Jyvaskyla 
+#######################################################################################
+
+# This file converts GPAW orbitals for QMCPACK using ESHDF Python classes
+
+import argparse
+
+from Eshdf_gpaw import EshdfFilePwGpaw
+
+parser = argparse.ArgumentParser(description='GPAW4QMCPACK')
+parser.add_argument('infile',help='input .gpw restart file')
+parser.add_argument('outfile',default='eshdf.h5',help='output .h5 orbital file')
+parser.add_argument('-d','--density',action='store_true',help='whether or not to convert density')
+
+if __name__=='__main__':
+    args = vars(parser.parse_args())
+    eshdf = EshdfFilePwGpaw(**args)
+    eshdf.write()
+    print('Converted GPAW orbitals for QMCPACK!')
+#end if
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 20f58437bf..8cfbfbea88 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -66,6 +66,7 @@ if(HAVE_PYSCF)
     endif()
   endif()
 endif()
+
 #
 #
 if(NOT ${QMCAPP_ERR} STREQUAL "")
diff --git a/tests/converter/CMakeLists.txt b/tests/converter/CMakeLists.txt
index 1d563ab41a..357032222f 100644
--- a/tests/converter/CMakeLists.txt
+++ b/tests/converter/CMakeLists.txt
@@ -63,3 +63,12 @@ add_converter_test(test_LiH_qp)
 add_converter_test(test_C4_MSD_excited)
 add_converter_test(test_diamond2_rmg)
 add_converter_test(test_O_rmg)
+
+if(HAVE_GPAW)
+    include("${qmcpack_SOURCE_DIR}/CMake/python.cmake")
+    check_python_reqs("numpy;h5py;ase;gpaw" gpaw4qmcpack_test add_tests)
+    if(add_tests)
+      message("Python dependencies met. Adding GPAW4QMCPACK tests")
+      add_subdirectory("gpaw")
+    endif()
+endif()
diff --git a/tests/converter/gpaw/CMakeLists.txt b/tests/converter/gpaw/CMakeLists.txt
new file mode 100644
index 0000000000..eb8fcf1ff6
--- /dev/null
+++ b/tests/converter/gpaw/CMakeLists.txt
@@ -0,0 +1,18 @@
+# add tests for GPAW4QMCPACK converter
+file(CREATE_LINK ${CMAKE_CURRENT_SOURCE_DIR}/gpaw4qmcpack_test.py ${CMAKE_CURRENT_BINARY_DIR}/gpaw4qmcpack_test.py SYMBOLIC)
+function(ADD_GPAW4QMCPACK_TEST test_name)
+  set(EXE_NAME "${PROJECT_SOURCE_DIR}/src/QMCTools/gpaw4qmcpack.py")
+  file(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/${test_name}")
+  file(CREATE_LINK "${CMAKE_CURRENT_SOURCE_DIR}/${test_name}/restart.gpw"
+	  "${CMAKE_CURRENT_BINARY_DIR}/${test_name}/restart.gpw" SYMBOLIC)
+  file(CREATE_LINK "${CMAKE_CURRENT_SOURCE_DIR}/${test_name}/gold.orbs.h5"
+	  "${CMAKE_CURRENT_BINARY_DIR}/${test_name}/gold.orbs.h5" SYMBOLIC)
+  add_test(NAME gpaw4qmcpack_test COMMAND "${CMAKE_CURRENT_BINARY_DIR}/gpaw4qmcpack_test.py"
+      "${CMAKE_CURRENT_BINARY_DIR}/${test_name}"
+       --exe "${EXE_NAME}" --h5diff "${HDF5_DIFF_EXECUTABLE}")
+  set_tests_properties(gpaw4qmcpack_test PROPERTIES
+                         FAIL_REGULAR_EXPRESSION "  FAIL"
+                         PASS_REGULAR_EXPRESSION "  pass"
+                         TIMEOUT 120 LABELS "converter;deterministic")
+endfunction()
+add_gpaw4qmcpack_test(test_Si_diamond)
diff --git a/tests/converter/gpaw/gpaw4qmcpack_test.py b/tests/converter/gpaw/gpaw4qmcpack_test.py
new file mode 100755
index 0000000000..8be1f5813c
--- /dev/null
+++ b/tests/converter/gpaw/gpaw4qmcpack_test.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+import subprocess
+import sys
+
+# unit test GPAW4QMCPACK
+
+# TODO: implement and test converter generated QMCPACK input
+
+def run_gpaw4qmcpack_test(test_name, g4q_exe, h5diff_exe):
+    okay = True
+
+    # Example invocation of converter
+    #gpaw4qmcpack infile.gpw outfile.h5 --density
+
+    infile = 'restart.gpw'
+    outfile = 'test.orbs.h5'
+    cmd = [g4q_exe,infile,outfile]
+
+    p = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, universal_newlines=True)
+    stdout, stderr = p.communicate()
+
+    file_out = open('stdout.txt', 'w')
+    file_out.write(stdout)
+    file_out.close()
+    if len(stderr) > 0 :
+        file_err = open('stderr.txt', 'w')
+        file_err.write(stderr)
+        file_err.close()
+
+    ret = p.returncode
+
+    if ret != 0:
+        print("Return code nonzero: ", ret)
+        okay = False
+
+    if len(stderr.strip()) != 0:
+        # some MPI output on stderr is okay
+        # TODO - more general way of checking acceptable stderr strings
+        if not stderr.startswith('Rank'):
+            print("Stderr not empty")
+            print(stderr)
+            okay = False
+
+    ret = os.system(h5diff_exe + ' -d 0.000001 gold.orbs.h5 test.orbs.h5')
+    # if it's okay up to this point
+    if ret==0 and okay:
+        print("  pass")
+        return True
+    else:
+        print("h5diff reported a difference")
+        print("  FAIL")
+        return False
+
+    return okay
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Test gpaw4qmcpck')
+    parser.add_argument('test_name',
+                        default='test_Si_diamond',
+                        help='Name of test to run (name of directory)')
+    parser.add_argument('--exe',
+                        default='gpaw4qmcpack',
+                        help='Location of gpaw4qmcpack executable')
+    parser.add_argument('--h5diff',
+                        default='h5diff',
+                        help='Location of h5diff executable')
+    args = parser.parse_args()
+
+    test_dir = args.test_name
+    if not os.path.exists(test_dir):
+        print("Test not found: ", test_dir)
+        sys.exit(1)
+
+    curr_dir = os.getcwd()
+    os.chdir(test_dir)
+
+    ret = run_gpaw4qmcpack_test(test_dir, args.exe, args.h5diff)
+
+    os.chdir(curr_dir)
+
+    if ret:
+        sys.exit(0)
+    sys.exit(1)
+
diff --git a/tests/converter/gpaw/test_Si_diamond/gold.orbs.h5 b/tests/converter/gpaw/test_Si_diamond/gold.orbs.h5
new file mode 100644
index 0000000000..f58a60b043
Binary files /dev/null and b/tests/converter/gpaw/test_Si_diamond/gold.orbs.h5 differ
diff --git a/tests/converter/gpaw/test_Si_diamond/restart.gpw b/tests/converter/gpaw/test_Si_diamond/restart.gpw
new file mode 100644
index 0000000000..228266c850
Binary files /dev/null and b/tests/converter/gpaw/test_Si_diamond/restart.gpw differ
diff --git a/tests/converter/gpaw/test_Si_diamond/run_si.py b/tests/converter/gpaw/test_Si_diamond/run_si.py
new file mode 100755
index 0000000000..92c2049409
--- /dev/null
+++ b/tests/converter/gpaw/test_Si_diamond/run_si.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+
+# This file provided just for reference, not executed in the test
+
+from ase.build import bulk
+from ase.units import Ry
+from gpaw import GPAW,PW,FermiDirac,setup_paths
+
+silicon = bulk('Si', 'diamond', a=5.459)
+
+econv = 1e-9*2*Ry
+ecut  = 30*Ry
+silicon.calc = GPAW(xc='PBE',
+        setups      = {'default':'sg15'}, # Based on standard library of soft SG15 ecps
+        mode        = PW(ecut,force_complex_dtype=True), # give all modes explicitly
+        txt         = 'silicon.txt',
+        kpts        = (2,2,2),
+        occupations = FermiDirac(0.1),
+        convergence = {'energy':econv},
+        symmetry    = 'off'  
+        )
+silicon.get_potential_energy()
+silicon.calc.write('restart.gpw',mode='all') # important to request mode='all'