diff --git a/src/cas_based_on_top/NEED b/src/cas_based_on_top/NEED
new file mode 100644
index 00000000..d9a14e4a
--- /dev/null
+++ b/src/cas_based_on_top/NEED
@@ -0,0 +1,2 @@
+two_body_rdm
+dft_utils_in_r
diff --git a/src/cas_based_on_top/README.rst b/src/cas_based_on_top/README.rst
new file mode 100644
index 00000000..3b8a1897
--- /dev/null
+++ b/src/cas_based_on_top/README.rst
@@ -0,0 +1,9 @@
+==============
+on_top_density
+==============
+
+This plugin proposes different routines/providers to compute the on-top pair density of a CAS-based wave function. 
+
+This means that all determinants in psi_det must belong to an active-space. 
+
+As usual, see the file "example.irp.f" to get introduced to the main providers/routines.
diff --git a/src/cas_based_on_top/c_i_a_v_mos.irp.f b/src/cas_based_on_top/c_i_a_v_mos.irp.f
new file mode 100644
index 00000000..58192f5a
--- /dev/null
+++ b/src/cas_based_on_top/c_i_a_v_mos.irp.f
@@ -0,0 +1,125 @@
+
+ BEGIN_PROVIDER[double precision, core_mos_in_r_array, (n_core_orb,n_points_final_grid)]
+&BEGIN_PROVIDER[double precision, core_mos_in_r_array_transp,(n_points_final_grid,n_core_orb)]
+ implicit none
+ BEGIN_DOC
+! all COREE  MOs on the grid points, arranged in two different ways
+ END_DOC
+ integer :: i,j,k
+ do i = 1, n_core_orb
+  j = list_core(i) 
+  do k = 1, n_points_final_grid
+   core_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
+  enddo
+ enddo
+
+ do k = 1, n_points_final_grid
+  do i = 1, n_core_orb
+   core_mos_in_r_array(i,k) = core_mos_in_r_array_transp(k,i)
+  enddo
+ enddo
+
+END_PROVIDER 
+
+
+ BEGIN_PROVIDER[double precision, inact_mos_in_r_array, (n_inact_orb,n_points_final_grid)]
+&BEGIN_PROVIDER[double precision, inact_mos_in_r_array_transp,(n_points_final_grid,n_inact_orb)]
+ implicit none
+ BEGIN_DOC
+! all INACTIVE MOs on the grid points, arranged in two different ways
+ END_DOC
+ integer :: i,j,k
+ do i = 1, n_inact_orb
+  j = list_inact(i) 
+  do k = 1, n_points_final_grid
+   inact_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
+  enddo
+ enddo
+
+ do k = 1, n_points_final_grid
+  do i = 1, n_inact_orb
+   inact_mos_in_r_array(i,k) = inact_mos_in_r_array_transp(k,i)
+  enddo
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER[double precision, act_mos_in_r_array, (n_act_orb,n_points_final_grid)]
+&BEGIN_PROVIDER[double precision, act_mos_in_r_array_transp,(n_points_final_grid,n_act_orb)]
+ implicit none
+ BEGIN_DOC
+! all ACTIVE MOs on the grid points, arranged in two different ways
+ END_DOC
+ integer :: i,j,k
+ do i = 1, n_act_orb
+  j = list_act(i) 
+  do k = 1, n_points_final_grid
+   act_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
+  enddo
+ enddo
+
+ do k = 1, n_points_final_grid
+  do i = 1, n_act_orb
+   act_mos_in_r_array(i,k) = act_mos_in_r_array_transp(k,i)
+  enddo
+ enddo
+
+END_PROVIDER 
+
+
+ BEGIN_PROVIDER[double precision, virt_mos_in_r_array, (n_virt_orb,n_points_final_grid)]
+&BEGIN_PROVIDER[double precision, virt_mos_in_r_array_transp,(n_points_final_grid,n_virt_orb)]
+ implicit none
+ BEGIN_DOC
+! all VIRTUAL MOs on the grid points, arranged in two different ways
+ END_DOC
+ integer :: i,j,k
+ do i = 1, n_virt_orb
+  j = list_virt(i) 
+  do k = 1, n_points_final_grid
+   virt_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
+  enddo
+ enddo
+
+ do k = 1, n_points_final_grid
+  do i = 1, n_virt_orb
+   virt_mos_in_r_array(i,k) = virt_mos_in_r_array_transp(k,i)
+  enddo
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER[double precision, core_inact_act_mos_in_r_array, (n_core_inact_act_orb,n_points_final_grid)]
+&BEGIN_PROVIDER[double precision, core_inact_act_mos_in_r_array_transp,(n_points_final_grid,n_core_inact_act_orb)]
+ implicit none
+ integer :: i,j,k
+ do i = 1, n_core_inact_act_orb
+  j = list_core_inact_act(i) 
+  do k = 1, n_points_final_grid
+   core_inact_act_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
+  enddo
+ enddo
+
+ do k = 1, n_points_final_grid
+  do i = 1, n_core_inact_act_orb
+   core_inact_act_mos_in_r_array(i,k) = core_inact_act_mos_in_r_array_transp(k,i)
+  enddo
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER[double precision, core_inact_act_mos_grad_in_r_array, (3,n_core_inact_act_orb,n_points_final_grid)]
+ implicit none
+ integer :: i,j,k,l
+ do i = 1, n_core_inact_act_orb
+  j = list_core_inact_act(i) 
+  do k = 1, n_points_final_grid
+   do l = 1, 3
+    core_inact_act_mos_grad_in_r_array(l,i,k) = mos_grad_in_r_array(j,k,l)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+
diff --git a/src/cas_based_on_top/cas_based_density.irp.f b/src/cas_based_on_top/cas_based_density.irp.f
new file mode 100644
index 00000000..bec4da15
--- /dev/null
+++ b/src/cas_based_on_top/cas_based_density.irp.f
@@ -0,0 +1,47 @@
+program cas_based_density
+  implicit none
+  BEGIN_DOC
+! TODO : Small example to use the different quantities in this plugin
+  END_DOC
+
+  !! You force QP2 to read the wave function in the EZFIO folder
+  !! It is assumed that all Slater determinants in the wave function 
+  !! belongs to an active space defined by core, inactive and active list of orbitals 
+  read_wf = .True.
+  touch read_wf 
+  call routine_test_cas_based_density
+  
+end
+
+subroutine routine_test_cas_based_density
+ implicit none
+ integer :: ipoint, istate
+ double precision :: accu_n_elec(N_states),accu_n_elec_2(N_states)
+
+ ! PROVIDERS 
+ accu_n_elec = 0.d0
+ do istate = 1, N_states
+  do ipoint = 1, n_points_final_grid
+   accu_n_elec(istate) += one_e_cas_total_density(ipoint,istate) * final_weight_at_r_vector(ipoint)
+  enddo
+  print*,'istate = ',istate
+  print*,'accu_n_elec = ',accu_n_elec(istate)
+ enddo
+
+ ! ROUTINES 
+ double precision :: r(3),core_dens,inact_dens,act_dens(2,N_states),total_cas_dens(N_states)
+ accu_n_elec_2 = 0.d0
+ do ipoint = 1, n_points_final_grid
+  r(:) = final_grid_points(:,ipoint)
+  call give_cas_density_in_r(core_dens,inact_dens,act_dens,total_cas_dens,r)
+  do istate = 1, N_states
+   accu_n_elec_2(istate) += total_cas_dens(istate) * final_weight_at_r_vector(ipoint)
+  enddo
+ enddo
+
+ do istate = 1, N_states
+  print*,'istate = ',istate
+  print*,'accu_n_elec = ',accu_n_elec_2(istate)
+ enddo
+
+end
diff --git a/src/cas_based_on_top/cas_based_on_top.irp.f b/src/cas_based_on_top/cas_based_on_top.irp.f
new file mode 100644
index 00000000..89516336
--- /dev/null
+++ b/src/cas_based_on_top/cas_based_on_top.irp.f
@@ -0,0 +1,19 @@
+program cas_based_on_top_density
+  implicit none
+  BEGIN_DOC
+! TODO : Small example to use the different quantities in this plugin
+  END_DOC
+
+  !! You force QP2 to read the wave function in the EZFIO folder
+  !! It is assumed that all Slater determinants in the wave function 
+  !! belongs to an active space defined by core, inactive and active list of orbitals 
+  read_wf = .True.
+  touch read_wf
+!  call routine_test_cas_based_on_top_density
+ call routine
+end
+
+subroutine routine
+ implicit none
+ provide total_cas_on_top_density
+end
diff --git a/src/cas_based_on_top/cas_dens_prov.irp.f b/src/cas_based_on_top/cas_dens_prov.irp.f
new file mode 100644
index 00000000..453356b3
--- /dev/null
+++ b/src/cas_based_on_top/cas_dens_prov.irp.f
@@ -0,0 +1,99 @@
+
+ BEGIN_PROVIDER [double precision, one_e_cas_total_density ,(n_points_final_grid,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! one_e_cas_total_density = TOTAL DENSITY FOR a CAS wave function 
+ !
+ ! WARNING : if "no_core_density" == .True. then the core part of density is ignored 
+ END_DOC
+ integer :: ipoint,i,j,istate
+ do istate = 1, N_states
+  do ipoint = 1, n_points_final_grid
+   one_e_cas_total_density(ipoint,istate)   = one_e_act_density_alpha(ipoint,istate) + one_e_act_density_beta(ipoint,istate) & 
+                                            + 2.d0 *  inact_density(ipoint) 
+   if(.not.no_core_density)then !!! YOU ADD THE CORE DENSITY 
+    one_e_cas_total_density(ipoint,istate) += 2.d0 * core_density(ipoint) 
+   endif
+  enddo
+ enddo
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, one_e_act_density_alpha,(n_points_final_grid,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! one_e_act_density_alpha = pure ACTIVE part of the DENSITY for ALPHA ELECTRONS 
+ END_DOC
+ one_e_act_density_alpha = 0.d0
+ integer :: ipoint,i,j,istate
+ do istate = 1, N_states
+  do ipoint = 1, n_points_final_grid
+   do i = 1, n_act_orb
+    do j = 1, n_act_orb
+     one_e_act_density_alpha(ipoint,istate) += one_e_act_dm_alpha_mo_for_dft(j,i,istate) * act_mos_in_r_array(j,ipoint) * act_mos_in_r_array(i,ipoint)
+    enddo
+   enddo
+  enddo
+ enddo
+ END_PROVIDER 
+
+
+
+
+ BEGIN_PROVIDER [double precision, one_e_act_density_beta,(n_points_final_grid,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! one_e_act_density_beta = pure ACTIVE part of the DENSITY for BETA ELECTRONS 
+ END_DOC
+ one_e_act_density_beta = 0.d0
+ integer :: ipoint,i,j,istate
+ do istate = 1, N_states
+  do ipoint = 1, n_points_final_grid
+   do i = 1, n_act_orb
+    do j = 1, n_act_orb
+     one_e_act_density_beta(ipoint,istate) += one_e_act_dm_beta_mo_for_dft(j,i,istate) * act_mos_in_r_array(j,ipoint) * act_mos_in_r_array(i,ipoint)
+    enddo
+   enddo
+  enddo
+ enddo
+ END_PROVIDER 
+
+
+
+ BEGIN_PROVIDER [double precision, inact_density, (n_points_final_grid) ]
+ implicit none
+ BEGIN_DOC
+! INACTIVE part of the density for alpha/beta. 
+!
+! WARNING :: IF YOU NEED THE TOTAL DENSITY COMING FROM THE INACTIVE, 
+!
+!            YOU MUST MULTIPLY BY TWO 
+ END_DOC
+ integer :: i,j
+ inact_density = 0.d0
+ do i = 1, n_points_final_grid
+  do j = 1, n_inact_orb
+   inact_density(i) += inact_mos_in_r_array(j,i) **2
+  enddo
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, core_density, (n_points_final_grid) ]
+ implicit none
+ BEGIN_DOC
+! CORE part of the density for alpha/beta. 
+!
+! WARNING :: IF YOU NEED THE TOTAL DENSITY COMING FROM THE CORE, 
+!
+!            YOU MUST MULTIPLY BY TWO 
+ END_DOC
+ integer :: i,j
+ core_density = 0.d0
+ do i = 1, n_points_final_grid
+  do j = 1, n_core_orb
+   core_density(i) += core_mos_in_r_array(j,i) **2
+  enddo
+ enddo
+
+END_PROVIDER 
+
diff --git a/src/cas_based_on_top/cas_dens_rout.irp.f b/src/cas_based_on_top/cas_dens_rout.irp.f
new file mode 100644
index 00000000..7e7d8d99
--- /dev/null
+++ b/src/cas_based_on_top/cas_dens_rout.irp.f
@@ -0,0 +1,65 @@
+subroutine give_cas_density_in_r(core_dens,inact_dens,act_dens,total_cas_dens,r)
+ implicit none
+ BEGIN_DOC
+ ! returns the different component of the density at grid point r(3) for a CAS wave function
+ !
+ ! core_dens      : density coming from the CORE orbitals
+ !                
+ ! inact_dens     : density coming from the INACT orbitals
+ !                
+ ! act_dens(1/2,1:N_states)  : active part of the alpha/beta electrons for all states
+ !
+ ! total_cas_dens : total density of the cas wave function 
+ !
+ ! WARNING : if "no_core_density" == .True. then the core part of density is ignored in total_cas_dens
+ END_DOC
+
+ double precision, intent(in)  :: r(3)
+ double precision, intent(out) :: core_dens, inact_dens, act_dens(2,N_states), total_cas_dens(N_states)
+
+ double precision, allocatable :: mos_array(:),act_mos(:)
+ allocate(mos_array(mo_num))
+ call give_all_mos_at_r(r,mos_array)
+ integer :: i,iorb,j,jorb,istate
+
+ ! core part of the density 
+ core_dens = 0.d0
+ do i = 1, n_core_orb
+  iorb = list_core(i)
+  core_dens += mos_array(iorb)*mos_array(iorb)
+ enddo
+ core_dens = core_dens * 2.d0
+
+ ! inactive part of the density 
+ inact_dens = 0.d0
+ do i = 1, n_inact_orb
+  iorb = list_inact(i)
+  inact_dens += mos_array(iorb)*mos_array(iorb)
+ enddo
+ inact_dens = inact_dens * 2.d0
+ 
+ allocate(act_mos(n_act_orb))
+ do i = 1, n_act_orb
+  iorb = list_act(i)
+  act_mos(i) = mos_array(iorb)
+ enddo
+
+ ! active part of the density for alpha/beta and all states
+ act_dens = 0.d0
+ do istate = 1, N_states
+  do i = 1, n_act_orb
+   do j = 1, n_act_orb
+    act_dens(1,istate) += one_e_act_dm_alpha_mo_for_dft(j,i,istate) * act_mos(j) * act_mos(i)
+    act_dens(2,istate) += one_e_act_dm_beta_mo_for_dft(j,i,istate) * act_mos(j) * act_mos(i)
+   enddo
+  enddo
+ enddo
+ 
+ ! TOTAL density for all states
+ do istate = 1, N_states
+  total_cas_dens(istate) = inact_dens + act_dens(1,istate) + act_dens(2,istate) 
+  if(.not.no_core_density)then !!! YOU ADD THE CORE DENSITY 
+   total_cas_dens(istate) += core_dens
+  endif
+ enddo
+end
diff --git a/src/cas_based_on_top/cas_one_e_rdm.irp.f b/src/cas_based_on_top/cas_one_e_rdm.irp.f
new file mode 100644
index 00000000..df4ca271
--- /dev/null
+++ b/src/cas_based_on_top/cas_one_e_rdm.irp.f
@@ -0,0 +1,37 @@
+
+ BEGIN_PROVIDER [double precision, one_e_act_dm_beta_mo_for_dft, (n_act_orb,n_act_orb,N_states)]
+ implicit none
+ BEGIN_DOC
+ ! one_e_act_dm_beta_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the BETA ELECTRONS 
+ END_DOC
+ integer :: i,j,ii,jj,istate
+ do istate = 1, N_states
+  do ii = 1, n_act_orb
+   i = list_act(ii)
+   do jj = 1, n_act_orb
+    j = list_act(jj)
+    one_e_act_dm_beta_mo_for_dft(jj,ii,istate) = one_e_dm_mo_beta_for_dft(j,i,istate)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, one_e_act_dm_alpha_mo_for_dft, (n_act_orb,n_act_orb,N_states)]
+ implicit none
+ BEGIN_DOC
+ ! one_e_act_dm_alpha_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the ALPHA ELECTRONS 
+ END_DOC
+ integer :: i,j,ii,jj,istate
+ do istate = 1, N_states
+  do ii = 1, n_act_orb
+   i = list_act(ii)
+   do jj = 1, n_act_orb
+    j = list_act(jj)
+    one_e_act_dm_alpha_mo_for_dft(jj,ii,istate) = one_e_dm_mo_alpha_for_dft(j,i,istate)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
diff --git a/src/cas_based_on_top/eff_spin_dens.irp.f b/src/cas_based_on_top/eff_spin_dens.irp.f
new file mode 100644
index 00000000..8a70846a
--- /dev/null
+++ b/src/cas_based_on_top/eff_spin_dens.irp.f
@@ -0,0 +1,68 @@
+ BEGIN_PROVIDER [double precision, effective_spin_dm, (n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, grad_effective_spin_dm, (3,n_points_final_grid,N_states) ]
+ implicit none
+ BEGIN_DOC
+! effective_spin_dm(r_i) = \sqrt( n(r)^2 - 4 * ontop(r) )
+! effective spin density obtained from the total density and on-top pair density
+! see equation (6) of Phys. Chem. Chem. Phys., 2015, 17, 22412--22422 | 22413
+ END_DOC
+ provide total_cas_on_top_density 
+ integer :: i_point,i_state,i
+ double precision :: n2,m2,thr
+ thr = 1.d-14
+ effective_spin_dm = 0.d0
+ grad_effective_spin_dm = 0.d0
+ do i_state = 1, N_states
+  do i_point = 1, n_points_final_grid
+   n2 = (one_e_dm_and_grad_alpha_in_r(4,i_point,i_state)  + one_e_dm_and_grad_beta_in_r(4,i_point,i_state))
+   ! density squared 
+   n2 = n2 * n2 
+   if(n2 - 4.D0 * total_cas_on_top_density(i_point,i_state).gt.thr)then
+    effective_spin_dm(i_point,i_state) = dsqrt(n2 - 4.D0 * total_cas_on_top_density(i_point,i_state))
+    if(isnan(effective_spin_dm(i_point,i_state)))then
+     print*,'isnan(effective_spin_dm(i_point,i_state)' 
+     stop
+    endif
+    m2 = effective_spin_dm(i_point,i_state) 
+    m2 = 0.5d0 / m2 ! 1/(2 * sqrt(n(r)^2 - 4 * ontop(r)) )
+    do i = 1, 3
+     grad_effective_spin_dm(i,i_point,i_state) = m2 * ( one_e_stuff_for_pbe(i,i_point,i_state) - 4.d0 * grad_total_cas_on_top_density(i,i_point,i_state) )
+    enddo
+   else
+    effective_spin_dm(i_point,i_state) = 0.d0
+    grad_effective_spin_dm(:,i_point,i_state) = 0.d0
+   endif
+  enddo
+ enddo
+
+END_PROVIDER 
+
+
+ BEGIN_PROVIDER [double precision, effective_alpha_dm, (n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, effective_beta_dm, (n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, grad_effective_alpha_dm, (3,n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, grad_effective_beta_dm, (3,n_points_final_grid,N_states) ]
+ implicit none
+ BEGIN_DOC
+! effective_alpha_dm(r_i) = 1/2 * (effective_spin_dm(r_i) + n(r_i))
+! effective_beta_dm(r_i) = 1/2 * (-effective_spin_dm(r_i) + n(r_i))
+ END_DOC
+ provide total_cas_on_top_density 
+ integer :: i_point,i_state,i
+ double precision :: n,grad_n
+ do i_state = 1, N_states
+  do i_point = 1, n_points_final_grid
+   n = (one_e_dm_and_grad_alpha_in_r(4,i_point,i_state)  + one_e_dm_and_grad_beta_in_r(4,i_point,i_state))
+   effective_alpha_dm(i_point,i_state) = 0.5d0 * (n + effective_spin_dm(i_point,i_state))
+   effective_beta_dm(i_point,i_state) = 0.5d0 * (n - effective_spin_dm(i_point,i_state))
+   do i = 1, 3
+    grad_n = (one_e_dm_and_grad_alpha_in_r(i,i_point,i_state)  + one_e_dm_and_grad_beta_in_r(i,i_point,i_state)) 
+    grad_effective_alpha_dm(i,i_point,i_state) = 0.5d0 * (grad_n + grad_effective_spin_dm(i,i_point,i_state) )
+    grad_effective_beta_dm(i,i_point,i_state) = 0.5d0 *  (grad_n - grad_effective_spin_dm(i,i_point,i_state) )
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+
diff --git a/src/cas_based_on_top/example.irp.f b/src/cas_based_on_top/example.irp.f
new file mode 100644
index 00000000..f00956e3
--- /dev/null
+++ b/src/cas_based_on_top/example.irp.f
@@ -0,0 +1,41 @@
+
+subroutine write_on_top_in_real_space
+ implicit none
+ BEGIN_DOC
+! This routines is a simple example of how to plot the on-top pair density on a simple 1D grid 
+ END_DOC
+ double precision :: zmax,dz,r(3),on_top_in_r,total_density,zcenter,dist
+
+ integer :: nz,i,istate
+ character*(128)                :: output
+ integer                        :: i_unit_output,getUnitAndOpen
+ PROVIDE ezfio_filename
+ output=trim(ezfio_filename)//'.on_top'
+ print*,'output = ',trim(output)                                                                                         
+ i_unit_output = getUnitAndOpen(output,'w')
+
+
+ zmax = 2.0d0
+ print*,'nucl_coord(1,3) = ',nucl_coord(1,3)
+ print*,'nucl_coord(2,3) = ',nucl_coord(2,3)
+ dist = dabs(nucl_coord(1,3) - nucl_coord(2,3))
+ zmax += dist 
+ zcenter = (nucl_coord(1,3) + nucl_coord(2,3))*0.5d0
+ print*,'zcenter = ',zcenter
+ print*,'zmax    = ',zmax
+ nz = 1000
+ dz = zmax / dble(nz)
+ r(:) = 0.d0 
+ r(3) = zcenter -zmax * 0.5d0 
+ print*,'r(3)    = ',r(3)
+ istate = 1
+ do i = 1, nz
+  call give_on_top_in_r_one_state(r,istate,on_top_in_r)
+  call give_cas_density_in_r(r,total_density)
+  write(i_unit_output,*)r(3),on_top_in_r,total_density
+  r(3) += dz
+ enddo
+
+
+end
+
diff --git a/src/cas_based_on_top/on_top_cas_prov.irp.f b/src/cas_based_on_top/on_top_cas_prov.irp.f
new file mode 100644
index 00000000..18a7b7d2
--- /dev/null
+++ b/src/cas_based_on_top/on_top_cas_prov.irp.f
@@ -0,0 +1,76 @@
+
+ subroutine act_on_top_on_grid_pt(ipoint,istate,pure_act_on_top_of_r)
+ implicit none
+ BEGIN_DOC
+ ! act_on_top_on_grid_pt returns the purely ACTIVE part of the on top pair density
+ !
+ ! at the grid point ipoint, for the state istate
+ END_DOC
+ integer, intent(in) :: ipoint,istate
+ double precision, intent(out) :: pure_act_on_top_of_r
+ double precision :: phi_i,phi_j,phi_k,phi_l
+ integer :: i,j,k,l
+ ASSERT (istate <= N_states)
+
+ pure_act_on_top_of_r = 0.d0
+  do l = 1, n_act_orb
+   phi_l = act_mos_in_r_array(l,ipoint)
+   do k = 1, n_act_orb
+    phi_k = act_mos_in_r_array(k,ipoint)
+     do j = 1, n_act_orb
+      phi_j = act_mos_in_r_array(j,ipoint)
+      do i = 1, n_act_orb
+       phi_i = act_mos_in_r_array(i,ipoint)
+       !                                       1 2 1 2
+       pure_act_on_top_of_r += act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i * phi_j * phi_k * phi_l
+     enddo
+    enddo
+   enddo
+  enddo
+ end
+
+
+ BEGIN_PROVIDER [double precision, total_cas_on_top_density,(n_points_final_grid,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! on top pair density :: n2(r,r) at each of the Becke's grid point of a CAS-BASED wf
+ !
+ ! Contains all core/inact/act contribution.
+ !
+ ! !!!!! WARNING !!!!! If no_core_density then you REMOVE ALL CONTRIBUTIONS COMING FROM THE CORE ORBITALS 
+ END_DOC
+ integer :: i_point,istate
+ double precision :: wall_0,wall_1,core_inact_dm,pure_act_on_top_of_r
+ logical :: no_core
+ print*,'providing the total_cas_on_top_density'
+ ! for parallelization
+ provide inact_density core_density one_e_act_density_beta one_e_act_density_alpha act_mos_in_r_array
+ i_point = 1
+ istate = 1
+ call act_on_top_on_grid_pt(i_point,istate,pure_act_on_top_of_r)
+ call wall_time(wall_0)
+ if(no_core_density)then
+  print*,'USING THE VALENCE ONLY TWO BODY DENSITY'
+ endif
+ do istate = 1, N_states
+  !$OMP PARALLEL DO &
+  !$OMP DEFAULT (NONE)  &
+  !$OMP PRIVATE (i_point,core_inact_dm,pure_act_on_top_of_r) & 
+  !$OMP SHARED(total_cas_on_top_density,n_points_final_grid,inact_density,core_density,one_e_act_density_beta,one_e_act_density_alpha,no_core_density,istate)
+  do i_point = 1, n_points_final_grid
+    call act_on_top_on_grid_pt(i_point,istate,pure_act_on_top_of_r)
+    if(no_core_density) then
+     core_inact_dm = inact_density(i_point) 
+    else 
+     core_inact_dm = (inact_density(i_point) + core_density(i_point))
+    endif
+    total_cas_on_top_density(i_point,istate) = pure_act_on_top_of_r + core_inact_dm * (one_e_act_density_beta(i_point,istate) + one_e_act_density_alpha(i_point,istate)) + core_inact_dm*core_inact_dm
+   enddo
+  !$OMP END PARALLEL DO
+ enddo
+ call wall_time(wall_1)
+ print*,'provided the total_cas_on_top_density'
+ print*,'Time to provide :',wall_1 - wall_0
+
+ END_PROVIDER 
+
diff --git a/src/cas_based_on_top/on_top_cas_rout.irp.f b/src/cas_based_on_top/on_top_cas_rout.irp.f
new file mode 100644
index 00000000..6a4c15b9
--- /dev/null
+++ b/src/cas_based_on_top/on_top_cas_rout.irp.f
@@ -0,0 +1,118 @@
+
+subroutine give_core_inact_act_density_in_r(r,mos_array,core_density_in_r,inact_density_in_r,act_density_in_r, total_density)
+ implicit none
+ double precision, intent(in) :: r(3),mos_array(mo_num)
+ double precision, intent(out):: core_density_in_r,inact_density_in_r,act_density_in_r(2,N_states),total_density(N_states)
+ BEGIN_DOC
+! core, inactive and active part of the density for alpha/beta electrons 
+! 
+! the density coming from the core and inactive are the same for alpha/beta electrons
+!
+! act_density(1/2, i) = alpha/beta density for the ith state
+!
+! total_density(i) = 2 * (core_density_in_r+inact_density_in_r) + act_density_in_r(1,i) + act_density_in_r(2,i)
+ END_DOC
+ integer :: i,j,istate
+
+ core_density_in_r = 0.d0
+ do i = 1, n_core_orb
+  j = list_core(i)
+  core_density_in_r += mos_array(j) * mos_array(j) 
+ enddo
+
+ inact_density_in_r = 0.d0
+ do i = 1, n_inact_orb
+  j = list_inact(i)
+  inact_density_in_r += mos_array(j) * mos_array(j)
+ enddo
+
+ double precision, allocatable :: act_mos(:)
+ double precision :: tmp
+ allocate(act_mos(n_act_orb))
+ do i = 1, n_act_orb
+  j = list_act(i)
+  act_mos(i) = mos_array(j)
+ enddo
+
+ act_density_in_r = 0.d0
+ do istate = 1, N_states
+  do i = 1, n_act_orb
+   do j = 1, n_act_orb
+    tmp = act_mos(i) * act_mos(j) 
+    act_density_in_r(1,istate) += tmp * one_e_act_dm_alpha_mo_for_dft(j,i,istate) 
+    act_density_in_r(2,istate) += tmp * one_e_act_dm_beta_mo_for_dft(j,i,istate) 
+   enddo
+  enddo
+  total_density(istate) = 2.d0 * (core_density_in_r + inact_density_in_r) + act_density_in_r(1,istate) + act_density_in_r(2,istate)
+ enddo
+
+end
+
+ subroutine give_active_on_top_in_r_one_state(r,istate,mos_array,act_on_top)
+ implicit none
+ BEGIN_DOC
+ ! gives the purely active on-top pair density for a given state 
+ END_DOC
+ integer, intent(in) :: istate
+ double precision, intent(in) :: r(3),mos_array(mo_num)
+ double precision, intent(out) :: act_on_top
+ double precision :: phi_i,phi_j,phi_k,phi_l
+ integer :: i,j,k,l
+
+ double precision, allocatable :: act_mos(:)
+ double precision :: tmp
+ allocate(act_mos(n_act_orb))
+ do i = 1, n_act_orb
+  j = list_act(i)
+  act_mos(i) = mos_array(j)
+ enddo
+
+ act_on_top = 0.d0
+  do l = 1, n_act_orb
+   phi_l = act_mos(l)
+   do k = 1, n_act_orb
+    phi_k = act_mos(k)
+     do j = 1, n_act_orb
+      phi_j = act_mos(j)
+      tmp = phi_l * phi_k * phi_j 
+      do i = 1, n_act_orb
+       phi_i = act_mos(i)
+       !                                                       1 2 1 2
+       act_on_top += act_2_rdm_ab_mo(i,j,k,l,istate) * tmp * phi_i
+     enddo
+    enddo
+   enddo
+  enddo
+ end
+
+ subroutine give_on_top_in_r_one_state(r,istate,on_top_in_r)
+ implicit none
+ integer, intent(in) :: istate
+ double precision, intent(in) :: r(3)
+ double precision, intent(out) :: on_top_in_r
+ BEGIN_DOC
+ ! on top pair density in r for the state istate a CAS-BASED wf 
+ !
+ ! note that if no_core_density .EQ. .True., all core contributions are excluded
+ END_DOC
+ double precision, allocatable :: mos_array(:)
+ provide act_2_rdm_ab_mo one_e_act_dm_alpha_mo_for_dft one_e_act_dm_beta_mo_for_dft
+ allocate(mos_array(mo_num))
+ call give_all_mos_at_r(r,mos_array)
+
+ double precision :: core_density_in_r, inact_density_in_r, act_density_in_r(2,N_states), total_density(N_states)
+ double precision :: act_on_top,core_inact_dm
+ ! getting the different part of the density in r
+ call give_core_inact_act_density_in_r(r,mos_array,core_density_in_r,inact_density_in_r,act_density_in_r, total_density)
+ ! getting the purely active part of the density in r
+ call give_active_on_top_in_r_one_state(r,istate,mos_array,act_on_top)
+
+ if(no_core_density) then
+  core_inact_dm = inact_density_in_r
+ else 
+  core_inact_dm = core_density_in_r + inact_density_in_r
+ endif
+ on_top_in_r = act_on_top + core_inact_dm * (act_density_in_r(1,istate) + act_density_in_r(2,istate)) + core_inact_dm*core_inact_dm
+
+ end
+
diff --git a/src/cas_based_on_top/on_top_grad.irp.f b/src/cas_based_on_top/on_top_grad.irp.f
new file mode 100644
index 00000000..4c581273
--- /dev/null
+++ b/src/cas_based_on_top/on_top_grad.irp.f
@@ -0,0 +1,73 @@
+ subroutine give_on_top_gradient(ipoint,istate,ontop_grad)
+ implicit none
+ BEGIN_DOC
+ ! on top pair density and its gradient evaluated at a given point of the grid 
+ ! ontop_grad(1:3) :: gradients of the on-top pair density 
+ ! ontop_grad(4)   :: on-top pair density 
+ END_DOC
+ double precision, intent(out) :: ontop_grad(4)
+ integer, intent(in) :: ipoint,istate
+ double precision :: phi_jkl,phi_ikl,phi_ijl,phi_ijk
+ integer :: i,j,k,l,m
+ 
+ ontop_grad = 0.d0
+ do l = 1, n_core_inact_act_orb
+  do k = 1, n_core_inact_act_orb
+    do j = 1, n_core_inact_act_orb
+     do i = 1, n_core_inact_act_orb
+      phi_jkl = core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(k,ipoint) * core_inact_act_mos_in_r_array(l,ipoint) 
+      phi_ikl = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(k,ipoint) * core_inact_act_mos_in_r_array(l,ipoint) 
+      phi_ijl = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(l,ipoint) 
+      phi_ijk = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(k,ipoint) 
+     !                                                                                          1 2 1 2 
+      ontop_grad(4) += phi_ijk * core_inact_act_mos_in_r_array(l,ipoint) * full_occ_2_rdm_ab_mo(i,j,k,l,istate)
+     do m = 1,3
+      ontop_grad (m) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * & 
+     ( core_inact_act_mos_grad_in_r_array(m,i,ipoint) * phi_jkl +  core_inact_act_mos_grad_in_r_array(m,j,ipoint) * phi_ikl + & 
+       core_inact_act_mos_grad_in_r_array(m,k,ipoint) * phi_ijl +  core_inact_act_mos_grad_in_r_array(m,l,ipoint) * phi_ijk )
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+ end
+
+ BEGIN_PROVIDER [double precision, grad_total_cas_on_top_density,(4,n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, wall_time_core_inact_act_on_top_of_r ]
+ implicit none
+ BEGIN_DOC
+ ! grad_total_cas_on_top_density(1:3,ipoint,istate) : provider for the on top pair density gradient (x,y,z) for the point 'ipoint' and state 'istate'
+ !
+ ! grad_total_cas_on_top_density(4,ipoint,istate)   : on top pair density for the point 'ipoint' and state 'istate'
+ END_DOC
+ integer :: i_point,i_state,i
+ double precision :: wall_0,wall_1
+ double precision :: core_inact_act_on_top_of_r_from_provider,ontop_grad(4)
+
+ print*,'providing the core_inact_act_on_top_of_r'
+ i_point = 1
+ provide full_occ_2_rdm_ab_mo 
+ i_state = 1
+ call give_on_top_gradient(i_point,i_state,ontop_grad)
+ call wall_time(wall_0)
+ !$OMP PARALLEL DO &
+ !$OMP DEFAULT (NONE)  &
+ !$OMP PRIVATE (i_point,i_state,ontop_grad) & 
+ !$OMP SHARED(grad_total_cas_on_top_density,n_points_final_grid,N_states)
+ do i_point = 1, n_points_final_grid
+  do i_state = 1, N_states
+   call give_on_top_gradient(i_point,i_state,ontop_grad)
+   do i = 1, 4
+    grad_total_cas_on_top_density(i,i_point,i_state) = ontop_grad(i)
+   enddo
+  enddo
+ enddo
+ !$OMP END PARALLEL DO
+ call wall_time(wall_1)
+ print*,'provided the core_inact_act_on_top_of_r'
+ print*,'Time to provide :',wall_1 - wall_0
+ wall_time_core_inact_act_on_top_of_r = wall_1 - wall_0
+
+ END_PROVIDER 
+
+