diff --git a/devel/qmcchem/EZFIO.cfg b/devel/qmcchem/EZFIO.cfg
new file mode 100644
index 0000000..c9d238f
--- /dev/null
+++ b/devel/qmcchem/EZFIO.cfg
@@ -0,0 +1,6 @@
+[ci_threshold]
+type: Threshold
+doc: Threshold on the CI coefficients as computed in QMCChem
+interface: ezfio,provider,ocaml
+default: 1.e-8
+
diff --git a/devel/qmcchem/NEED b/devel/qmcchem/NEED
new file mode 100644
index 0000000..afb53ca
--- /dev/null
+++ b/devel/qmcchem/NEED
@@ -0,0 +1,4 @@
+determinants
+davidson_undressed
+fci
+zmq
diff --git a/devel/qmcchem/README.rst b/devel/qmcchem/README.rst
new file mode 100644
index 0000000..0c72b8e
--- /dev/null
+++ b/devel/qmcchem/README.rst
@@ -0,0 +1,12 @@
+==============
+QmcChem Module
+==============
+
+For multi-state calculations, to extract state 2 use:
+
+``
+QP_STATE=2 qp_run save_for_qmcchem x.ezfio
+``
+(state 1 is the ground state).
+
+
diff --git a/devel/qmcchem/expand_spindets_qmcchem.irp.f b/devel/qmcchem/expand_spindets_qmcchem.irp.f
new file mode 100644
index 0000000..2e1894b
--- /dev/null
+++ b/devel/qmcchem/expand_spindets_qmcchem.irp.f
@@ -0,0 +1,8 @@
+program densify
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call generate_all_alpha_beta_det_products()
+ call diagonalize_ci
+ call save_wavefunction
+end
diff --git a/devel/qmcchem/pot_ao_pseudo_ints.irp.f b/devel/qmcchem/pot_ao_pseudo_ints.irp.f
new file mode 100644
index 0000000..97ac748
--- /dev/null
+++ b/devel/qmcchem/pot_ao_pseudo_ints.irp.f
@@ -0,0 +1,131 @@
+
+BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ]
+ implicit none
+ BEGIN_DOC
+! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
+!
+!
+ END_DOC
+ ! l,m : Y(l,m) parameters
+ ! c(3) : pseudopotential center
+ ! a(3) : Atomic Orbital center
+ ! n_a(3) : Powers of x,y,z in the Atomic Orbital
+ ! g_a : Atomic Orbital exponent
+ ! r : Distance between the Atomic Orbital center and the considered point
+ double precision, external :: ylm_orb
+ integer :: n_a(3)
+ double precision :: a(3), c(3), g_a
+ integer :: i,j,k,l,m,n,p
+ double precision :: dr, Ulc
+ double precision :: y
+ double precision, allocatable :: r(:)
+
+ allocate (r(pseudo_grid_size))
+ dr = pseudo_grid_rmax/dble(pseudo_grid_size)
+ r(1) = 0.d0
+ do j=2,pseudo_grid_size
+ r(j) = r(j-1) + dr
+ enddo
+
+ ao_pseudo_grid = 0.d0
+ do j=1,pseudo_grid_size
+ do k=1,nucl_num
+ c(1:3) = nucl_coord(k,1:3)
+ do l=0,pseudo_lmax
+ do i=1,ao_num
+ a(1:3) = nucl_coord(ao_nucl(i),1:3)
+ n_a(1:3) = ao_power(i,1:3)
+ do p=1,ao_prim_num(i)
+ g_a = ao_expo_ordered_transp(p,i)
+ do m=-l,l
+ y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
+ ao_pseudo_grid(i,m,l,k,j) = ao_pseudo_grid(i,m,l,k,j) + &
+ ao_coef_normalized_ordered_transp(p,i)*y
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ deallocate(r)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ]
+ implicit none
+ BEGIN_DOC
+! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
+!
+!
+ END_DOC
+ ! l,m : Y(l,m) parameters
+ ! c(3) : pseudopotential center
+ ! a(3) : Atomic Orbital center
+ ! n_a(3) : Powers of x,y,z in the Atomic Orbital
+ ! g_a : Atomic Orbital exponent
+ ! r : Distance between the Atomic Orbital center and the considered point
+ double precision, external :: ylm_orb
+ integer :: n_a(3)
+ double precision :: a(3), c(3), g_a
+ integer :: i,j,k,l,m,n,p
+ double precision :: dr, Ulc
+ double precision :: y
+ double precision, allocatable :: r(:)
+
+ allocate (r(pseudo_grid_size))
+
+ dr = pseudo_grid_rmax/dble(pseudo_grid_size)
+ r(1) = 0.d0
+ do j=2,pseudo_grid_size
+ r(j) = r(j-1) + dr
+ enddo
+
+ mo_pseudo_grid = 0.d0
+ do n=1,pseudo_grid_size
+ do k=1,nucl_num
+ do l=0,pseudo_lmax
+ do m=-l,l
+ do i=1,ao_num
+ do j=1,mo_num
+ if (dabs(ao_pseudo_grid(i,m,l,k,n)) < 1.e-12) then
+ cycle
+ endif
+ if (dabs(mo_coef(i,j)) < 1.e-8) then
+ cycle
+ endif
+ mo_pseudo_grid(j,m,l,k,n) = mo_pseudo_grid(j,m,l,k,n) + &
+ ao_pseudo_grid(i,m,l,k,n) * mo_coef(i,j)
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ deallocate(r)
+
+END_PROVIDER
+
+double precision function test_pseudo_grid_ao(i,j)
+ implicit none
+ integer, intent(in) :: i,j
+ integer :: k,l,m,n
+ double precision :: r, dr,u
+ dr = pseudo_grid_rmax/dble(pseudo_grid_size)
+
+ test_pseudo_grid_ao = 0.d0
+ r = 0.d0
+ do k=1,pseudo_grid_size
+ do n=1,nucl_num
+ do l = 0,pseudo_lmax
+ u = pseudo_v_kl(n,l,1) * exp(-pseudo_dz_kl(n,l,1)*r*r)* r*r*dr
+ do m=-l,l
+ test_pseudo_grid_ao += ao_pseudo_grid(i,m,l,n,k) * ao_pseudo_grid(j,m,l,n,k) * u
+ enddo
+ enddo
+ enddo
+ r = r+dr
+ enddo
+end
diff --git a/devel/qmcchem/pseudo.irp.f b/devel/qmcchem/pseudo.irp.f
new file mode 100644
index 0000000..1398100
--- /dev/null
+++ b/devel/qmcchem/pseudo.irp.f
@@ -0,0 +1,8 @@
+subroutine write_pseudopotential
+ implicit none
+ BEGIN_DOC
+! Write the pseudo_potential into the EZFIO file
+ END_DOC
+ call ezfio_set_pseudo_mo_pseudo_grid(mo_pseudo_grid)
+end
+
diff --git a/devel/qmcchem/qmc_e_curve.irp.f b/devel/qmcchem/qmc_e_curve.irp.f
new file mode 100644
index 0000000..c0a3b7d
--- /dev/null
+++ b/devel/qmcchem/qmc_e_curve.irp.f
@@ -0,0 +1,113 @@
+program e_curve
+ use bitmasks
+ implicit none
+ integer :: i,j,k, kk, nab, m, l
+ double precision :: norm, E, hij, num, ci, cj
+ integer, allocatable :: iorder(:)
+ double precision , allocatable :: norm_sort(:)
+ double precision :: e_0(N_states)
+ PROVIDE mo_two_e_integrals_in_map
+
+ nab = n_det_alpha_unique+n_det_beta_unique
+ allocate ( norm_sort(0:nab), iorder(0:nab) )
+
+ double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
+ double precision, allocatable :: u_0(:,:), v_0(:,:)
+ allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
+ allocate(u_0(N_states,N_det),v_0(N_states,N_det))
+
+
+
+ norm_sort(0) = 0.d0
+ iorder(0) = 0
+ do i=1,n_det_alpha_unique
+ norm_sort(i) = det_alpha_norm(i)
+ iorder(i) = i
+ enddo
+
+ do i=1,n_det_beta_unique
+ norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+ iorder(i+n_det_alpha_unique) = -i
+ enddo
+
+ call dsort(norm_sort(1),iorder(1),nab)
+
+ if (.not.read_wf) then
+ stop 'Please set read_wf to true'
+ endif
+
+ PROVIDE psi_bilinear_matrix_values nuclear_repulsion
+ print *, ''
+ print *, '=============================='
+ print *, 'Energies at different cut-offs'
+ print *, '=============================='
+ print *, ''
+ print *, '=========================================================='
+ print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)', 'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
+ print *, '=========================================================='
+ double precision :: thresh
+ integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+ thresh = 1.d-10
+ do j=0,nab
+ i = iorder(j)
+ if (i<0) then
+ do k=1,n_det
+ if (psi_bilinear_matrix_columns(k) == -i) then
+ psi_bilinear_matrix_values(k,1) = 0.d0
+ endif
+ enddo
+ else
+ do k=1,n_det
+ if (psi_bilinear_matrix_rows(k) == i) then
+ psi_bilinear_matrix_values(k,1) = 0.d0
+ endif
+ enddo
+ endif
+ if (thresh > norm_sort(j)) then
+ cycle
+ endif
+
+ u_0 = psi_bilinear_matrix_values(1:N_det,1:N_states)
+ v_t = 0.d0
+ s_t = 0.d0
+ call dtranspose( &
+ u_0, &
+ size(u_0, 1), &
+ u_t, &
+ size(u_t, 1), &
+ N_det, N_states)
+ call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
+ call dtranspose( &
+ v_t, &
+ size(v_t, 1), &
+ v_0, &
+ size(v_0, 1), &
+ N_states, N_det)
+
+ double precision, external :: u_dot_u, u_dot_v
+ do i=1,N_states
+ e_0(i) = u_dot_v(v_t(1,i),u_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
+ enddo
+
+ m = 0
+ do k=1,n_det
+ if (psi_bilinear_matrix_values(k,1) /= 0.d0) then
+ m = m+1
+ endif
+ enddo
+
+ if (m == 0) then
+ exit
+ endif
+ E = E_0(1) + nuclear_repulsion
+ norm = u_dot_u(u_0(1,1),N_det)
+ print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)', thresh, m, &
+ dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
+ dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, E
+ thresh = thresh * dsqrt(10.d0)
+ enddo
+ print *, '=========================================================='
+
+ deallocate (iorder, norm_sort)
+end
+
diff --git a/devel/qmcchem/save_for_qmcchem.irp.f b/devel/qmcchem/save_for_qmcchem.irp.f
new file mode 100644
index 0000000..c7bcf46
--- /dev/null
+++ b/devel/qmcchem/save_for_qmcchem.irp.f
@@ -0,0 +1,45 @@
+program save_for_qmc
+
+ integer :: iunit
+ logical :: exists
+ double precision :: e_ref
+
+ ! Determinants
+ read_wf = .True.
+ TOUCH read_wf
+ print *, "N_det = ", N_det
+ call write_spindeterminants
+
+ ! Reference Energy
+ if (do_pseudo) then
+ call write_pseudopotential
+ endif
+ call system( &
+ 'mkdir -p '//trim(ezfio_filename)//'/simulation ;' // &
+ 'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/simulation/.version ; ' // &
+ 'mkdir -p '//trim(ezfio_filename)//'/properties ;' // &
+ 'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/properties/.version ; ' // &
+ 'echo T > '//trim(ezfio_filename)//'/properties/e_loc' &
+ )
+ iunit = 13
+ open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write')
+ call ezfio_has_fci_energy_pt2(exists)
+ if (exists) then
+ call ezfio_get_fci_energy_pt2(e_ref)
+ else
+ call ezfio_has_fci_energy(exists)
+ if (exists) then
+ call ezfio_get_fci_energy(e_ref)
+ else
+ call ezfio_has_hartree_fock_energy(exists)
+ if (exists) then
+ call ezfio_get_hartree_fock_energy(e_ref)
+ else
+ e_ref = 0.d0
+ endif
+ endif
+ endif
+ write(iunit,*) e_ref
+ close(iunit)
+
+end
diff --git a/devel/qmcchem/truncate_wf_qmcchem.irp.f b/devel/qmcchem/truncate_wf_qmcchem.irp.f
new file mode 100644
index 0000000..bc33499
--- /dev/null
+++ b/devel/qmcchem/truncate_wf_qmcchem.irp.f
@@ -0,0 +1,116 @@
+program truncate
+ read_wf = .True.
+ SOFT_TOUCH read_wf
+ call run
+end
+
+subroutine run
+ use bitmasks
+ implicit none
+ integer :: i,j,k, kk, nab, m, l
+ double precision :: norm, E, hij, num, ci, cj
+ integer, allocatable :: iorder(:)
+ double precision , allocatable :: norm_sort(:)
+ double precision :: e_0(N_states)
+
+ PROVIDE mo_two_e_integrals_in_map H_apply_buffer_allocated
+
+ nab = n_det_alpha_unique+n_det_beta_unique
+ allocate ( norm_sort(0:nab), iorder(0:nab) )
+
+ integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+ double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
+ double precision, allocatable :: u_0(:,:), v_0(:,:)
+ allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
+ allocate(u_0(N_det,N_states),v_0(N_det,N_states))
+
+ norm_sort(0) = 0.d0
+ iorder(0) = 0
+ do i=1,n_det_alpha_unique
+ norm_sort(i) = det_alpha_norm(i)
+ iorder(i) = i
+ enddo
+
+ do i=1,n_det_beta_unique
+ norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+ iorder(i+n_det_alpha_unique) = -i
+ enddo
+
+ call dsort(norm_sort(1),iorder(1),nab)
+
+
+ PROVIDE psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns
+ PROVIDE nuclear_repulsion
+ print *, ''
+ do j=0,nab
+ i = iorder(j)
+ if (i<0) then
+ !$OMP PARALLEL DO PRIVATE(k)
+ do k=1,n_det
+ if (psi_bilinear_matrix_columns(k) == -i) then
+ do l=1,N_states
+ psi_bilinear_matrix_values(k,l) = 0.d0
+ enddo
+ endif
+ enddo
+ !$OMP END PARALLEL DO
+ else
+ !$OMP PARALLEL DO PRIVATE(k)
+ do k=1,n_det
+ if (psi_bilinear_matrix_rows(k) == i) then
+ do l=1,N_states
+ psi_bilinear_matrix_values(k,l) = 0.d0
+ enddo
+ endif
+ enddo
+ !$OMP END PARALLEL DO
+ endif
+ if (ci_threshold <= norm_sort(j)) then
+ exit
+ endif
+ enddo
+
+ m = 0
+ do k=1,n_det
+ if (sum(psi_bilinear_matrix_values(k,1:N_states)) /= 0.d0) then
+ m = m+1
+ endif
+ enddo
+
+ do k=1,N_states
+ E = E_0(k) + nuclear_repulsion
+ enddo
+ print *, 'Number of determinants:', m
+ call wf_of_psi_bilinear_matrix(.True.)
+ call save_wavefunction
+
+ u_0(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
+ v_0(1:N_det,1:N_states) = 0.d0
+ u_t(1:N_states,1:N_det) = 0.d0
+ v_t(1:N_states,1:N_det) = 0.d0
+ s_t(1:N_states,1:N_det) = 0.d0
+ call dtranspose( &
+ u_0, &
+ size(u_0, 1), &
+ u_t, &
+ size(u_t, 1), &
+ N_det, N_states)
+ print *, 'Computing H|Psi> ...'
+ call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
+ print *, 'Done'
+ call dtranspose( &
+ v_t, &
+ size(v_t, 1), &
+ v_0, &
+ size(v_0, 1), &
+ N_states, N_det)
+
+ double precision, external :: u_dot_u, u_dot_v
+ do i=1,N_states
+ e_0(i) = u_dot_v(u_0(1,i),v_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
+ print *, 'E(',i,') = ', e_0(i) + nuclear_repulsion
+ enddo
+
+ deallocate (iorder, norm_sort)
+end
+