diff --git a/src/mol_properties/EZFIO.cfg b/src/mol_properties/EZFIO.cfg
new file mode 100644
index 00000000..35a095fb
--- /dev/null
+++ b/src/mol_properties/EZFIO.cfg
@@ -0,0 +1,23 @@
+[print_all_transitions]
+type: logical
+doc: If true, print the transition between all the states
+interface: ezfio,provider,ocaml
+default: false
+
+[calc_dipole_moment]
+type: logical
+doc: If true, the electric dipole moment will be computed
+interface: ezfio,provider,ocaml
+default: false
+
+[calc_tr_dipole_moment]
+type: logical
+doc: If true and N_states > 1, the transition electric dipole moment will be computed
+interface: ezfio,provider,ocaml
+default: false
+
+[calc_osc_str]
+type: logical
+doc: If true and N_states > 1, the oscillator strength will be computed
+interface: ezfio,provider,ocaml
+default: false
diff --git a/src/mol_properties/NEED b/src/mol_properties/NEED
new file mode 100644
index 00000000..8d89a452
--- /dev/null
+++ b/src/mol_properties/NEED
@@ -0,0 +1,2 @@
+determinants
+davidson_undressed
diff --git a/src/mol_properties/README.md b/src/mol_properties/README.md
new file mode 100644
index 00000000..583e79ea
--- /dev/null
+++ b/src/mol_properties/README.md
@@ -0,0 +1,17 @@
+# Molecular properties
+
+Available quantities:  
+- Electric dipole moment  
+- Electric transition dipole moment  
+- Oscillator strength  
+
+They are not computed by default. To compute them:
+```
+qp set mol_properties calc_dipole_moment true  
+qp set mol_properties calc_tr_dipole_moment true  
+qp set mol_properties calc_osc_str true  
+```
+If you are interested in transitions between two excited states:  
+```
+qp set mol_properties print_all_transitions true
+```
diff --git a/src/mol_properties/ci_energy_no_diag.irp.f b/src/mol_properties/ci_energy_no_diag.irp.f
new file mode 100644
index 00000000..a4407d3b
--- /dev/null
+++ b/src/mol_properties/ci_energy_no_diag.irp.f
@@ -0,0 +1,13 @@
+BEGIN_PROVIDER [double precision,  ci_energy_no_diag, (N_states) ]
+
+  implicit none
+
+  BEGIN_DOC
+  ! CI energy from density matrices and integrals
+  ! Avoid the rediagonalization for ci_energy
+  END_DOC
+
+  ci_energy_no_diag = psi_energy + nuclear_repulsion
+
+END_PROVIDER
+
diff --git a/src/mol_properties/mo_deriv_1.irp.f b/src/mol_properties/mo_deriv_1.irp.f
new file mode 100644
index 00000000..cfe6f789
--- /dev/null
+++ b/src/mol_properties/mo_deriv_1.irp.f
@@ -0,0 +1,30 @@
+ BEGIN_PROVIDER [double precision, mo_deriv_1_x , (mo_num,mo_num)]
+&BEGIN_PROVIDER [double precision, mo_deriv_1_y , (mo_num,mo_num)]
+&BEGIN_PROVIDER [double precision, mo_deriv_1_z , (mo_num,mo_num)]
+ BEGIN_DOC
+ ! array of the integrals of MO_i * d/dx MO_j
+ ! array of the integrals of MO_i * d/dy MO_j
+ ! array of the integrals of MO_i * d/dz MO_j
+ END_DOC
+ implicit none
+
+  call ao_to_mo(                                                     &
+      ao_deriv_1_x,                                                   &
+      size(ao_deriv_1_x,1),                                           &
+      mo_deriv_1_x,                                                   &
+      size(mo_deriv_1_x,1)                                            &
+      )
+  call ao_to_mo(                                                     &
+      ao_deriv_1_y,                                                   &
+      size(ao_deriv_1_y,1),                                           &
+      mo_deriv_1_y,                                                   &
+      size(mo_deriv_1_y,1)                                            &
+      )
+  call ao_to_mo(                                                     &
+      ao_deriv_1_z,                                                   &
+      size(ao_deriv_1_z,1),                                           &
+      mo_deriv_1_z,                                                   &
+      size(mo_deriv_1_z,1)                                            &
+      )
+
+END_PROVIDER
diff --git a/src/mol_properties/multi_s_deriv_1.irp.f b/src/mol_properties/multi_s_deriv_1.irp.f
new file mode 100644
index 00000000..b30130b7
--- /dev/null
+++ b/src/mol_properties/multi_s_deriv_1.irp.f
@@ -0,0 +1,78 @@
+ BEGIN_PROVIDER [double precision, multi_s_deriv_1, (N_states, N_states)]
+&BEGIN_PROVIDER [double precision, multi_s_x_deriv_1, (N_states, N_states)]
+&BEGIN_PROVIDER [double precision, multi_s_y_deriv_1, (N_states, N_states)]
+&BEGIN_PROVIDER [double precision, multi_s_z_deriv_1, (N_states, N_states)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! Providers for :
+  ! <Psi_m|v_x|Psi_n>
+  ! <Psi_m|v_y|Psi_n>
+  ! <Psi_m|v_z|Psi_n>
+  ! ||v|| = sqrt(v_x^2 + v_y^2 + v_z^2)
+  ! v_x = d/dx
+  ! Cf. multi_s_dipole_moment for the equations
+  END_DOC
+
+  integer          :: istate,jstate ! States
+  integer          :: i,j           ! general spatial MOs
+  double precision :: nuclei_part_x, nuclei_part_y, nuclei_part_z
+ 
+  multi_s_x_deriv_1 = 0.d0
+  multi_s_y_deriv_1 = 0.d0
+  multi_s_z_deriv_1 = 0.d0
+ 
+  do jstate = 1, N_states
+    do istate = 1, N_states
+ 
+      do i = 1, mo_num  
+        ! Diag part
+        multi_s_x_deriv_1(istate,jstate) -= one_e_tr_dm_mo(i,i,istate,jstate) * mo_deriv_1_x(i,i)
+        multi_s_y_deriv_1(istate,jstate) -= one_e_tr_dm_mo(i,i,istate,jstate) * mo_deriv_1_y(i,i)
+        multi_s_z_deriv_1(istate,jstate) -= one_e_tr_dm_mo(i,i,istate,jstate) * mo_deriv_1_z(i,i)
+ 
+        do j = 1, mo_num  
+          if (i == j) then
+           cycle
+          endif
+          ! Extra diag part
+          multi_s_x_deriv_1(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_deriv_1_x(j,i)  
+          multi_s_y_deriv_1(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_deriv_1_y(j,i) 
+          multi_s_z_deriv_1(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_deriv_1_z(j,i) 
+        enddo
+      enddo
+ 
+    enddo
+  enddo
+ 
+  ! Nuclei part
+  nuclei_part_x = 0.d0
+  nuclei_part_y = 0.d0
+  nuclei_part_z = 0.d0
+ 
+  do i = 1,nucl_num 
+    nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) 
+    nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) 
+    nuclei_part_z += nucl_charge(i) * nucl_coord(i,3) 
+  enddo
+ 
+  ! Only if istate = jstate, otherwise 0 by the orthogonality of the states
+  do istate = 1, N_states
+    multi_s_x_deriv_1(istate,istate) += nuclei_part_x
+    multi_s_y_deriv_1(istate,istate) += nuclei_part_y
+    multi_s_z_deriv_1(istate,istate) += nuclei_part_z
+  enddo
+ 
+  ! d = <Psi|r|Psi>
+  do jstate = 1, N_states
+    do istate = 1, N_states
+      multi_s_deriv_1(istate,jstate) = &
+        dsqrt(multi_s_x_deriv_1(istate,jstate)**2 & 
+            + multi_s_y_deriv_1(istate,jstate)**2 &
+            + multi_s_z_deriv_1(istate,jstate)**2) 
+    enddo
+  enddo
+
+END_PROVIDER
+
diff --git a/src/mol_properties/multi_s_dipole_moment.irp.f b/src/mol_properties/multi_s_dipole_moment.irp.f
new file mode 100644
index 00000000..d5e62799
--- /dev/null
+++ b/src/mol_properties/multi_s_dipole_moment.irp.f
@@ -0,0 +1,93 @@
+! Providers for the dipole moments along x,y,z and the total dipole
+! moments.
+
+! The dipole moment along the x axis is:
+! \begin{align*}
+! \mu_x = < \Psi_m | \sum_i x_i + \sum_A Z_A R_A | \Psi_n >
+! \end{align*}
+! where $i$ is used for the electrons and $A$ for the nuclei.
+! $Z_A$ the charge of the nucleus $A$ and $R_A$ its position in the
+! space.
+
+! And it can be computed using the (transition, if n /= m) density
+! matrix as a expectation value
+! \begin{align*}
+! <\Psi_n|x| \Psi_m > = \sum_p \gamma_{pp}^{nm} < \phi_p | x | \phi_p >
+!      + \sum_{pq, p \neq q} \gamma_{pq}^{nm} < \phi_p | x | \phi_q > +  < \Psi_m | \sum_A Z_A R_A | \Psi_n >
+! \end{align*}
+
+
+
+BEGIN_PROVIDER [double precision, multi_s_dipole_moment, (N_states, N_states)]
+&BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment, (N_states, N_states)]
+&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment, (N_states, N_states)]
+&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment, (N_states, N_states)]
+
+  implicit none
+
+  BEGIN_DOC
+  ! Providers for :
+  ! <\Psi_m|\mu_x|\Psi_n>
+  ! <\Psi_m|\mu_y|\Psi_n>
+  ! <\Psi_m|\mu_z|\Psi_n>
+  ! ||\mu|| = \sqrt{\mu_x^2 + \mu_y^2 + \mu_z^2}
+  !
+  ! <\Psi_n|x| \Psi_m > = \sum_p \gamma_{pp}^{nm} \bra{\phi_p} x \ket{\phi_p} 
+  !   + \sum_{pq, p \neq q} \gamma_{pq}^{nm} \bra{\phi_p} x \ket{\phi_q}
+  ! \Psi: wf
+  ! n,m indexes for the states
+  ! p,q: general spatial MOs 
+  ! gamma^{nm}: density matrix \bra{\Psi^n} a^{\dagger}_a a_i \ket{\Psi^m}
+  END_DOC
+
+  integer          :: istate,jstate ! States
+  integer          :: i,j           ! general spatial MOs
+  double precision :: nuclei_part_x, nuclei_part_y, nuclei_part_z
+ 
+  multi_s_x_dipole_moment = 0.d0
+  multi_s_y_dipole_moment = 0.d0
+  multi_s_z_dipole_moment = 0.d0
+ 
+  do jstate = 1, N_states
+    do istate = 1, N_states
+ 
+      do i = 1, mo_num  
+        do j = 1, mo_num  
+          multi_s_x_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_x(j,i)  
+          multi_s_y_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_y(j,i) 
+          multi_s_z_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_z(j,i) 
+        enddo
+      enddo
+ 
+    enddo
+  enddo
+ 
+  ! Nuclei part
+  nuclei_part_x = 0.d0
+  nuclei_part_y = 0.d0
+  nuclei_part_z = 0.d0
+ 
+  do i = 1,nucl_num 
+    nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) 
+    nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) 
+    nuclei_part_z += nucl_charge(i) * nucl_coord(i,3) 
+  enddo
+ 
+  ! Only if istate = jstate, otherwise 0 by the orthogonality of the states
+  do istate = 1, N_states
+    multi_s_x_dipole_moment(istate,istate) += nuclei_part_x
+    multi_s_y_dipole_moment(istate,istate) += nuclei_part_y
+    multi_s_z_dipole_moment(istate,istate) += nuclei_part_z
+  enddo
+ 
+  ! d = <Psi|r|Psi>
+  do jstate = 1, N_states
+    do istate = 1, N_states
+      multi_s_dipole_moment(istate,jstate) = &
+        dsqrt(multi_s_x_dipole_moment(istate,jstate)**2 & 
+            + multi_s_y_dipole_moment(istate,jstate)**2 &
+            + multi_s_z_dipole_moment(istate,jstate)**2) 
+    enddo
+  enddo
+
+END_PROVIDER
diff --git a/src/mol_properties/print_mol_properties.irp.f b/src/mol_properties/print_mol_properties.irp.f
new file mode 100644
index 00000000..3753a3dd
--- /dev/null
+++ b/src/mol_properties/print_mol_properties.irp.f
@@ -0,0 +1,24 @@
+subroutine print_mol_properties()
+
+  implicit none
+
+  BEGIN_DOC
+  ! Run the propertie calculations
+  END_DOC
+
+  ! Electric dipole moment
+  if (calc_dipole_moment) then
+    call print_dipole_moment
+  endif
+
+  ! Transition electric dipole moment
+  if (calc_tr_dipole_moment .and. N_states > 1) then
+    call print_transition_dipole_moment
+  endif
+
+  ! Oscillator strength
+  if (calc_osc_str .and. N_states > 1) then
+    call print_oscillator_strength 
+  endif
+
+end
diff --git a/src/mol_properties/print_properties.irp.f b/src/mol_properties/print_properties.irp.f
new file mode 100644
index 00000000..4c0a9f38
--- /dev/null
+++ b/src/mol_properties/print_properties.irp.f
@@ -0,0 +1,194 @@
+! Dipole moments
+
+! Provided
+! | N_states                                   | integer          | Number of states                         |
+! | multi_s_x_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along x axis |
+! | multi_s_y_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along y axis |
+! | multi_s_z_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along z axis |
+! | multi_s_dipole_moment(N_states,N_states)   | double precision | Total (transition) dipole moments        |
+
+
+subroutine print_dipole_moment
+
+  implicit none
+
+  BEGIN_DOC
+  ! To print the dipole moment ||<\Psi_i|µ|\Psi_i>|| and its x,y,z components
+  END_DOC
+
+  integer :: istate
+  double precision, allocatable :: d(:), d_x(:), d_y(:), d_z(:)
+
+  allocate(d(N_states),d_x(N_states),d_y(N_states),d_z(N_states))
+ 
+  do istate = 1, N_states 
+    d_x(istate) = multi_s_x_dipole_moment(istate,istate)
+    d_y(istate) = multi_s_y_dipole_moment(istate,istate)
+    d_z(istate) = multi_s_z_dipole_moment(istate,istate)
+    d(istate)   = multi_s_dipole_moment(istate,istate)
+  enddo
+
+  ! Atomic units
+  print*,''
+  print*,'# Dipoles:'
+  print*,'=============================================='
+  print*,' Dipole moments (au)'
+  print*,' State      X           Y           Z         ||µ||' 
+
+  do istate = 1, N_states 
+    write(*,'(I5,4(F12.6))') (istate-1), d_x(istate), d_y(istate), d_z(istate), d(istate)
+  enddo
+
+  ! Debye
+  print*,''
+  print*,' Dipole moments (D)'
+  print*,' State      X           Y           Z         ||µ||' 
+
+  do istate = 1, N_states 
+    write(*,'(I5,4(F12.6))') (istate-1), d_x(istate)*au_to_D, d_y(istate)*au_to_D, d_z(istate)*au_to_D, d(istate)*au_to_D
+  enddo
+
+  print*,'=============================================='
+  print*,''
+
+  deallocate(d,d_x,d_y,d_z)
+
+ end
+
+! Transition dipole moments
+
+! Provided
+! | N_states                                   | integer          | Number of states                         |
+! | multi_s_x_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along x axis |
+! | multi_s_y_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along y axis |
+! | multi_s_z_dipole_moment(N_states,N_states) | double precision | (transition) dipole moments along z axis |
+! | multi_s_dipole_moment(N_states,N_states)   | double precision | Total (transition) dipole moments        |
+
+
+subroutine print_transition_dipole_moment
+
+  implicit none
+
+  BEGIN_DOC
+  ! To print the transition dipole moment ||<\Psi_i|µ|\Psi_j>|| and its components along x, y and z
+  END_DOC
+
+  integer          :: istate,jstate, n_states_print
+  double precision :: f, d, d_x, d_y, d_z, dip_str
+
+  if (N_states == 1 .or. N_det == 1) then
+    return
+  endif
+
+  print*,''
+  print*,'# Transition dipoles:'
+  print*,'=============================================='
+  print*,' Transition dipole moments (au)'
+  write(*,'(A89)') '   #  Transition       X           Y           Z         ||µ||     Dip. str.   Osc. str.'
+ 
+  if (print_all_transitions) then
+    n_states_print = N_states
+  else
+   n_states_print = 1
+  endif
+
+  do jstate = 1, n_states_print !N_states
+    do istate = jstate + 1, N_states
+      d_x = multi_s_x_dipole_moment(istate,jstate)
+      d_y = multi_s_y_dipole_moment(istate,jstate)
+      d_z = multi_s_z_dipole_moment(istate,jstate)
+      dip_str = d_x**2 + d_y**2 + d_z**2
+      d = multi_s_dipole_moment(istate,jstate)
+      f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
+      write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), '  ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
+    enddo
+  enddo
+
+  print*,''
+  print*,' Transition dipole moments (D)'
+  write(*,'(A89)') '   #  Transition       X           Y           Z         ||µ||     Dip. str.   Osc. str.'
+  
+  do jstate = 1, n_states_print !N_states
+    do istate = jstate + 1, N_states
+      d_x = multi_s_x_dipole_moment(istate,jstate) * au_to_D
+      d_y = multi_s_y_dipole_moment(istate,jstate) * au_to_D
+      d_z = multi_s_z_dipole_moment(istate,jstate) * au_to_D
+      d = multi_s_dipole_moment(istate,jstate)
+      dip_str = d_x**2 + d_y**2 + d_z**2
+      f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
+      d = multi_s_dipole_moment(istate,jstate) * au_to_D
+      write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), '  ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
+    enddo
+  enddo
+  print*,'=============================================='
+  print*,''
+
+end
+
+! Oscillator strengths
+
+! Provided
+! | N_states                                 | integer          | Number of states                  |
+! | multi_s_dipole_moment(N_states,N_states) | double precision | Total (transition) dipole moments |
+! | multi_s_deriv1_moment(N_states,N_states) | double precision | Total (transition) ...            |
+! | ci_energy_no_diag(N_states)              | double precision | CI energy of each state |
+
+! Internal
+! | f_l            | double precision | Oscillator strength in length gauge   |
+! | f_v            | double precision | Oscillator strength in velocity gauge |
+! | f_m            | double precision | Oscillator strength in mixed gauge    |
+! | n_states_print | integer          | Number of printed states              |
+
+
+subroutine print_oscillator_strength
+
+  implicit none
+
+  BEGIN_DOC
+  ! https://doi.org/10.1016/j.cplett.2004.03.126
+  ! Oscillator strength in:
+  ! - length gauge, f^l_{ij} = 2/3 (E_i - E_j) <\Psi_i|r|\Psi_j> <\Psi_j|r|\Psi_i>
+  ! - velocity gauge, f^v_{ij} = 2/3 (E_i - E_j)^(-1) <\Psi_i|v|\Psi_j> <\Psi_j|v|\Psi_i>
+  ! - mixed gauge, f^m_{ij} = -2i/3 <\Psi_i|r|\Psi_j> <\Psi_j|v|\Psi_i> 
+  END_DOC
+  
+  integer :: istate,jstate,k, n_states_print
+  double precision :: f_l,f_v,f_m,d,v
+
+  if (N_states == 1 .or. N_det == 1) then
+    return
+  endif
+
+  print*,''
+  print*,'# Oscillator strength:'
+  print*,'=============================================='
+
+  if (print_all_transitions) then
+    n_states_print = N_states
+  else
+   n_states_print = 1
+  endif
+
+  write(*,'(A103)') '  Oscillator strength in length gauge (f_l), velocity gauge (f_v) and mixed length-velocity gauge (f_m)'
+  do jstate = 1, n_states_print !N_states
+    do istate = jstate + 1, N_states
+      d = multi_s_dipole_moment(istate,jstate)
+      v = multi_s_deriv_1(istate,jstate)
+      ! Length gauge
+      f_l = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
+      ! Velocity gauge
+      f_v = 2d0/3d0 * v * v * 1d0/dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
+      ! Mixed gauge
+      f_m = 2d0/3d0 * d * v
+
+      write(*,'(A19,I3,A9,F10.6,A5,F7.1,A10,F9.6,A6,F9.6,A6,F9.6,A8,F7.3)') '   #  Transition n.', (istate-1), ': Excit.=', dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*ha_to_ev, &
+      ' eV ( ',dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*Ha_to_nm,' nm), f_l=',f_l, ', f_v=', f_v, ', f_m=', f_m, ', <S^2>=', s2_values(istate)
+      !write(*,'(I4,I4,A4,I3,A6,F6.1,A6,F6.1)') (istate-1), (jstate-1), '  ->', (istate-1), ', %T1=', percent_exc(2,istate), ', %T2=',percent_exc(3,istate)
+  
+    enddo
+  enddo
+
+  print*,'=============================================='
+  print*,''
+
+end
diff --git a/src/mol_properties/properties.irp.f b/src/mol_properties/properties.irp.f
new file mode 100644
index 00000000..c781c723
--- /dev/null
+++ b/src/mol_properties/properties.irp.f
@@ -0,0 +1,14 @@
+program mol_properties
+
+  implicit none
+
+  BEGIN_DOC
+  ! Run the propertie calculations
+  END_DOC
+
+  read_wf = .True.
+  touch read_wf
+
+  call print_mol_properties()
+
+end