diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 03dffc89..15fb6ab4 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -496,6 +496,13 @@ add_test(NAME DavidsonReplicated
          ${CMAKE_CURRENT_SOURCE_DIR}/DavidsonReplicated/davidson.cfg
          ${CMAKE_CURRENT_SOURCE_DIR}/Davidson/coords.in
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+ add_test(NAME DavidsonMixing
+         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/DavidsonMixing/test.py
+         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+         ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt
+         ${CMAKE_CURRENT_SOURCE_DIR}/DavidsonMixing/davidson.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/DavidsonMixing/sih4.xyz
+         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 add_test(NAME testSpinO2
          COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/SpinO2/test.py
          ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
diff --git a/tests/DavidsonMixing/davidson.cfg b/tests/DavidsonMixing/davidson.cfg
new file mode 100644
index 00000000..39c7d457
--- /dev/null
+++ b/tests/DavidsonMixing/davidson.cfg
@@ -0,0 +1,35 @@
+verbosity=2
+xcFunctional=PBE
+FDtype=4th
+[Mesh]
+nx=40
+ny=40
+nz=40
+[Domain]
+ox=-6.75
+oy=-6.75
+oz=-6.75
+lx=13.5
+ly=13.5
+lz=13.5
+[Potentials]
+pseudopotential=pseudo.Si_ONCV_PBE_SG15
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+[Run]
+type=QUENCH
+[Quench]
+solver=Davidson
+max_steps=100
+atol=1.e-8
+[Orbitals]
+nempty=1
+initial_type=random
+temperature=300.
+[ProjectedMatrices]
+solver=exact
+replicated=true
+[DensityMatrix]
+nb_inner_it=1
+mixing=0.8
+[Restart]
+output_level=0
diff --git a/tests/DavidsonMixing/sih4.xyz b/tests/DavidsonMixing/sih4.xyz
new file mode 100644
index 00000000..1294ac62
--- /dev/null
+++ b/tests/DavidsonMixing/sih4.xyz
@@ -0,0 +1,8 @@
+5
+SiH4 molecule (coordinates in Angstrom)
+Si       0.0   0.0   0.0
+H        0.863   0.863   0.863
+H       -0.863  -0.863   0.863
+H       -0.863   0.863  -0.863
+H        0.863  -0.863  -0.863 
+
diff --git a/tests/DavidsonMixing/test.py b/tests/DavidsonMixing/test.py
new file mode 100755
index 00000000..16383c01
--- /dev/null
+++ b/tests/DavidsonMixing/test.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+
+print("Test Davidson solver with mixing...")
+
+nargs=len(sys.argv)
+
+mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3]
+for i in range(4,nargs-4):
+  mpicmd = mpicmd + " "+sys.argv[i]
+print("MPI run command: {}".format(mpicmd)) 
+
+exe = sys.argv[nargs-4]
+inp = sys.argv[nargs-3]
+coords = sys.argv[nargs-2]
+print("coordinates file: %s"%coords)
+
+#create links to potentials files
+dst1 = 'pseudo.H_ONCV_PBE_SG15'
+dst2 = 'pseudo.Si_ONCV_PBE_SG15'
+src1 = sys.argv[nargs-1] + '/' + dst1
+src2 = sys.argv[nargs-1] + '/' + dst2
+
+cwd = os.getcwd()
+if not os.path.exists(cwd+'/'+dst1):
+  print("Create link to %s"%dst1)
+  os.symlink(src1, dst1)
+if not os.path.exists(cwd+'/'+dst2):
+  print("Create link to %s"%dst2)
+  os.symlink(src2, dst2)
+
+#run mgmol
+command = "{} {} -c {} -i {}".format(mpicmd,exe,inp,coords)
+print("Run command: {}".format(command))
+
+output = subprocess.check_output(command,shell=True)
+
+#analyse mgmol standard output
+#make sure force is below tolerance
+lines=output.split(b'\n')
+
+convergence=0
+for line in lines:
+  if line.count(b'DavidsonSolver') and line.count(b'convergence'):
+    convergence=1
+    break
+
+if convergence==0:
+  print("DavidsonSolver did not converge")
+  sys.exit(1)
+
+tol = 2.e-4
+energies=[]
+print("Check forces are smaller than tol = {}".format(tol))
+for line in lines:
+  if line.count(b'%%'):
+    print(line)
+    words=line.split()
+    energy=(words[5].split(b','))[0]
+    energies.append(energy)
+  if line.count(b'##'):
+    words=line.split()
+    if len(words)==8:
+      print(line)
+      for i in range(5,8):
+        if abs(eval(words[i]))>tol:
+          print("Error on force component {}".format(i-5))
+          sys.exit(1)
+
+niterations = len(energies)
+print("Davidson solver ran for {} iterations".format(niterations))
+if niterations>50:
+  print("Davidson test FAILED for taking too many iterations")
+  sys.exit(1)
+
+print("Check energy...")
+ref_energy=-6.26479
+last_energy = eval(energies[-1])
+print("Energy = {}".format(last_energy))
+if last_energy>ref_energy:
+  print("Last energy = {}".format(last_energy))
+  sys.exit(1)
+
+print("Test PASSED")
+sys.exit(0)